From a813c65d3a0f3e33928d41a34d2a1d9ec0363296 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Fri, 23 Mar 2012 05:10:57 +0000 Subject: [PATCH] another pass at freq offset estimation, this time using DBPSK pilot git-svn-id: https://svn.code.sf.net/p/freetel/code@351 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fdmdv.m | 74 +++++++++++++++++++++--------------- codec2-dev/octave/fdmdv_ut.m | 21 +++++----- 2 files changed, 56 insertions(+), 39 deletions(-) diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m index 82d3bd37..d9a2a442 100644 --- a/codec2-dev/octave/fdmdv.m +++ b/codec2-dev/octave/fdmdv.m @@ -28,7 +28,7 @@ global Fsep = 75; % Separation between carriers (Hz) global Fcentre = 1200; % Centre frequency, Nc/2 below this, N/c above (Hz) global Nt = 5; % number of symbols we estimate timing over global P = 4; % oversample factor used for rx symbol filtering -global pilot_phase = pi/4;% current phase of symbol used to make pilot carrier +global pilot_bit = 0; % current phase of symbol used to make pilot carrier % Generate root raised cosine (Root Nyquist) filter --------------- @@ -73,7 +73,7 @@ global gt_alpha5_root = real((ifft_GF_alpha5_root(1:Nfilter))); function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation) global Nc; global Nb; - global pilot_phase; + global pilot_bit; % re-arrange as (Nc,Nb) matrix @@ -108,8 +108,16 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation) % +1 -1 +1 -1 BPSK sync carrier, once filtered becomes two spectral % lines at +/- Rs/2 - pilot_phase *= -1; - tx_symbols(c+1) = pilot_phase; + if pilot_bit + tx_symbols(Nc+1) = -prev_tx_symbols(c+1); + else + tx_symbols(Nc+1) = prev_tx_symbols(c+1); + end + if pilot_bit + pilot_bit = 0; + else + pilot_bit = 1; + end endfunction @@ -141,7 +149,7 @@ endfunction % Construct FDM signal by frequency shifting each filtered symbol % stream -function tx_fdm = fdm_upconvert(tx_filt) +function [tx_fdm pilot] = fdm_upconvert(tx_filt) global Fs; global M; global Nc; @@ -174,7 +182,8 @@ function tx_fdm = fdm_upconvert(tx_filt) c = Nc+1; for i=1:M phase_tx(c) = phase_tx(c) * freq(c); - tx_fdm(i) = tx_fdm(i) + 4*tx_filt(c,i)*phase_tx(c); + pilot(i) = sqrt(2)*2*tx_filt(c,i)*phase_tx(c); + tx_fdm(i) = tx_fdm(i) + pilot(i); end % Scale such that total Carrier power C of real(tx_fdm) = Nc. This @@ -215,6 +224,15 @@ function rx_baseband = fdm_downconvert(rx_fdm) rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)'; end end + + % Pilot + + c = Nc+1; + for i=1:M + phase_rx(c) = phase_rx(c) * freq(c); + rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)'; + end + endfunction @@ -229,7 +247,7 @@ function rx_filt = rx_filter(rx_baseband) global gt_alpha5_root; global Fsep; - rx_filt = zeros(Nc,P); + rx_filt = zeros(Nc+1,P); % rx filter each symbol, generate P filtered output samples for each symbol. % Note we keep memory at rate M, it's just the filter output at rate P @@ -245,10 +263,10 @@ function rx_filt = rx_filter(rx_baseband) endfunction -% Estimate frequency offset of FDM signal by peak picking pilot -% tone at Fcentre +% Estimate frequency offset of FDM signal using BPSK pilot. This is quite +% sensitive to pilot tone level wrt other carriers -function foff = rx_est_freq_offset(rx_fdm) +function foff = rx_est_freq_offset(rx_fdm, pilot) global M; global Nc; global Fs; @@ -265,20 +283,15 @@ function foff = rx_est_freq_offset(rx_fdm) global pilot_coeff; % down convert latest M samples of pilot by multiplying by - % ideal two-tone BPSK signal + % ideal two-tone BPSK pilot signal pilot_baseband(1:Npilotbaseband-M) = pilot_baseband(M+1:Npilotbaseband); c = Nc+1; for i=1:M - phase_rx_pilot(1) *= freq_rx_pilot(1); - phase_rx_pilot(2) *= freq_rx_pilot(2); - % phase_rx_pilot(1) *= freq(c); - pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i) * (phase_rx_pilot(1) + phase_rx_pilot(2)); - %pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i) * phase_rx_pilot(1); - %pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i); + pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i)*conj(pilot(i)); end - % LPF cutoff 200Hz + % LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf); j = 1; @@ -295,12 +308,11 @@ function foff = rx_est_freq_offset(rx_fdm) S = abs(fft(s, Mpilotfft)); figure(3) - clf - subplot(211) - plot(real(pilot_baseband)); - subplot(212) - plot(abs(fft(pilot_baseband))); - %plot(S); + %plot(real(pilot_baseband)) + plot(real(s)) + figure(4) + %plot(abs(fft(pilot_baseband))) + plot(S) % peak pick and convert to Hz @@ -321,6 +333,7 @@ endfunction function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband) global M; global Nt; + global Nc; global rx_filter_mem_timing; global rx_baseband_mem_timing; global P; @@ -331,7 +344,7 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband) % update buffer of Nt rate P filtered symbols rx_filter_mem_timing(:,1:(Nt-1)*P) = rx_filter_mem_timing(:,P+1:Nt*P); - rx_filter_mem_timing(:,(Nt-1)*P+1:Nt*P) = rx_filt; + rx_filter_mem_timing(:,(Nt-1)*P+1:Nt*P) = rx_filt(1:Nc,:); % sum envelopes of all carriers @@ -372,7 +385,7 @@ endfunction % Phase estimation function - probably won't work over a HF channel -% Tries to operate over a sinle symbol but uses phase information from +% Tries to operate over a single symbol but uses phase information from % all Nc carriers which should increase the SNR of phase estimate. % Maybe phase is coherent over a couple of symbols in HF channel,not % sure but it's worth 3dB so worth experimenting or using coherent as @@ -419,6 +432,7 @@ function rx_bits = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation) rx_bits(2*(c-1)+1) = msb; rx_bits(2*(c-1)+2) = lsb; end + else % map (Nc,1) QPSK symbols back into an (1,Nc*Nb) array of bits @@ -478,7 +492,7 @@ endfunction % Initialise ---------------------------------------------------- global tx_filter_memory = zeros(Nc+1, Nfilter); -global rx_filter_memory = zeros(Nc, Nfilter); +global rx_filter_memory = zeros(Nc+1, Nfilter); % phasors used for up and down converters @@ -501,17 +515,17 @@ global phase_rx = ones(Nc+1,1); global Npilotcoeff = 30; % number of pilot LPF coeffs global pilot_coeff = fir1(Npilotcoeff, 200/(Fs/2))'; % 200Hz LPF -global Npilotbaseband = Npilotcoeff + 4*M; % number of pilot baseband samples reqd for pilot LPF +global Npilotbaseband = Npilotcoeff + 4*M; % number of pilot baseband samples reqd for pilot LPF global Npilotlpf = 4*M; % number of samples we DFT pilot over, pilot est window global pilot_baseband = zeros(1, Npilotbaseband); % pilot baseband samples global pilot_lpf = zeros(1, Npilotlpf); % LPC pilot samples global phase_rx_pilot = [1 1]; -global freq_rx_pilot = [ exp(-j*2*pi*(Fcentre-Rs/2)/Fs) exp(-j*2*pi*(Fcentre+Rs/2)/Fs) ]; +global freq_rx_pilot = [ exp(-j*2*pi*(Fcentre-Rs/4)/Fs) exp(-j*2*pi*(Fcentre+Rs/4)/Fs) ]; % Timing estimator states global rx_filter_mem_timing = zeros(Nc, Nt*P); -global rx_baseband_mem_timing = zeros(Nc, Nfiltertiming); +global rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming); % Test bit stream state variables diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m index 4e8588d3..7f0f4215 100644 --- a/codec2-dev/octave/fdmdv_ut.m +++ b/codec2-dev/octave/fdmdv_ut.m @@ -15,13 +15,13 @@ fdmdv; % load modem code frames = 50; EbNo_dB = 7.3; -Foff_hz = 0; +Foff_hz = -100; modulation = 'dqpsk'; % ------------------------------------------------------------ tx_filt = zeros(Nc,M); -rx_symbols_log = zeros(Nc,1); +rx_symbols_log = []; rx_phase_log = 0; rx_timing_log = 0; tx_pwr = 0; @@ -31,8 +31,9 @@ total_bits = 0; rx_fdm_log = []; rx_baseband_log = []; rx_bits_offset = zeros(Nc*Nb*2); -prev_tx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4); -prev_rx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4); +prev_tx_symbols = sqrt(2)*ones(Nc+1,1)*exp(j*pi/4); +prev_tx_symbols(Nc+1) = 1; +prev_rx_symbols = sqrt(2)*ones(Nc+1,1)*exp(j*pi/4); foff_log = []; tx_baseband_log = []; @@ -84,7 +85,7 @@ for i=1:frames prev_tx_symbols = tx_symbols; tx_baseband = tx_filter(tx_symbols); tx_baseband_log = [tx_baseband_log tx_baseband]; - tx_fdm = fdm_upconvert(tx_baseband); + [tx_fdm pilot] = fdm_upconvert(tx_baseband); tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M); % ------------------- @@ -94,8 +95,10 @@ for i=1:frames % frequency offset for i=1:M - Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs)); - t++; + % Time varying freq offset + % Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs)); + % t++; + Foff = Foff_hz; freq_offset = exp(j*2*pi*Foff/Fs); phase_offset *= freq_offset; rx_fdm(i) = phase_offset*real(tx_fdm(i)); @@ -119,7 +122,7 @@ for i=1:frames % frequency offset estimation and correction - foff = rx_est_freq_offset(rx_fdm); + foff = rx_est_freq_offset(rx_fdm, pilot); foff_log = [ foff_log foff ]; %foff = 0; foff_rect = exp(j*2*pi*foff/Fs); @@ -168,7 +171,7 @@ printf("Eb/No (meas): %2.2f (%2.2f) dB %d bits %d errors QPSK BER (meas): %1. figure(1) clf; [n m] = size(rx_symbols_log); -plot(real(rx_symbols_log(:,20:m)),imag(rx_symbols_log(:,20:m)),'+') +plot(real(rx_symbols_log(1:Nc,20:m)),imag(rx_symbols_log(1:Nc,20:m)),'+') figure(2) clf; -- 2.25.1