From 9900b10532fd505b46a48951853ea5de27a91bb4 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sat, 9 May 2015 23:32:59 +0000 Subject: [PATCH] cohpsk freq tracking unit test, integrated into tcohpsk and doing sensible things with +/- 0.5Hz/s freq drift git-svn-id: https://svn.code.sf.net/p/freetel/code@2127 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/cohpsk.m | 51 ++++++++-- codec2-dev/octave/tcohpsk.m | 52 ++++++++-- codec2-dev/octave/test_ftrack.m | 168 ++++++++++++++++++++++++++++++++ 3 files changed, 255 insertions(+), 16 deletions(-) create mode 100644 codec2-dev/octave/test_ftrack.m diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index f5ccd433..4d06ba80 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -589,17 +589,56 @@ function [next_sync cohpsk] = coarse_freq_offset_est(cohpsk, fdmdv, ch_fdm_frame endfunction -function [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, nsymb, nin) +function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, f_est, nsymb, nin, freq_track) M = afdmdv.M; + % figure(10); clf; hold on; for r=1:nsymb + % downconvert entire signal to nominal baseband + + [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame(1+(r-1)*M:r*M), -f_est, afdmdv.Fs, afdmdv.fbb_phase_rx); + % downconvert each FDM carrier to Nc separate baseband signals - [rx_baseband afdmdv] = fdm_downconvert(afdmdv, rx_fdm_frame_bb(1+(r-1)*M:r*M), nin); + [rx_baseband afdmdv] = fdm_downconvert(afdmdv, rx_fdm_frame_bb, nin); [rx_filt afdmdv] = rx_filter(afdmdv, rx_baseband, nin); [rx_onesym rx_timing env afdmdv] = rx_est_timing(afdmdv, rx_filt, nin); ch_symb(r,:) = rx_onesym; + + % freq tracking, see test_ftrack.m for unit test. Placed in + % this function as it needs to work on a symbol by symbol basis + % rather than frame by frame. This means the control loop + % operates at a sample rate of Rs = 50Hz for say 1 Hz/s drift. + + if freq_track + beta = 0.005; + g = 0.2; + + % combine difference on phase from last symbol over Nc carriers + + mod_strip = 0; + for c=1:afdmdv.Nc+1 + adiff = rx_onesym(c) .* conj(afdmdv.prev_rx_symb(c)); + afdmdv.prev_rx_symb(c) = rx_onesym(c); + amod_strip = adiff.^4; + amod_strip = abs(real(amod_strip)) + j*imag(amod_strip); + mod_strip += amod_strip; + end + %plot(mod_strip) + + % 4th power strips QPSK modulation, by multiplying phase by 4 + % Using the abs value of the real coord was found to help + % non-linear issues when noise power was large + + + % loop filter made up of 1st order IIR plus integrator. Integerator + % was found to be reqd + + afdmdv.filt = ((1-beta)*afdmdv.filt + beta*angle(mod_strip)); + f_est += g*afdmdv.filt; + + end end endfunction @@ -683,7 +722,7 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne max_corr = 0; st = cohpsk.f_fine_est - 1; en = cohpsk.f_fine_est + 1; - for f_fine = st:0.25*en + for f_fine = st:0.25:en f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)'; corr = 0; mag = 0; for c=1:Nc*Nd @@ -764,10 +803,10 @@ function [sync cohpsk] = sync_state_machine(cohpsk, sync, next_sync) else cohpsk.sync_timer = 0; end - %printf(" ratio: %f sync timer: %d\n", cohpsk.ratio, cohpsk.sync_timer); + % printf(" ratio: %f sync timer: %d\n", cohpsk.ratio, cohpsk.sync_timer); - if cohpsk.sync_timer == 5 - printf(" lost sync ....\n"); + if cohpsk.sync_timer == 10 + printf(" [%d] lost sync ....\n", cohpsk.frame); next_sync = 0; end end diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 6fdc4f24..2a9cb7cd 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -22,9 +22,9 @@ randn('state',1); % test parameters ---------------------------------------------------------- -frames = 100; -foff = 0; -dfoff = 0; +frames = 200; +foff = -20; +dfoff = 0/Fs; EsNodB = 12; fading_en = 1; hf_delay_ms = 2; @@ -78,6 +78,9 @@ afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M; afdmdv.rx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter); +afdmdv.filt = 0; +afdmdv.prev_rx_symb = zeros(1,afdmdv.Nc+1); + % COHPSK Init -------------------------------------------------------- acohpsk = standard_init(); @@ -120,6 +123,8 @@ rx_fdm_frame_log = []; ct_symb_ff_log = []; rx_timing_log = []; ratio_log = []; +foff_log = []; +fest_log = []; % Channel modeling and BER measurement ---------------------------------------- @@ -177,7 +182,9 @@ for f=1:frames phase_ch *= foff_rect; ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch; end + foff_log = [foff_log foff]; phase_ch /= abs(phase_ch); + % printf("foff: %f ", foff); if fading_en ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M); @@ -221,8 +228,8 @@ for f=1:frames % we are out of sync so reset f_est and process two frames to clean out memories 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] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, Nsw*acohpsk.Nsymbrowpilot, nin); + %[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); end @@ -237,8 +244,8 @@ for f=1:frames 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); + %[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); end @@ -261,12 +268,17 @@ for f=1:frames % If in sync just do sample rate processing on latest frame if sync == 1 - [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame, -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, acohpsk.Nsymbrowpilot, nin); - acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb, acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot); + %[rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame, -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, acohpsk.Nsymbrowpilot, nin); + [ch_symb rx_timing rx_filt rx_baseband afdmdv acohpsk.f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, acohpsk.f_est, acohpsk.Nsymbrowpilot, nin, 1); + [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb, sync, next_sync); acohpsk.ct_symb_ff_buf(1:2,:) = acohpsk.ct_symb_ff_buf(acohpsk.Nsymbrowpilot+1:acohpsk.Nsymbrowpilot+2,:); acohpsk.ct_symb_ff_buf(3:acohpsk.Nsymbrowpilot+2,:) = acohpsk.ct_symb_buf(acohpsk.ct+3:acohpsk.ct+acohpsk.Nsymbrowpilot+2,:); + + ch_symb_log = [ch_symb_log; ch_symb]; + rx_timing_log = [rx_timing_log rx_timing]; + fest_log = [fest_log acohpsk.f_est]; end % if we are in sync complete demodulation with symbol rate processing @@ -373,12 +385,32 @@ else figure(1) plot(rx_symb_log*exp(j*pi/4),'+') + title('Scatter'); + figure(2) plot(rx_timing_log) + title('rx timing'); + figure(3) stem(nerr_log) + title('Bit Errors'); + figure(4) stem(ratio_log) + title('Sync ratio'); + + figure(5); + subplot(211) + plot(foff_log); + hold on; + plot(fest_log - Fcentre,'g'); + hold off; + + title('freq offset error'); + + subplot(212) + plot(foff_log(1:length(fest_log)) - fest_log + Fcentre) + title('freq offset estimation error'); end diff --git a/codec2-dev/octave/test_ftrack.m b/codec2-dev/octave/test_ftrack.m new file mode 100644 index 00000000..5b02e523 --- /dev/null +++ b/codec2-dev/octave/test_ftrack.m @@ -0,0 +1,168 @@ +% test_ftrack.m +% +% David Rowe May 2015 +% +% 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 + +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; + +phase_ch = 1; +prev_rx_symb = 1; + +EsNodB = 8; +EsNo = 10^(EsNodB/10); +variance = 1/EsNo; + +Nsym = 5; + +beta = 0.005; +g = 0.2; +ftrack = 0; +filt = 0; + +% Second order IIR filter coeffs, derived with pencil and paper + +a1 = (2-beta)/(1+g*beta); +a2 = (1-beta)/(1+g*beta); + +% Which can be used to work out the polar coordinates of the pole + +gamma = sqrt(a2); % radius (distance from origin) +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) + case (0) symbol = 1; + case (1) symbol = j; + case (2) symbol = -j; + case (3) symbol = -1; + 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 + +tx_bits = rand(1, Nbits) > 0.5; +tx_symb = zeros(1, Nsymb); +for i=1:Nsymb + tx_symb(i) = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); +end + +% run symbol by symbol + +rx_symb = zeros(1, Nsymb); +rx_filt = zeros(1,Nsym); + +diff_log = []; +foff_log = []; +ftrack_log = []; +mod_strip_log = []; + +for i=1:Nsymb + + % channel + + foff_rect = exp(j*2*pi*(foff-ftrack)*M/Fs); + foff += M*dfoff; + phase_ch *= foff_rect; + rx_symb(i) = tx_symb(i) * phase_ch; + noise = sqrt(variance*0.5)*(randn(1,1) + j*randn(1, 1)); + rx_symb(i) += noise; + + % simulate delay of root nyquist filter + + rx_filt(1:Nsym-1) = rx_filt(2:Nsym); + rx_filt(Nsym) = rx_symb(i); + + % freq tracking loop + + diff = rx_filt(1) .* conj(prev_rx_symb); + prev_rx_symb = rx_filt(1); + + % 4th power strips QPSK modulation, by multiplying phase by 4 + % Using the abs value of the real coord was found to help + % non-linear issues when noise power was large + + 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 + filt = ((1-beta)*filt + beta*angle(mod_strip)); + ftrack += g*filt; + + diff_log = [diff_log diff]; + foff_log = [foff_log foff]; + ftrack_log = [ftrack_log ftrack]; + mod_strip_log = [mod_strip_log mod_strip]; +end + +% plot results + +figure(1) +clf +plot(angle(mod_strip_log)) +title('mod stripping phase'); + +figure(2); +clf +plot(mod_strip_log,'+') +axis([-2 2 -2 2]) +title('mod stripping scatter'); + +figure(3) +clf +subplot(211) +plot(real(mod_strip_log)) +title('mod stripping real and imag'); +subplot(212) +plot(imag(mod_strip_log)) + +figure(4) +clf +subplot(211) +plot(foff_log,';foff;') +hold on +plot(ftrack_log,'g+;ftrack;') +hold off; +legend("boxoff"); +ylabel('Freq Offset Hz'); + +subplot(212) +plot(foff_log-ftrack_log,';foff-ftrack;' ) +ylabel('Freq Offset Tracking Error Hz'); +xlabel('symbols') +legend("boxoff"); + +figure(5) +freqz(g*beta/(1+g*beta), [1 -2*gamma*cos(w) gamma*gamma]) -- 2.25.1