From 588f84926f9d6846d7ab3c945793e4c54c973504 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sat, 7 Apr 2018 20:38:06 +0000 Subject: [PATCH] characterising freq offset estimation metrics git-svn-id: https://svn.code.sf.net/p/freetel/code@3447 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/ofdm_dev.m | 138 ++++++++++++++++++++++++----------- codec2-dev/octave/ofdm_lib.m | 51 +++++++++---- 2 files changed, 132 insertions(+), 57 deletions(-) diff --git a/codec2-dev/octave/ofdm_dev.m b/codec2-dev/octave/ofdm_dev.m index d39bbef8..40f17ac6 100644 --- a/codec2-dev/octave/ofdm_dev.m +++ b/codec2-dev/octave/ofdm_dev.m @@ -424,26 +424,28 @@ function [sim_out rx states] = run_sim(sim_in) % print results of this simulation point to the console - if ldpc_en || diversity_en - printf("Coded EbNodB: % -4.1f BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpackets: %5d Tpacket_errs: %5d\n", - EbNodB(nn), Terrs_coded/(Nbits*rate), Nbits*rate, Terrs_coded, - Tpacketerrs_coded/Tpackets_coded, Tpackets_coded, Tpacketerrs_coded); - end - EbNodB_raw = EbNodB(nn) + 10*log10(rate); - printf("Raw EbNodB..: % -4.1f BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpackets: %5d Tpacket_errs: %5d\n", - EbNodB_raw, Terrs/Nbits, Nbits, Terrs, - Tpacketerrs/Tpackets, Tpackets, Tpacketerrs); - - % returns results for plotting curves - - if ldpc_en || diversity_en - sim_out.ber(nn) = Terrs_coded/(Nbits*rate); - sim_out.per(nn) = Tpacketerrs_coded/Tpackets_coded; - else - sim_out.ber(nn) = Terrs/Nbits; - sim_out.per(nn) = Tpacketerrs/Tpackets; - end + if verbose + if ldpc_en || diversity_en + printf("Coded EbNodB: % -4.1f BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpackets: %5d Tpacket_errs: %5d\n", + EbNodB(nn), Terrs_coded/(Nbits*rate), Nbits*rate, Terrs_coded, + Tpacketerrs_coded/Tpackets_coded, Tpackets_coded, Tpacketerrs_coded); + end + EbNodB_raw = EbNodB(nn) + 10*log10(rate); + printf("Raw EbNodB..: % -4.1f BER: %5.4f Tbits: %5d Terrs: %5d PER: %5.4f Tpackets: %5d Tpacket_errs: %5d\n", + EbNodB_raw, Terrs/Nbits, Nbits, Terrs, + Tpacketerrs/Tpackets, Tpackets, Tpacketerrs); + + % returns results for plotting curves + if ldpc_en || diversity_en + sim_out.ber(nn) = Terrs_coded/(Nbits*rate); + sim_out.per(nn) = Tpacketerrs_coded/Tpackets_coded; + else + sim_out.ber(nn) = Terrs/Nbits; + sim_out.per(nn) = Tpacketerrs/Tpackets; + end + end + % Optional plots, mostly used with run-single if verbose @@ -862,7 +864,7 @@ function [delta_ct delta_foff timing_mx_log] = acquisition_test(Ntests=10, EbNod [sim_out rx states] = run_sim(sim_in); - states.verbose = 0; + states.verbose = 2; % set up acquistion @@ -899,9 +901,9 @@ function [delta_ct delta_foff timing_mx_log] = acquisition_test(Ntests=10, EbNod ct_target = Nsamperframe/2; % actual known position of correct coarse timing for w=1:Nsamperframe:length(rx)-4*Nsamperframe - [ct_est foff_est timing_valid timing_mx1 timing_mx2] = coarse_sync(states, rx(w+st:w+en), rate_fs_pilot_samples); + [ct_est foff_est timing_valid timing_mx1 timing_mx2] = coarse_sync(states, real(rx(w+st:w+en)), rate_fs_pilot_samples); if states.verbose - printf("w: %d ct_est: %4d foff_est: %3.1f\n", w, ct_est, foff_est); + printf("w: %5d ct_est: %4d foff_est: %5.1f timing_mx1: %3.2f timing_mx2: %3.2f\n", w, ct_est, foff_est, timing_mx1, timing_mx2); end % valid coarse timing ests are modulo Nsamperframe @@ -989,37 +991,87 @@ function acquisition_histograms(fine_en = 0, foff) endfunction -% Used to develop sync state machine - in particular metric to show we -% are out of sync of have lost nodem signal +% Used to develop sync state machine - in particular a metric to show +% we are out of sync, or have sync with a bad freq offset est, or have +% lost nodem signal -function sync_metrics() +function sync_metrics(x_axis = 'EbNo') Fs = 8000; - Ntests = 10; - f_offHz = [0 0.5 1 2 5 10]; - EbNodB = [0 2 4 6 10 20]; + Ntests = 4; + f_offHz = [-25:25]; + EbNodB = [-10 4 20]; - figure(1); clf; - + mean_mx1_log = mean_mx2_log = mean_dfoff_log = []; for f = 1:length(f_offHz) af_offHz = f_offHz(f); - mean_mx1_log = mean_mx2_log = []; + mean_mx1_row = mean_mx2_row = mean_dfoff_row = []; for e = 1:length(EbNodB) aEbNodB = EbNodB(e); [dct dfoff timing_mx_log] = acquisition_test(Ntests, aEbNodB, af_offHz); mean_mx1 = mean(timing_mx_log(:,1)); mean_mx2 = mean(timing_mx_log(:,2)); - printf("f_offHz: %3.2f EbNodB: %3.2f mx1: %3.2f mx2: %3.2f\n", af_offHz, aEbNodB, mean_mx1, mean_mx2); - mean_mx1_log = [mean_mx1_log mean_mx1]; mean_mx2_log = [mean_mx2_log mean_mx2]; + printf("f_offHz: %5.2f EbNodB: % 6.2f mx1: %3.2f mx2: %3.2f\n", af_offHz, aEbNodB, mean_mx1, mean_mx2); + mean_mx1_row = [mean_mx1_row mean_mx1]; mean_mx2_row = [mean_mx2_row mean_mx2]; + mean_dfoff_row = [mean_dfoff_row mean(dfoff)]; end - if f == 2, hold on, end; - leg1 = sprintf("b+-;mx1 f_offHz %3.2f;", af_offHz); - leg2 = sprintf("g*-;mx2 f_offHz %3.2f;", af_offHz); - plot(EbNodB, mean_mx1_log, leg1) - plot(EbNodB, mean_mx2_log, leg2) + mean_mx1_log = [mean_mx1_log; mean_mx1_row]; mean_mx2_log = [mean_mx2_log; mean_mx2_row]; + mean_dfoff_log = [mean_dfoff_log; mean_dfoff_row]; end - hold off; - xlabel('Eb/No (dB'); - ylabel('Coefficient') + + figure(1); clf; + if strcmp(x_axis,'EbNo') + for f = 1:length(f_offHz) + if f == 2, hold on, end; + leg1 = sprintf("b+-;mx1 %4.1f Hz;", f_offHz(f)); + plot(EbNodB, mean_mx1_log(f,:), leg1) + end + for f = 1:length(f_offHz) + leg2 = sprintf("g*-;mx2 %4.1f Hz;", f_offHz(f)); + plot(EbNodB, mean_mx2_log(f,:), leg2) + end + hold off; + xlabel('Eb/No (dB)'); + ylabel('Coefficient') + title('Pilot Correlation Metric 1 and 2 against Eb/No for different Freq Offsets'); + legend("location", "northwest"); legend("boxoff"); + axis([min(EbNodB) max(EbNodB) 0 1.2]) + print('-deps', '-color', "ofdm_dev_pilot_correlation_ebno.eps") + end + + if strcmp(x_axis,'freq') + % x axis is freq + for e = 1:length(EbNodB) + if e == 2, hold on, end; + leg1 = sprintf("b+-;mx1 %3.0f dB;", EbNodB(e)); + plot(f_offHz, mean_mx1_log(:,e), leg1) + end + for e = 1:length(EbNodB) + leg2 = sprintf("g*-;mx2 %3.0f dB;", EbNodB(e)); + plot(f_offHz, mean_mx2_log(:,e), leg2) + end + hold off; + xlabel('freq offset (Hz)'); + ylabel('Coefficient') + title('Pilot Correlation Metric 1 and 2 against Freq Offset for different Eb/No dB'); + legend("location", "northwest"); legend("boxoff"); + axis([min(f_offHz) max(f_offHz) 0 1.2]) + print('-deps', '-color', "ofdm_dev_pilot_correlation_freq.eps") + + mean_dfoff_log + + figure(2); + for e = 1:length(EbNodB) + if e == 2, hold on, end; + leg1 = sprintf("+-;mx1 %3.0f dB;", EbNodB(e)); + plot(f_offHz, mean_dfoff_log(:,e), leg1) + end + hold off; + xlabel('freq offset (Hz)'); + ylabel('Mean Freq Est Error') + title('Freq Est Error against Freq Offset for different Eb/No dB'); + axis([min(f_offHz) max(f_offHz) -5 5]) + end + endfunction @@ -1036,5 +1088,5 @@ more off; %run_curves %run_curves_estimators %acquisition_histograms(0, 0) -%acquisition_test(10, 4, 5) -sync_metrics +acquisition_test(Ntests=3, EbNodB=100, foff_hz=0) +%sync_metrics('freq') diff --git a/codec2-dev/octave/ofdm_lib.m b/codec2-dev/octave/ofdm_lib.m index 5495a465..b41b93c3 100644 --- a/codec2-dev/octave/ofdm_lib.m +++ b/codec2-dev/octave/ofdm_lib.m @@ -70,10 +70,10 @@ function [t_est foff_est timing_valid timing_mx1 timing_mx2] = coarse_sync(state end [timing_mx1 t_est] = max(corr1); - timing_mx2 = max(corr2); - timing_valid = timing_mx1 > timing_mx_thresh; - - if verbose > 1 + [timing_mx2 t_est] = max(corr2); + timing_valid = timing_mx2 > timing_mx_thresh; + + if verbose > 2 printf(" mx1: %f mx2: %f timing_est: %d timing_valid: %d\n", timing_mx1, timing_mx2, timing_est, timing_valid); end @@ -93,10 +93,9 @@ function [t_est foff_est timing_valid timing_mx1 timing_mx2] = coarse_sync(state foff_est = foff_est_pos - 1; else foff_est = foff_est_neg - fmax - 1; - end - + end #} - + p1 = rx(t_est:t_est+Npsam/2-1) * rate_fs_pilot_samples(1:Npsam/2)'; p2 = rx(t_est+Npsam/2:t_est+Npsam-1) * rate_fs_pilot_samples(Npsam/2+1:Npsam)'; p3 = rx(t_est+Nsamperframe:t_est+Nsamperframe+Npsam/2-1) * rate_fs_pilot_samples(1:Npsam/2)'; @@ -104,22 +103,44 @@ function [t_est foff_est timing_valid timing_mx1 timing_mx2] = coarse_sync(state Fs1 = Fs/(Npsam/2); foff_est = Fs1*angle( (conj(p1)*p2 + conj(p3)*p4))/(2*pi); + if verbose > 1 - %printf("t_est: %d\n", t_est); figure(7); clf; - plot(abs(corr)) + plot(abs(corr2)) + + figure(8) + subplot(211) + plot(real(rate_fs_pilot_samples)) + hold on; plot(real(rx(t_est:t_est+Npsam-1)),'g'); hold off + subplot(212) + plot(imag(rate_fs_pilot_samples)) + hold on; plot(imag(rx(t_est:t_est+Npsam-1)),'g'); hold off + + figure(9) + dp = rx(t_est:t_est+Npsam-1) .* conj(rate_fs_pilot_samples); + subplot(211); plot(real(dp)); + subplot(212); plot(imag(dp)); + + figure(10) + plot(dp,'+') + + figure(11) + plot([p1 p2; p3 p4], '+') + %axis([1 Ncorr 0 2]) #{ figure(8) plot(C) axis([0 200 0 0.4]) #} + #{ figure(9) %rate_fs_pilot_samples plot([p1 p2 p3 p4] ,'b+') %plot(rate_fs_pilot_samples * rate_fs_pilot_samples' ,'b+') axis([-0.2 0.2 -0.2 0.2]) %hold on; plot(rx(t_est+Nsamperframe:t_est+Npsam+Nsamperframe-1) .* rate_fs_pilot_samples','g+'); hold off; + #} end endfunction @@ -186,7 +207,7 @@ function states = ofdm_init(bps, Rs, Tcp, Ns, Nc) states.foff_est_en = 1; states.phase_est_en = 1; - states.foff_est_gain = 0.01; + states.foff_est_gain = 0.05; states.foff_est_hz = 0; states.sample_point = states.timing_est = 1; states.nin = states.Nsamperframe; @@ -300,8 +321,8 @@ function [timing_valid states] = ofdm_sync_search(states, rxbuf_in) st = M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe; [ct_est foff_est timing_valid timing_mx1 timing_mx2] = coarse_sync(states, states.rxbuf(st:en), states.rate_fs_pilot_samples); - if states.verbose - printf(" ct_est: %4d foff_est: %3.1f timing_valid: %d timing_mx1: %f timing_mx2: %f\\n", ct_est, foff_est, timing_valid, timing_mx1, timing_mx2); + if verbose + printf(" mx1: %3.2f mx2: %3.2f coarse_foff: %4.1f\n", timing_mx1, timing_mx2, foff_est); end if timing_valid @@ -318,6 +339,7 @@ function [timing_valid states] = ofdm_sync_search(states, rxbuf_in) else states.nin = Nsamperframe; end + states.timing_valid = timing_valid; states.timing_mx1 = timing_mx1; states.timing_mx2 = timing_mx2; @@ -379,7 +401,7 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, end if verbose > 1 - printf(" ft_est: %2d timing_est: %2d mx1: %3.2f mx2: %3.2f sample_point: %2d\n", ft_est, timing_est, timing_mx1, timing_mx2, sample_point); + printf(" mx1: %3.2f mx2: %3.2f coarse_foff: %4.1f foff: %4.1f\n", timing_mx1, timing_mx2, coarse_foff_est_hz, foff_est_hz); end end @@ -504,7 +526,8 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, states.rxbuf = rxbuf; states.nin = nin; states.timing_valid = timing_valid; - states.timing_mx = timing_mx; + states.timing_mx1 = timing_mx1; + states.timing_mx2 = timing_mx2; states.timing_est = timing_est; states.sample_point = sample_point; states.delta_t = delta_t; -- 2.25.1