From 3bd220da500f4695b841fc89048d17457c89eec2 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Wed, 21 Mar 2012 05:00:27 +0000 Subject: [PATCH] first pass at BPSK pilot frequency offset sestimation, works OK at Eb/No=7.3 operating point, tracks severl Hz/sec drift git-svn-id: https://svn.code.sf.net/p/freetel/code@350 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fdmdv.m | 60 ++++++++++++++++++++++++++---------- codec2-dev/octave/fdmdv_ut.m | 22 ++++++++++--- 2 files changed, 62 insertions(+), 20 deletions(-) diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m index 7f5927bd..82d3bd37 100644 --- a/codec2-dev/octave/fdmdv.m +++ b/codec2-dev/octave/fdmdv.m @@ -28,6 +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 % Generate root raised cosine (Root Nyquist) filter --------------- @@ -66,11 +67,13 @@ global gt_alpha5_root = real((ifft_GF_alpha5_root(1:Nfilter))); % Functions ---------------------------------------------------- -% generate Nc QPSK symbols from vector of (1,Nc*Nb) input bits +% generate Nc+1 QPSK symbols from vector of (1,Nc*Nb) input bits. The +% Nc+1 symbol is the +1 -1 +1 .... BPSK sync carrier function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation) global Nc; global Nb; + global pilot_phase; % re-arrange as (Nc,Nb) matrix @@ -102,6 +105,12 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation) tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2); endif + % +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; + endfunction @@ -115,7 +124,7 @@ function tx_baseband = tx_filter(tx_symbols) global Nfilter; global gt_alpha5_root; - tx_baseband = zeros(Nc,M); + tx_baseband = zeros(Nc+1,M); % tx filter each symbol, generate M filtered output samples for each symbol. % Efficient polyphase filter techniques used as tx_filter_memory is sparse @@ -125,7 +134,7 @@ function tx_baseband = tx_filter(tx_symbols) tx_baseband(:,i) = M*tx_filter_memory(:,M:M:Nfilter) * gt_alpha5_root(M-i+1:M:Nfilter)'; end tx_filter_memory(:,1:Nfilter-M) = tx_filter_memory(:,M+1:Nfilter); - tx_filter_memory(:,Nfilter-M+1:Nfilter) = zeros(Nc,M); + tx_filter_memory(:,Nfilter-M+1:Nfilter) = zeros(Nc+1,M); endfunction @@ -151,7 +160,7 @@ function tx_fdm = fdm_upconvert(tx_filt) end end - % Nc/2 tones above centre freq + % Nc/2 tones above centre freq for c=Nc/2+1:Nc for i=1:M @@ -160,16 +169,17 @@ function tx_fdm = fdm_upconvert(tx_filt) end end - % add centre pilot tone at twice amplitude of other tones + % add centre pilot tone c = Nc+1; for i=1:M phase_tx(c) = phase_tx(c) * freq(c); - tx_fdm(i) = tx_fdm(i) + 1/sqrt(2)*phase_tx(c); + tx_fdm(i) = tx_fdm(i) + 4*tx_filt(c,i)*phase_tx(c); end - % Scale such that total Carrier power C of real(tx_fdm) = Nc - % we return the complex (single sided) signal to make frequency + % Scale such that total Carrier power C of real(tx_fdm) = Nc. This + % excludes the power of the pilot tone. + % We return the complex (single sided) signal to make frequency % shifting for the purpose of testing easier tx_fdm = sqrt(2)*tx_fdm; @@ -242,8 +252,11 @@ function foff = rx_est_freq_offset(rx_fdm) global M; global Nc; global Fs; + global Rs; + global Fcentre; global freq; - global phase_rx; + global freq_rx_pilot; + global phase_rx_pilot; global Npilotbaseband; global Npilotlpf; global Npilotcoeff; @@ -251,13 +264,18 @@ function foff = rx_est_freq_offset(rx_fdm) global pilot_lpf; global pilot_coeff; - % down convert latest M samples of pilot - + % down convert latest M samples of pilot by multiplying by + % ideal two-tone BPSK signal + pilot_baseband(1:Npilotbaseband-M) = pilot_baseband(M+1:Npilotbaseband); c = Nc+1; for i=1:M - phase_rx(c) = phase_rx(c) * freq(c); - pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i) * phase_rx(c)'; + 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); end % LPF cutoff 200Hz @@ -276,11 +294,19 @@ function foff = rx_est_freq_offset(rx_fdm) s = pilot_lpf(1:Mpilot:Npilotlpf) .* hanning(Npilotlpf/Mpilot)'; S = abs(fft(s, Mpilotfft)); + figure(3) + clf + subplot(211) + plot(real(pilot_baseband)); + subplot(212) + plot(abs(fft(pilot_baseband))); + %plot(S); + % peak pick and convert to Hz [x ix] = max(S); r = 2*200/Mpilotfft; % maps FFT bin to frequency in Hz - + if ix > Mpilotfft/2 foff = (ix - Mpilotfft - 1)*r; else @@ -451,7 +477,7 @@ endfunction % Initialise ---------------------------------------------------- -global tx_filter_memory = zeros(Nc, Nfilter); +global tx_filter_memory = zeros(Nc+1, Nfilter); global rx_filter_memory = zeros(Nc, Nfilter); % phasors used for up and down converters @@ -475,10 +501,12 @@ 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 + 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) ]; % Timing estimator states diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m index 96652504..4e8588d3 100644 --- a/codec2-dev/octave/fdmdv_ut.m +++ b/codec2-dev/octave/fdmdv_ut.m @@ -13,9 +13,9 @@ fdmdv; % load modem code % Simulation Parameters -------------------------------------- -frames = 100; +frames = 50; EbNo_dB = 7.3; -Foff_hz = 124; +Foff_hz = 0; modulation = 'dqpsk'; % ------------------------------------------------------------ @@ -34,6 +34,7 @@ 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); foff_log = []; +tx_baseband_log = []; Ndelay = M+20; rx_fdm_delay = zeros(Ndelay,1); @@ -68,6 +69,7 @@ SNR = CNo_dB - 10*log10(B); phase_offset = 1; freq_offset = exp(j*2*pi*Foff_hz/Fs); foff_phase = 1; +t = 0; % Main loop ---------------------------------------------------- @@ -81,7 +83,8 @@ for i=1:frames tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation); prev_tx_symbols = tx_symbols; tx_baseband = tx_filter(tx_symbols); - tx_fdm = fdm_upconvert(tx_baseband) + 2*cos(2*pi*(0:M-1)*Fcentre/Fs); + tx_baseband_log = [tx_baseband_log tx_baseband]; + tx_fdm = fdm_upconvert(tx_baseband); tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M); % ------------------- @@ -91,6 +94,9 @@ for i=1:frames % frequency offset for i=1:M + Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs)); + t++; + freq_offset = exp(j*2*pi*Foff/Fs); phase_offset *= freq_offset; rx_fdm(i) = phase_offset*real(tx_fdm(i)); end @@ -115,6 +121,7 @@ for i=1:frames foff = rx_est_freq_offset(rx_fdm); foff_log = [ foff_log foff ]; + %foff = 0; foff_rect = exp(j*2*pi*foff/Fs); for i=1:M @@ -167,13 +174,18 @@ figure(2) clf; subplot(211) plot(rx_timing_log) +title('timing offset (samples)'); subplot(212) +plot(foff_log) +title('Freq offset (Hz)'); + +%figure(3) +%clf; %Nfft=Fs; %S=fft(rx_fdm_log,Nfft); %SdB=20*log10(abs(S)); %plot(-Fs/2+1:Fs/2,fftshift(SdB)) %plot(SdB(1:Fs/4)) -plot(foff_log) @@ -192,6 +204,8 @@ plot(foff_log) % dump file type plotting & instrumentation % determine if error pattern is bursty % HF channel simulation +% Offset or pi/4 QPSK and tests with real tx HPA +% real time SNR get function % % phase estimator not working too well and would need a UW % to resolve ambiguity. But this is probably worth it for -- 2.25.1