From 8e7bf122a2173c8da67d8f28c8f5683dc69c6efb Mon Sep 17 00:00:00 2001 From: drowe67 Date: Tue, 19 May 2015 07:12:16 +0000 Subject: [PATCH] fine progress on reducing PAPR down to just 6dB in Octave sim. Es/No=12dB ading BER seems a little high at 0.03, need to look into that git-svn-id: https://svn.code.sf.net/p/freetel/code@2137 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/cohpsk.m | 22 ++++++++------- codec2-dev/octave/tcohpsk.m | 51 +++++++++++++++++++++++++---------- codec2-dev/octave/test_foff.m | 44 +++++++++++++++++++++++------- codec2-dev/octave/tpapr.m | 35 ++++++++++++++++++++++++ 4 files changed, 120 insertions(+), 32 deletions(-) create mode 100644 codec2-dev/octave/tpapr.m diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m index 31069972..01aabba8 100644 --- a/codec2-dev/octave/cohpsk.m +++ b/codec2-dev/octave/cohpsk.m @@ -217,7 +217,7 @@ function [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param) end tx_symb = tmp; - % ensures energy/symbol is normalised with diveristy + % ensures energy/symbol is normalised with diversity tx_symb = tx_symb/sqrt(Nd); end @@ -609,7 +609,7 @@ function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_proce rx_baseband = [rx_baseband arx_baseband]; rx_filt = [rx_filt arx_filt]; - rx_timng = [rx_timing arx_timing]; + rx_timing = [rx_timing arx_timing]; ch_symb(r,:) = rx_onesym; @@ -691,15 +691,17 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne max_corr = 0; for f_fine=-20:0.25:20 - f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)'; + f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)'; % note: this could be pre-computed at init or compile time for t=0:cohpsk.Nsymbrowpilot-1 corr = 0; mag = 0; for c=1:Nc*Nd - f_corr_vec = f_fine_rect .* ct_symb_buf(t+sampling_points,c); + f_corr_vec = f_fine_rect .* ct_symb_buf(t+sampling_points,c); % note: this could be pre-computed at init or compile time + acorr = 0.0; for p=1:length(sampling_points) - corr += pilot2(p,c-Nc*floor((c-1)/Nc)) * f_corr_vec(p); - mag += abs(f_corr_vec(p)); + acorr += pilot2(p,c-Nc*floor((c-1)/Nc)) * f_corr_vec(p); + mag += abs(f_corr_vec(p)); end + corr += abs(acorr); end %printf(" f: %f t: %d corr: %f %f\n", f_fine, t, real(corr), imag(corr)); if corr >= max_corr @@ -712,9 +714,9 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne end end - printf(" [%d] fine freq f: %f max_corr: %f max_mag: %f ct: %d\n", cohpsk.frame, cohpsk.f_fine_est, abs(max_corr), max_mag, cohpsk.ct); + printf(" [%d] fine freq f: %f max_ratio: %f ct: %d\n", cohpsk.frame, cohpsk.f_fine_est, abs(max_corr)/max_mag, cohpsk.ct); if abs(max_corr/max_mag) > 0.7 - printf(" [%d] encouraging sync word!\n", cohpsk.frame); + printf(" [%d] encouraging sync word! ratio: %f\n", cohpsk.frame, abs(max_corr/max_mag)); cohpsk.sync_timer = 0; next_sync = 1; else @@ -731,10 +733,12 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne f_fine_rect = exp(-j*cohpsk.f_fine_est*2*pi*sampling_points/Rs)'; for c=1:Nc*Nd f_corr_vec = f_fine_rect .* ct_symb_buf(cohpsk.ct+sampling_points,c); + acorr = 0; for p=1:length(sampling_points) - corr += pilot2(p, c-Nc*floor((c-1)/Nc)) * f_corr_vec(p); + acorr += pilot2(p, c-Nc*floor((c-1)/Nc)) * f_corr_vec(p); mag += abs(f_corr_vec(p)); end + corr += abs(acorr); end cohpsk.ratio = abs(corr)/mag; %printf("f_fine_est: %f ratio: %f\n", cohpsk.f_fine_est, cohpsk.ratio); diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m index 36a1e370..fe121053 100644 --- a/codec2-dev/octave/tcohpsk.m +++ b/codec2-dev/octave/tcohpsk.m @@ -22,6 +22,8 @@ % [X] linear interp of phase for better fading perf % [X] freq offset/drift feedback loop % [ ] PAPR measurement and reduction +% [ ] check it doesn't sync on noise +% [ ] check "unsync" graphics_toolkit ("gnuplot"); more off; @@ -35,11 +37,10 @@ randn('state',1); % test parameters ---------------------------------------------------------- -% TODO -% [ ] set up various tests we use to characterise modem for easy -% repeating when we change modem -test = 'compare to c'; +%test = 'compare to c'; +%test = 'awgn'; +test = 'fading'; if strcmp(test, 'compare to c') frames = 35; @@ -68,7 +69,7 @@ end if strcmp(test, 'fading'); frames = 100; foff = -40; - dfoff = -0.5/Fs; + dfoff = 0.5/Fs; EsNodB = 12; fading_en = 1; hf_delay_ms = 2; @@ -99,7 +100,7 @@ afdmdv.Nfilter = Nfilter; afdmdv.gt_alpha5_root = gen_rn_coeffs(0.5, 1/Fs, Rs, afdmdv.Nsym, afdmdv.M); afdmdv.Fsep = 75; afdmdv.phase_tx = ones(afdmdv.Nc+1,1); -freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd) ); +freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd) ) afdmdv.freq_pol = 2*pi*freq_hz/Fs; afdmdv.freq = exp(j*afdmdv.freq_pol); afdmdv.Fcentre = 1500; @@ -211,6 +212,11 @@ for f=1:frames [tx_fdm afdmdv] = fdm_upconvert(afdmdv, tx_baseband); tx_fdm_frame = [tx_fdm_frame tx_fdm]; end + + %clip = 50; + %ind = find(abs(tx_fdm_frame) > clip); + %tx_fdm_frame(ind) = clip*exp(j*angle(tx_fdm_frame(ind))); + tx_fdm_frame_log = [tx_fdm_frame_log tx_fdm_frame]; % @@ -323,11 +329,14 @@ for f=1:frames next_sync = 0; end + if acohpsk.ratio < 0.9 + next_sync = 0 + end if next_sync == 1 % OK we are in sync! % demodulate first frame (demod completed below) - printf(" [%d] in sync!\n", f); + printf(" [%d] in sync! f_est: %f ratio: %f \n", f, f_est, acohpsk.ratio); acohpsk.ct_symb_ff_buf(1:acohpsk.Nsymbrowpilot+2,:) = acohpsk.ct_symb_buf(acohpsk.ct+1:acohpsk.ct+acohpsk.Nsymbrowpilot+2,:); end end @@ -344,7 +353,7 @@ for f=1:frames rx_baseband_log = [rx_baseband_log rx_baseband]; rx_filt_log = [rx_filt_log rx_filt]; - ch_symb_log = [ch_symb_log; ch_symb]; + ch_symb_log = [ch_symb_log; ch_symb]; rx_timing_log = [rx_timing_log rx_timing]; f_est_log = [f_est_log acohpsk.f_est]; %printf("%f\n", acohpsk.f_est); @@ -455,30 +464,42 @@ if compare_with_c else + + papr = max(tx_fdm_frame_log.*conj(tx_fdm_frame_log)) / mean(tx_fdm_frame_log.*conj(tx_fdm_frame_log)); + papr_dB = 10*log10(papr); + printf("PAPR: %4.2f dB\n", papr_dB); % some other useful plots figure(1) + clf + subplot(211) + plot(real(tx_fdm_frame_log)) + subplot(212) + plot(imag(tx_fdm_frame_log)) + title('tx fdm'); + + figure(2) clf; plot(rx_symb_log*exp(j*pi/4),'+') title('Scatter'); - figure(2) + figure(3) clf; plot(rx_timing_log) title('rx timing'); - figure(3) + figure(4) clf; stem(nerr_log) title('Bit Errors'); - figure(4) + figure(5) clf; stem(ratio_log) title('Sync ratio'); - figure(5); + figure(6); clf; subplot(211) plot(foff_log,';freq offset;'); @@ -487,11 +508,13 @@ else hold off; title('freq offset'); legend("boxoff"); - subplot(212) - plot(foff_log(1:length(f_est_log)) - fest_log + Fcentre) + plot(foff_log(1:length(f_est_log)) - f_est_log + Fcentre) title('freq offset estimation error'); + figure(7) + clf + plot(tx_fdm_frame_log) end diff --git a/codec2-dev/octave/test_foff.m b/codec2-dev/octave/test_foff.m index 6a424ee1..25dacb9f 100644 --- a/codec2-dev/octave/test_foff.m +++ b/codec2-dev/octave/test_foff.m @@ -42,7 +42,7 @@ function sim_out = freq_off_est_test(sim_in) afdmdv.gt_alpha5_root = gen_rn_coeffs(0.5, 1/Fs, Rs, afdmdv.Nsym, afdmdv.M); afdmdv.Fsep = 75; afdmdv.phase_tx = ones(afdmdv.Nc+1,1); - freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd) ); + freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd).^1.1 ); afdmdv.freq_pol = 2*pi*freq_hz/Fs; afdmdv.freq = exp(j*afdmdv.freq_pol); afdmdv.Fcentre = 1500; @@ -50,8 +50,8 @@ function sim_out = freq_off_est_test(sim_in) afdmdv.fbb_rect = exp(j*2*pi*Fcentre/Fs); afdmdv.fbb_phase_tx = 1; afdmdv.fbb_phase_rx = 1; - afdmdv.phase_rx = ones(afdmdv.Nc+1,1); - + %afdmdv.phase_rx = ones(afdmdv.Nc+1,1); + afdmdv.phase_rx = 1 - 2*(mod(1:Nc*Nd,2)); nin = M; P = afdmdv.P = 4; @@ -103,6 +103,7 @@ function sim_out = freq_off_est_test(sim_in) ch_fdm_frame_log = []; ch_symb_log = []; tx_fdm_frame_log = []; + ratio_log = []; for f=1:frames @@ -181,10 +182,10 @@ function sim_out = freq_off_est_test(sim_in) max_ratio = 0; for acohpsk.f_est = Fcentre-40:40:Fcentre+40 + %for acohpsk.f_est = Fcentre 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 --------------------------------------------- @@ -194,7 +195,7 @@ function sim_out = freq_off_est_test(sim_in) [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); + 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; @@ -213,6 +214,7 @@ function sim_out = freq_off_est_test(sim_in) acohpsk.f_est = f_est; printf(" [%d] trying sync and f_est: %f\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; 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 @@ -228,6 +230,7 @@ function sim_out = freq_off_est_test(sim_in) printf(" [%d] in sync!\n", f); freq_offset_log = [freq_offset_log Fcentre+foff-acohpsk.f_est]; sync_time_log = [sync_time_log f-sync_start]; + ratio_log = [ratio_log max_ratio]; next_sync = 0; sync_start = f; end end @@ -244,13 +247,14 @@ function sim_out = freq_off_est_test(sim_in) sim_out.ch_fdm_frame_log = ch_fdm_frame_log; sim_out.ch_symb_log = ch_symb_log; sim_out.tx_fdm_frame_log = tx_fdm_frame_log; + sim_out.ratio_log = ratio_log; endfunction function freq_off_est_test_single - sim_in.frames = 100; + sim_in.frames = 35; sim_in.EsNodB = 12; - sim_in.foff = -59; + sim_in.foff = -55.5; sim_in.dfoff = 0; sim_in.fading_en = 1; @@ -259,14 +263,26 @@ function freq_off_est_test_single figure(1) subplot(211) plot(sim_out.freq_offset_log) + ylabel('f est error') + xlabel('time') + subplot(212) - hist(sim_out.freq_offset_log) + if length(sim_out.freq_offset_log) > 0 + hist(sim_out.freq_offset_log) + xlabel('f est error') + ylabel('histogram'); + end figure(2) subplot(211) plot(sim_out.sync_time_log) + ylabel('time to sync (frames)') subplot(212) - hist(sim_out.sync_time_log) + if length(sim_out.sync_time_log) > 0 + hist(sim_out.sync_time_log) + xlabel('time to sync (frames)') + ylabel('histogram') + end figure(3) subplot(211) @@ -275,6 +291,16 @@ function freq_off_est_test_single subplot(212) plot(real(sim_out.ch_symb_log),'+') grid; + + figure(4) + clf; + plot(sim_out.ch_symb_log,'+') + + figure(5) + clf; + plot(sim_out.ratio_log) + xlabel('time') + ylabel('ratio for sync') endfunction diff --git a/codec2-dev/octave/tpapr.m b/codec2-dev/octave/tpapr.m new file mode 100644 index 00000000..194d9a7d --- /dev/null +++ b/codec2-dev/octave/tpapr.m @@ -0,0 +1,35 @@ +% tpapr.m +% David Rowe +% 18 May 2015 + +graphics_toolkit ("gnuplot"); +rand('state',1); + +Fs = 8000; +Rs = 50; +Nc = 8; +Fc = 1500; +Fsep = ((1:Nc).^1.2)*75; +%Fsep = (1:Nc)*75 + 5 - 20*rand(1,Nc) +n = 80000; +t = 1:n; +tx = zeros(1,n); +phi = ones(1,Nc); + +figure(1) +clf + +for m=1:Nc + s = cos(phi(m)+t*2*pi*(Rs/2 + Fc + Fsep(m))/Fs); + tx += s; + subplot(Nc,1,m); + plot(s) +end + +figure(2) +plot(tx) + +tx = tx(length(tx)*0.5:length(tx)); +papr = max(tx.*conj(tx)) / mean(tx.*conj(tx)); +papr_dB = 10*log10(papr); +printf("PAPR: %4.2f dB\n", papr_dB); -- 2.25.1