From: drowe67 Date: Mon, 11 May 2015 07:18:26 +0000 (+0000) Subject: extended freq offset estimation to +/-60Hz which is I think acceptable in the real... X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=e3362cda7f06329f64d317ee05a1440def020b91;p=freetel-svn-tracking.git extended freq offset estimation to +/-60Hz which is I think acceptable in the real world. gd results at AWGN and fading ch set points git-svn-id: https://svn.code.sf.net/p/freetel/code@2128 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index 4d06ba80..e742f652 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -710,7 +710,7 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne next_sync = 1; else next_sync = 0; - printf(" [%d] back to coarse freq offset est...\n", cohpsk.frame) ; + %printf(" [%d] back to coarse freq offset est...\n", cohpsk.frame) ; end cohpsk.ratio = abs(max_corr/max_mag); end diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 2a9cb7cd..f313a93a 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -22,9 +22,13 @@ randn('state',1); % test parameters ---------------------------------------------------------- -frames = 200; -foff = -20; -dfoff = 0/Fs; +% TODO +% [ ] set up various tests we use to characterise modem for easy +% repeating when we change modem + +frames = 100; +foff = -40; +dfoff = -0.5/Fs; EsNodB = 12; fading_en = 1; hf_delay_ms = 2; @@ -225,26 +229,40 @@ for f=1:frames if (sync == 0) && (f>1) - % we are out of sync so reset f_est and process two frames to clean out memories + % we can test +/- 20Hz, so we break this up into 3 tests to cover +/- 60Hz + + max_ratio = 0; + for acohpsk.f_est = Fcentre-40:40:Fcentre+40 + + printf(" [%d] acohpsk.f_est: %f +/- 20\n", f, acohpsk.f_est); - acohpsk.f_est = Fcentre; - %[rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame_buf, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx); - [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame_buf, acohpsk.f_est, Nsw*acohpsk.Nsymbrowpilot, nin, 0); - for i=1:Nsw-1 - acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb((i-1)*acohpsk.Nsymbrowpilot+1:i*acohpsk.Nsymbrowpilot,:), acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot); + % we are out of sync so reset f_est and process two frames to clean out memories + + [ch_symb rx_timing rx_filt rx_baseband afdmdv acohpsk.f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame_buf, acohpsk.f_est, Nsw*acohpsk.Nsymbrowpilot, nin, 0); + for i=1:Nsw-1 + acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb((i-1)*acohpsk.Nsymbrowpilot+1:i*acohpsk.Nsymbrowpilot,:), acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot); + end + [anext_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:), sync, next_sync); + + if anext_sync == 1 + %printf(" [%d] acohpsk.ratio: %f\n", f, acohpsk.ratio); + if acohpsk.ratio > max_ratio + max_ratio = acohpsk.ratio; + f_est = acohpsk.f_est - acohpsk.f_fine_est; + next_sync = anext_sync; + end + end end - [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:), sync, next_sync); if next_sync == 1 % we've found a sync candidate! % re-process last two frames with adjusted f_est then check again - acohpsk.f_est -= acohpsk.f_fine_est; + acohpsk.f_est = f_est; printf(" [%d] trying sync and f_est: %f\n", f, acohpsk.f_est); - %[rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame_buf, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx); [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame_buf, acohpsk.f_est, Nsw*acohpsk.Nsymbrowpilot, nin, 0); for i=1:Nsw-1 acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb((i-1)*acohpsk.Nsymbrowpilot+1:i*acohpsk.Nsymbrowpilot,:), acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot); @@ -384,29 +402,34 @@ else % some other useful plots figure(1) + clf; plot(rx_symb_log*exp(j*pi/4),'+') title('Scatter'); figure(2) + clf; plot(rx_timing_log) title('rx timing'); figure(3) + clf; stem(nerr_log) title('Bit Errors'); figure(4) + clf; stem(ratio_log) title('Sync ratio'); figure(5); + clf; subplot(211) - plot(foff_log); + plot(foff_log,';freq offset;'); hold on; - plot(fest_log - Fcentre,'g'); + plot(fest_log - Fcentre,'g;freq offset est;'); hold off; - - title('freq offset error'); + title('freq offset'); + legend("boxoff"); subplot(212) plot(foff_log(1:length(fest_log)) - fest_log + Fcentre) diff --git a/codec2-dev/octave/test_foff.m b/codec2-dev/octave/test_foff.m index 6958090b..6a424ee1 100644 --- a/codec2-dev/octave/test_foff.m +++ b/codec2-dev/octave/test_foff.m @@ -174,20 +174,35 @@ function sim_out = freq_off_est_test(sim_in) next_sync = sync; sync = 0; + if sync == 0 - acohpsk.f_est = Fcentre; - end - [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame_buf, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx); - [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, Nsw*acohpsk.Nsymbrowpilot, nin); - ch_symb_log = [ch_symb_log; ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:)]; + % we can test +/- 20Hz, so we break this up into 3 tests to cover +/- 60Hz - % coarse timing (frame sync) and initial fine freq est --------------------------------------------- + max_ratio = 0; + for acohpsk.f_est = Fcentre-40:40:Fcentre+40 + + printf(" [%d] acohpsk.f_est: %f +/- 20\n", f, acohpsk.f_est); + [ch_symb rx_timing rx_filt rx_baseband afdmdv acohpsk.f_est ] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame_buf, acohpsk.f_est, Nsw*acohpsk.Nsymbrowpilot, nin, 0); + ch_symb_log = [ch_symb_log; ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:)]; + + % coarse timing (frame sync) and initial fine freq est --------------------------------------------- - for i=1:Nsw-1 - acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb((i-1)*acohpsk.Nsymbrowpilot+1:i*acohpsk.Nsymbrowpilot,:), acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot); + for i=1:Nsw-1 + acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb((i-1)*acohpsk.Nsymbrowpilot+1:i*acohpsk.Nsymbrowpilot,:), acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot); + end + [anext_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:), sync, next_sync); + + if anext_sync == 1 + %printf(" [%d] acohpsk.ratio: %f\n", f, acohpsk.ratio); + if acohpsk.ratio > max_ratio + max_ratio = acohpsk.ratio; + f_est = acohpsk.f_est - acohpsk.f_fine_est; + next_sync = anext_sync; + end + end + end end - [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:), sync, next_sync); % we've found a sync candidate @@ -195,10 +210,9 @@ function sim_out = freq_off_est_test(sim_in) % rewind and re-process last few frames with f_est - acohpsk.f_est -= acohpsk.f_fine_est; + acohpsk.f_est = f_est; printf(" [%d] trying sync and f_est: %f\n", f, acohpsk.f_est); - [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame_buf, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx); - [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, Nsw*acohpsk.Nsymbrowpilot, nin); + [ch_symb rx_timing rx_filt rx_baseband afdmdv acohpsk.f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame_buf, acohpsk.f_est, Nsw*acohpsk.Nsymbrowpilot, nin, 0); for i=1:Nsw-1 acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb((i-1)*acohpsk.Nsymbrowpilot+1:i*acohpsk.Nsymbrowpilot,:), acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot); end @@ -236,7 +250,7 @@ endfunction function freq_off_est_test_single sim_in.frames = 100; sim_in.EsNodB = 12; - sim_in.foff = -20; + sim_in.foff = -59; sim_in.dfoff = 0; sim_in.fading_en = 1; @@ -268,7 +282,7 @@ function [freq_off_log EsNodBSet] = freq_off_est_test_curves EsNodBSet = [20 12 8]; sim_in.frames = 100; - sim_in.foff = -20; + sim_in.foff = -40; sim_in.dfoff = 0; sim_in.fading_en = 1; freq_off_log = 1E6*ones(sim_in.frames, length(EsNodBSet) ); @@ -354,8 +368,8 @@ function [freq_off_log EsNodBSet] = freq_off_est_test_curves endfunction -%freq_off_est_test_single; -freq_off_est_test_curves; +freq_off_est_test_single; +%freq_off_est_test_curves; % 1. start with +/- 20Hz offset % 2. Measure frames to sync. How to define sync? Foff to withn 1 Hz. Sync state diff --git a/codec2-dev/octave/test_ftrack.m b/codec2-dev/octave/test_ftrack.m index 5b02e523..5082c443 100644 --- a/codec2-dev/octave/test_ftrack.m +++ b/codec2-dev/octave/test_ftrack.m @@ -5,38 +5,75 @@ % Octave script that to test and develop symbol rate frequency offset % tracking % -% [ ] add Es/No noise and measure BER -% [ ] QPSK mod of random symbols -% [ ] tracking 1 Hz/s freq drift with negl impl loss -% [ ] simulate delay through filter, or implement R/N filter +% [X] add Es/No noise and measure BER +% [X] QPSK mod of random symbols +% [X] tracking 1 Hz/s freq drift shiel staying in lock +% [X] simulate R/N filter delay through filter rand('state',1); randn('state',1); graphics_toolkit ("gnuplot"); -Fs = 8000; -Rs = 50; -Nbits = 5000; -Nsymb = Nbits/2; -foff = 1; -dfoff = -1/Fs; -M = Fs/Rs; +Fs = 8000; % sample rate +Rs = 50; % QPSK modem symbol rate +Nbits = 5000; % number of bits to simulute over +Nsymb = Nbits/2; +foff = 1; % initial freq offset +dfoff = -1/Fs; % rate of change of frequency offset in Hz/sample +M = Fs/Rs; % oversampling rate -phase_ch = 1; -prev_rx_symb = 1; - -EsNodB = 8; +EsNodB = 8; % noise level EsNo = 10^(EsNodB/10); variance = 1/EsNo; -Nsym = 5; +Nsym = 5; % simulated delay of root nyquist filter in symbols + +% control system constants, these can be tweaked beta = 0.005; -g = 0.2; +g = 0.2; + +% init state variables + +phase_ch = 1; +prev_rx_symb = 1; ftrack = 0; filt = 0; -% Second order IIR filter coeffs, derived with pencil and paper +% QPSK symbols ------> x -------> x ----------G(z) ------------\ +% | | | +% exp(j*w_off*n) | | +% | | +% exp(-j*w_track*n) | +% \____________________________/ +% +% for small phase/freq errors, where phase/freq detector is linear: +% +% -1 beta 1 1 +% G(z) = 1 - z g --------------- -------- ------- +% -1 -1 -1 +% 1 - (1-beta)z 1 - z 1 - z +% +% (1) (2) (3) (4) +% +% (1) - differentiator that gives us an error based on frequency difference +% (2) - first order IIR loop filter +% (3) - integrator I added so loop could track ramp inputs (i.e. freq drift) +% (4) - frequency to phase integrator implicit in VCO +% +% Note (1) and (3) cancel out, H(z) = 1, so we get a 2nd order loop. We need to keep (1) +% in the time domain implementation to make the modulation stripper work, it likes angles +% near 0. +% +% Transfer function, derived with pencil and paper: +% +% G(z) g*beta +% ------------ = ---------- * ------------------------------------ +% 1 + G(z)H(z) 1 + g*beta 2 - beta -1 1 - beta -2 +% 1 - ---------- z + ----------- z +% 1 + g*beta 1 + g*beta + +% Second order IIR filter coeffs a1 = (2-beta)/(1+g*beta); a2 = (1-beta)/(1+g*beta); @@ -48,6 +85,7 @@ w = acos(a1/(2*gamma)); % angular frequency printf("2nd order loop w: %f gamma: %f\n", w, gamma); + function symbol = qpsk_mod(two_bits) two_bits_decimal = sum(two_bits .* [2 1]); switch(two_bits_decimal) @@ -58,17 +96,6 @@ function symbol = qpsk_mod(two_bits) endswitch endfunction -% Gray coded QPSK demodulation function - -function two_bits = qpsk_demod(symbol) - if isscalar(symbol) == 0 - printf("only works with scalars\n"); - return; - end - bit0 = real(symbol*exp(j*pi/4)) < 0; - bit1 = imag(symbol*exp(j*pi/4)) < 0; - two_bits = [bit1 bit0]; -endfunction % modulator @@ -78,7 +105,7 @@ for i=1:Nsymb tx_symb(i) = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); end -% run symbol by symbol +% run rest of simulation symbol by symbol rx_symb = zeros(1, Nsymb); rx_filt = zeros(1,Nsym); @@ -116,8 +143,10 @@ for i=1:Nsymb mod_strip = diff.^4; mod_strip = abs(real(mod_strip)) + j*imag(mod_strip); - % loop filter made up of 1st order IIR plus integrator. Integerator - % was found to be reqd + % loop filter made up of 1st order IIR plus integrator plus gain + % term. Integrator was found to be reqd to track ramp inputs + % (i.e. freq drift) + filt = ((1-beta)*filt + beta*angle(mod_strip)); ftrack += g*filt;