fine progress on reducing PAPR down to just 6dB in Octave sim. Es/No=12dB ading...

git-svn-id: https://svn.code.sf.net/p/freetel/code@2137 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m
index 3106997208bd3fdbefc02351cd7ea6781d3feb48..01aabba877e0d1b5492fabb7fa8c7b07a64308a5 100644 (file)
--- 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 36a1e37095969766d42c4bd338342e3ae4d1bfc6..fe121053fb1a195108f6ee86a3e7ba983c08d656 100644 (file)
--- 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 6a424ee1ad0bdbc15bbd6f4924585a054288af24..25dacb9ff5d2b1a3657a9ae0ddb6e0fb555e1355 100644 (file)
--- 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);\r
   afdmdv.Fsep = 75;\r
   afdmdv.phase_tx = ones(afdmdv.Nc+1,1);\r
-  freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd) );\r
+  freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd).^1.1 );\r
   afdmdv.freq_pol = 2*pi*freq_hz/Fs;\r
   afdmdv.freq = exp(j*afdmdv.freq_pol);\r
   afdmdv.Fcentre = 1500;\r
@@ -50,8 +50,8 @@ function sim_out = freq_off_est_test(sim_in)
   afdmdv.fbb_rect = exp(j*2*pi*Fcentre/Fs);\r
   afdmdv.fbb_phase_tx = 1;\r
   afdmdv.fbb_phase_rx = 1;\r
-  afdmdv.phase_rx = ones(afdmdv.Nc+1,1);\r
-\r
+  %afdmdv.phase_rx = ones(afdmdv.Nc+1,1);\r
+  afdmdv.phase_rx = 1 - 2*(mod(1:Nc*Nd,2));\r
   nin = M;\r
 \r
   P = afdmdv.P = 4;\r
@@ -103,6 +103,7 @@ function sim_out = freq_off_est_test(sim_in)
   ch_fdm_frame_log = [];\r
   ch_symb_log = [];\r
   tx_fdm_frame_log = [];\r
+  ratio_log = [];\r
 \r
   for f=1:frames\r
 \r
@@ -181,10 +182,10 @@ function sim_out = freq_off_est_test(sim_in)
 \r
       max_ratio = 0;\r
       for acohpsk.f_est = Fcentre-40:40:Fcentre+40\r
+      %for acohpsk.f_est = Fcentre\r
         \r
         printf("  [%d] acohpsk.f_est: %f +/- 20\n", f, acohpsk.f_est);\r
         [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);\r
-        ch_symb_log = [ch_symb_log; ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:)];\r
 \r
         % coarse timing (frame sync) and initial fine freq est ---------------------------------------------\r
   \r
@@ -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);\r
 \r
         if anext_sync == 1\r
-          %printf("  [%d] acohpsk.ratio: %f\n", f, acohpsk.ratio);\r
+          printf("  [%d] acohpsk.ratio: %f\n", f, acohpsk.ratio);\r
           if acohpsk.ratio > max_ratio\r
             max_ratio   = acohpsk.ratio;\r
             f_est       = acohpsk.f_est - acohpsk.f_fine_est;\r
@@ -213,6 +214,7 @@ function sim_out = freq_off_est_test(sim_in)
        acohpsk.f_est = f_est;\r
        printf("  [%d] trying sync and f_est: %f\n", f, acohpsk.f_est);\r
        [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);\r
+       ch_symb_log = ch_symb;\r
        for i=1:Nsw-1\r
          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);\r
        end\r
@@ -228,6 +230,7 @@ function sim_out = freq_off_est_test(sim_in)
          printf("  [%d] in sync!\n", f);\r
          freq_offset_log = [freq_offset_log Fcentre+foff-acohpsk.f_est];\r
          sync_time_log = [sync_time_log f-sync_start];\r
+         ratio_log = [ratio_log max_ratio];\r
          next_sync = 0; sync_start = f;\r
        end\r
     end\r
@@ -244,13 +247,14 @@ function sim_out = freq_off_est_test(sim_in)
   sim_out.ch_fdm_frame_log = ch_fdm_frame_log;\r
   sim_out.ch_symb_log = ch_symb_log;\r
   sim_out.tx_fdm_frame_log = tx_fdm_frame_log;\r
+  sim_out.ratio_log = ratio_log;\r
 endfunction\r
 \r
 \r
 function freq_off_est_test_single\r
-  sim_in.frames    = 100;\r
+  sim_in.frames    = 35;\r
   sim_in.EsNodB    = 12;\r
-  sim_in.foff      = -59;\r
+  sim_in.foff      = -55.5;\r
   sim_in.dfoff     = 0;\r
   sim_in.fading_en = 1;\r
 \r
@@ -259,14 +263,26 @@ function freq_off_est_test_single
   figure(1)\r
   subplot(211)\r
   plot(sim_out.freq_offset_log)\r
+  ylabel('f est error')\r
+  xlabel('time')\r
+\r
   subplot(212)\r
-  hist(sim_out.freq_offset_log)\r
+  if length(sim_out.freq_offset_log) > 0\r
+    hist(sim_out.freq_offset_log)\r
+    xlabel('f est error')\r
+    ylabel('histogram');\r
+  end\r
 \r
   figure(2)\r
   subplot(211)\r
   plot(sim_out.sync_time_log)\r
+  ylabel('time to sync (frames)')\r
   subplot(212)\r
-  hist(sim_out.sync_time_log)\r
+  if length(sim_out.sync_time_log) > 0\r
+    hist(sim_out.sync_time_log)\r
+    xlabel('time to sync (frames)')\r
+    ylabel('histogram')\r
+  end\r
 \r
   figure(3)\r
   subplot(211)\r
@@ -275,6 +291,16 @@ function freq_off_est_test_single
   subplot(212)\r
   plot(real(sim_out.ch_symb_log),'+')\r
   grid;\r
+\r
+  figure(4)\r
+  clf;\r
+  plot(sim_out.ch_symb_log,'+')\r
+\r
+  figure(5)\r
+  clf;\r
+  plot(sim_out.ratio_log)\r
+  xlabel('time')\r
+  ylabel('ratio for sync')\r
 endfunction\r
 \r
 \r

diff --git a/codec2-dev/octave/tpapr.m b/codec2-dev/octave/tpapr.m
new file mode 100644 (file)
index 0000000..194d9a7
--- /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);