changed filter order and increased dealy of new sync alg, good results on test_off...
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 6 May 2015 03:44:09 +0000 (03:44 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 6 May 2015 03:44:09 +0000 (03:44 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2125 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/tcohpsk.m
codec2-dev/octave/test_foff.m

index bc039750107136d12545509438afa26dcae7b8f2..6fdc4f246b6de1421313d2c7434c7352372a4c6d 100644 (file)
@@ -22,11 +22,11 @@ randn('state',1);
 
 % test parameters ----------------------------------------------------------
 
-frames = 10;
+frames = 100;
 foff = 0;
 dfoff = 0;
-EsNodB = 8;
-fading_en = 0;
+EsNodB = 12;
+fading_en = 1;
 hf_delay_ms = 2;
 compare_with_c = 0;
 
@@ -39,8 +39,9 @@ Nc = 4;
 Nd = 2;
 framesize = 32;
 
-afdmdv.Nsym = 2;
-afdmdv.Nt = 3;
+Nsw = 3;
+afdmdv.Nsym = 5;
+afdmdv.Nt = 5;
 
 % FDMDV init ---------------------------------------------------------------
 
@@ -94,7 +95,7 @@ acohpsk = symbol_rate_init(acohpsk);
 acohpsk.Ndft = 1024;
 acohpsk.f_est = afdmdv.Fcentre;
 
-ch_fdm_frame_buf = zeros(1, 2*acohpsk.Nsymbrowpilot*M);
+ch_fdm_frame_buf = zeros(1, Nsw*acohpsk.Nsymbrowpilot*M);
 
 % -----------------------------------------------------------
 
@@ -207,23 +208,25 @@ for f=1:frames
 
   % store two frames of received samples so we can rewind if we get a good candidate
 
-  ch_fdm_frame_buf(1:acohpsk.Nsymbrowpilot*M) = ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:2*acohpsk.Nsymbrowpilot*M);
-  ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:2*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
+  ch_fdm_frame_buf(1:(Nsw-1)*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M);
+  ch_fdm_frame_buf((Nsw-1)*acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
 
   next_sync = sync;
   nin = M;
 
   % if out of sync do Initial Freq offset estimation over last two frames
 
-  if (sync == 0) & (f>1)
+  if (sync == 0) && (f>1)
 
     % 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, 2*acohpsk.Nsymbrowpilot, nin);
-    acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb, acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot);
-    [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb, sync, next_sync);
+    [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, Nsw*acohpsk.Nsymbrowpilot, nin);
+    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
+    [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
 
@@ -235,9 +238,15 @@ 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, 2*acohpsk.Nsymbrowpilot, nin);
-      acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb, acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot);
-      [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb(acohpsk.Nsymbrowpilot+1:2*acohpsk.Nsymbrowpilot,:), sync, next_sync);
+      [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, Nsw*acohpsk.Nsymbrowpilot, nin);
+      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
+      [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:), sync, next_sync);
+      if abs(acohpsk.f_fine_est) > 2
+        printf("  [%d] Hmm %f is a bit big so back to coarse est ...\n", f, acohpsk.f_fine_est);
+        next_sync = 0;
+      end
 
       if next_sync == 1
         % OK we are in sync!
@@ -272,13 +281,14 @@ for f=1:frames
     ratio_log = [ratio_log acohpsk.ratio];
 
     % BER stats
-    %size(rx_bits)
-    %size(prev_tx_bits2)
-    %error_positions = xor(prev_tx_bits2, rx_bits);
-    error_positions = xor(prev_tx_bits, rx_bits);
-    Nerrs  += sum(error_positions);
-    nerr_log = [nerr_log sum(error_positions)];
-    Tbits += length(error_positions);
+
+    if f > 2
+      error_positions = xor(prev_tx_bits2, rx_bits);
+      %error_positions = xor(prev_tx_bits, rx_bits);
+      Nerrs  += sum(error_positions);
+      nerr_log = [nerr_log sum(error_positions)];
+      Tbits += length(error_positions);
+    end
     printf("\r  [%d]", f);
   end
 
index 242bad4e208c0a2d96d6c26f6c3631c1e1f138cf..6958090be819c1a4b4878e2557a1397cc5a0fea8 100644 (file)
-% test_foff.m
-% David Rowe April 2015
-%
-% Octave script for testing the cohpsk freq offset estimator
-
-graphics_toolkit ("gnuplot");
-more off;
-
-cohpsk;
-fdmdv;
-autotest;
-
-rand('state',1); 
-randn('state',1);
-
-% Core function for testing frequency offset estimator.  Performs one test
-
-function sim_out = freq_off_est_test(sim_in)
-  global Nfilter;
-  global M;
-
-  Rs = 50;
-  Nc = 4;
-  Nd = 2;
-  framesize = 32;
-  Fs = 8000;
-  Fcentre = 1500;
-
-  afdmdv.Nsym = 2;
-  afdmdv.Nt = 3;
-
-  afdmdv.Fs = 8000;
-  afdmdv.Nc = Nd*Nc-1;
-  afdmdv.Rs = Rs;
-  afdmdv.M  = afdmdv.Fs/afdmdv.Rs;
-  afdmdv.tx_filter_memory = zeros(afdmdv.Nc+1, Nfilter);
-  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) );
-  afdmdv.freq_pol = 2*pi*freq_hz/Fs;
-  afdmdv.freq = exp(j*afdmdv.freq_pol);
-  afdmdv.Fcentre = 1500;
-
-  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);
-
-  nin = M;
-
-  P = afdmdv.P = 4;
-  afdmdv.Nfilter = afdmdv.Nsym*afdmdv.M;
-  afdmdv.rx_filter_mem_timing = zeros(afdmdv.Nc+1, afdmdv.Nt*afdmdv.P);
-  afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M;
-  afdmdv.rx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter);
-
-  acohpsk = standard_init();
-  acohpsk.framesize        = framesize;
-  acohpsk.ldpc_code        = 0;
-  acohpsk.ldpc_code_rate   = 1;
-  acohpsk.Nc               = Nc;
-  acohpsk.Rs               = Rs;
-  acohpsk.Ns               = 4;
-  acohpsk.Nd               = Nd;
-  acohpsk.modulation       = 'qpsk';
-  acohpsk.do_write_pilot_file = 0;
-  acohpsk = symbol_rate_init(acohpsk);
-  acohpsk.Ncm  = 10*acohpsk.Nsymbrowpilot*M;
-  acohpsk.coarse_mem  = zeros(1,acohpsk.Ncm);
-  acohpsk.Ndft = 2^(ceil(log2(acohpsk.Ncm)));
-  ch_fdm_frame_buf = zeros(1, 2*acohpsk.Nsymbrowpilot*M);
-
-  frames    = sim_in.frames;
-  EsNodB    = sim_in.EsNodB;
-  foff      = sim_in.foff;
-  dfoff     = sim_in.dfoff;
-  fading_en = sim_in.fading_en;
-
-  EsNo = 10^(EsNodB/10);
-  hf_delay_ms = 2;
-  phase_ch = 1;
-
-  rand('state',1); 
-  tx_bits_coh = round(rand(1,framesize*10));
-  ptx_bits_coh = 1;
-  [spread spread_2ms hf_gain] = init_hf_model(Fs, Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M);
-
-  hf_n = 1;
-  nhfdelay = floor(hf_delay_ms*Fs/1000);
-  ch_fdm_delay = zeros(1, acohpsk.Nsymbrowpilot*M + nhfdelay);
-  
-  sync = 0; next_sync = 1;
-  sync_start = 1;
-  freq_offset_log = [];
-  sync_time_log = [];
-  ch_fdm_frame_log = [];
-  ch_symb_log = [];
-  tx_fdm_frame_log = [];
-
-  for f=1:frames
-
-    acohpsk.frame = f;
-
-    %
-    % Mod --------------------------------------------------------------------
-    %
-
-    tx_bits = tx_bits_coh(ptx_bits_coh:ptx_bits_coh+framesize-1);
-    ptx_bits_coh += framesize;
-    if ptx_bits_coh > length(tx_bits_coh)
-      ptx_bits_coh = 1;
-    end 
-
-    [tx_symb tx_bits] = bits_to_qpsk_symbols(acohpsk, tx_bits, [], []);
-
-    tx_fdm_frame = [];
-    for r=1:acohpsk.Nsymbrowpilot
-      tx_onesymb = tx_symb(r,:);
-      [tx_baseband afdmdv] = tx_filter(afdmdv, tx_onesymb);
-      [tx_fdm afdmdv] = fdm_upconvert(afdmdv, tx_baseband);
-      tx_fdm_frame = [tx_fdm_frame tx_fdm];
-    end
-    tx_fdm_frame_log = [tx_fdm_frame_log tx_fdm_frame];
-
-    %
-    % Channel --------------------------------------------------------------------
-    %
-
-    ch_fdm_frame = zeros(1,acohpsk.Nsymbrowpilot*M);
-    for i=1:acohpsk.Nsymbrowpilot*M
-      foff_rect = exp(j*2*pi*foff/Fs);
-      foff += dfoff;
-      phase_ch *= foff_rect;
-      ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch;
-    end
-    phase_ch /= abs(phase_ch);
-
-    if fading_en
-      ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M);
-      ch_fdm_delay(nhfdelay+1:nhfdelay+acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
-
-      for i=1:acohpsk.Nsymbrowpilot*M
-        ahf_model = hf_gain*(spread(hf_n)*ch_fdm_frame(i) + spread_2ms(hf_n)*ch_fdm_delay(i));
-        ch_fdm_frame(i) = ahf_model;
-        hf_n++;
-      end
-    end
-
-    % each carrier has power = 2, total power 2Nc, total symbol rate NcRs, noise BW B=Fs
-    % Es/No = (C/Rs)/(N/B), N = var = 2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No)
-
-    variance = 2*Fs/(acohpsk.Rs*EsNo);
-    uvnoise = sqrt(0.5)*(randn(1,acohpsk.Nsymbrowpilot*M) + j*randn(1,acohpsk.Nsymbrowpilot*M));
-    noise = sqrt(variance)*uvnoise;
-
-    ch_fdm_frame += noise;
-    ch_fdm_frame_log = [ch_fdm_frame_log ch_fdm_frame];
-
-    %
-    % Try to achieve sync --------------------------------------------------------------------
-    %
-
-    % store two frames so we can rewind if we get a good candidate
-
-    ch_fdm_frame_buf(1:acohpsk.Nsymbrowpilot*M) = ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:2*acohpsk.Nsymbrowpilot*M);
-    ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:2*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
-
-    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, 2*acohpsk.Nsymbrowpilot, nin);
-    ch_symb_log = [ch_symb_log; ch_symb(acohpsk.Nsymbrowpilot+1:2*acohpsk.Nsymbrowpilot,:)];
-
-    % coarse timing (frame sync) and initial fine freq est ---------------------------------------------
-  
-    acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb, acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot);
-    [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb(acohpsk.Nsymbrowpilot+1:2*acohpsk.Nsymbrowpilot,:), sync, next_sync);
-
-    % we've found a sync candidate
-
-    if (next_sync == 1)
-
-       % rewind and re-process last few frames with f_est
-
-       acohpsk.f_est -= acohpsk.f_fine_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, 2*acohpsk.Nsymbrowpilot, nin);
-       acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb, acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot);
-       [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb(acohpsk.Nsymbrowpilot+1:2*acohpsk.Nsymbrowpilot,:), sync, next_sync);
-
-       % candidate checks out so log stats
-
-       if (next_sync == 1)
-         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];
-         next_sync = 0; sync_start = f;
-       end
-    end
-
-    %printf("f: %d sync: %d next_sync: %d\n", f, sync, next_sync);
-    [sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync);
-
-  end
-
-  % ftx=fopen("coarse_tx.raw","wb"); fwrite(ftx, 1000*ch_fdm_frame_log, "short"); fclose(ftx);
-
-  sim_out.freq_offset_log = freq_offset_log;
-  sim_out.sync_time_log = sync_time_log;
-  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;
-endfunction
-
-
-function freq_off_est_test_single
-  sim_in.frames    = 10;
-  sim_in.EsNodB    = 100;
-  sim_in.foff      = 10;
-  sim_in.dfoff     = 0;
-  sim_in.fading_en = 0;
-
-  sim_out = freq_off_est_test(sim_in);
-
-  figure(1)
-  subplot(211)
-  plot(sim_out.freq_offset_log)
-  subplot(212)
-  hist(sim_out.freq_offset_log)
-
-  figure(2)
-  subplot(211)
-  plot(sim_out.sync_time_log)
-  subplot(212)
-  hist(sim_out.sync_time_log)
-
-  figure(3)
-  subplot(211)
-  plot(real(sim_out.tx_fdm_frame_log(1:2*960)))
-  subplot(212)
-  plot(real(sim_out.ch_symb_log),'+')
-endfunction
-
-
-function [freq_off_log EsNodBSet] = freq_off_est_test_curves
-  EsNodBSet = [20 12 8];
-
-  sim_in.frames    = 100;
-  sim_in.foff      = -20;
-  sim_in.dfoff     = 0;
-  sim_in.fading_en = 1;
-  freq_off_log = 1E6*ones(sim_in.frames, length(EsNodBSet) );
-  sync_time_log = 1E6*ones(sim_in.frames, length(EsNodBSet) );
-
-  for i=1:length(EsNodBSet)
-    sim_in.EsNodB = EsNodBSet(i);
-    printf("%f\n", sim_in.EsNodB);
-
-    sim_out = freq_off_est_test(sim_in);
-    freq_off_log(1:length(sim_out.freq_offset_log),i) = sim_out.freq_offset_log;
-    sync_time_log(1:length(sim_out.sync_time_log),i) = sim_out.sync_time_log;
-  end
-
-  figure(1)
-  clf
-  hold on;
-  for i=1:length(EsNodBSet)
-    data = freq_off_log(find(freq_off_log(:,i) < 1E6),i);
-    s = std(data);
-    m = mean(data);
-    stdbar = [m-s; m+s];
-    plot(EsNodBSet(i), data, '+')
-    plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-')
-  end
-  hold off
-
-  axis([6 22 -25 25])
-  if sim_in.fading_en
-    title_str = sprintf('foff = %d Hz Fading', sim_in.foff);
-  else
-    title_str = sprintf('foff = %d Hz AWGN', sim_in.foff);
-  end
-  title(title_str);
-  xlabel('Es/No (dB)')
-  ylabel('freq offset error (Hz)');
-  figure(2)
-  clf
-  hold on;
-  for i=1:length(EsNodBSet)
-    leg = sprintf("%d;%d dB;", i, EsNodBSet(i));
-    plot(freq_off_log(find(freq_off_log(:,i) < 1E6),i),leg)
-  end
-  hold off;
-  title(title_str);
-  xlabel('test');
-  ylabel('freq offset error (Hz)');
-  legend('boxoff');
-
-  figure(3)
-  clf
-  hold on;
-  for i=1:length(EsNodBSet)
-    data = sync_time_log(find(sync_time_log(:,i) < 1E6),i);
-    if length(data) 
-      s = std(data);
-      m = mean(data);
-      stdbar = [m-s; m+s];
-      plot(EsNodBSet(i), data, '+')
-      plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-')
-    end
-  end 
-  hold off;
-  axis([6 22 0 10])
-  ylabel('sync time (frames)')
-  xlabel('Es/No (dB)');
-  title(title_str);
-
-  figure(4)
-  clf
-  hold on;
-  for i=1:length(EsNodBSet)
-    leg = sprintf("%d;%d dB;", i, EsNodBSet(i));
-    plot(sync_time_log(find(sync_time_log(:,i) < 1E6),i),leg)
-  end
-  hold off;
-  title(title_str);
-  xlabel('test');
-  ylabel('sync time (frames)');
-  legend('boxoff');
-
-endfunction
-
-
-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
-%    Need to see if we get false sync
-% 3. Try shortened filter
-% 4. Extend to parallel demods at +/- 
+% test_foff.m\r
+% David Rowe April 2015\r
+%\r
+% Octave script for testing the cohpsk freq offset estimator\r
+\r
+graphics_toolkit ("gnuplot");\r
+more off;\r
+\r
+cohpsk;\r
+fdmdv;\r
+autotest;\r
+\r
+rand('state',1); \r
+randn('state',1);\r
+\r
\r
+% Core function for testing frequency offset estimator.  Performs one test\r
+\r
+function sim_out = freq_off_est_test(sim_in)\r
+  global Nfilter;\r
+  global M;\r
+\r
+  Rs = 50;\r
+  Nc = 4;\r
+  Nd = 2;\r
+  framesize = 32;\r
+  Fs = 8000;\r
+  Fcentre = 1500;\r
+  \r
+  Nsw = 3; % numbers of sync window frames we process over.  Set based\r
+           % on Nsym to flush filter memory by final frame in windw\r
+\r
+  afdmdv.Nsym = 5;\r
+  afdmdv.Nt = 5;\r
+\r
+  afdmdv.Fs = 8000;\r
+  afdmdv.Nc = Nd*Nc-1;\r
+  afdmdv.Rs = Rs;\r
+  afdmdv.M  = afdmdv.Fs/afdmdv.Rs;\r
+  afdmdv.tx_filter_memory = zeros(afdmdv.Nc+1, Nfilter);\r
+  afdmdv.Nfilter =  Nfilter;\r
+  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
+  afdmdv.freq_pol = 2*pi*freq_hz/Fs;\r
+  afdmdv.freq = exp(j*afdmdv.freq_pol);\r
+  afdmdv.Fcentre = 1500;\r
+\r
+  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
+  nin = M;\r
+\r
+  P = afdmdv.P = 4;\r
+  afdmdv.Nfilter = afdmdv.Nsym*afdmdv.M;\r
+  afdmdv.rx_filter_mem_timing = zeros(afdmdv.Nc+1, afdmdv.Nt*afdmdv.P);\r
+  afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M;\r
+  afdmdv.rx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter);\r
+\r
+  acohpsk = standard_init();\r
+  acohpsk.framesize        = framesize;\r
+  acohpsk.ldpc_code        = 0;\r
+  acohpsk.ldpc_code_rate   = 1;\r
+  acohpsk.Nc               = Nc;\r
+  acohpsk.Rs               = Rs;\r
+  acohpsk.Ns               = 4;\r
+  acohpsk.Nd               = Nd;\r
+  acohpsk.modulation       = 'qpsk';\r
+  acohpsk.do_write_pilot_file = 0;\r
+  acohpsk = symbol_rate_init(acohpsk);\r
+  acohpsk.Ncm  = 10*acohpsk.Nsymbrowpilot*M;\r
+  acohpsk.coarse_mem  = zeros(1,acohpsk.Ncm);\r
+  acohpsk.Ndft = 2^(ceil(log2(acohpsk.Ncm)));\r
\r
+  ch_fdm_frame_buf = zeros(1, Nsw*acohpsk.Nsymbrowpilot*M);\r
+\r
+  frames    = sim_in.frames;\r
+  EsNodB    = sim_in.EsNodB;\r
+  foff      = sim_in.foff;\r
+  dfoff     = sim_in.dfoff;\r
+  fading_en = sim_in.fading_en;\r
+\r
+  EsNo = 10^(EsNodB/10);\r
+  hf_delay_ms = 2;\r
+  phase_ch = 1;\r
+\r
+  rand('state',1); \r
+  tx_bits_coh = round(rand(1,framesize*10));\r
+  ptx_bits_coh = 1;\r
+  [spread spread_2ms hf_gain] = init_hf_model(Fs, Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M);\r
+\r
+  hf_n = 1;\r
+  nhfdelay = floor(hf_delay_ms*Fs/1000);\r
+  ch_fdm_delay = zeros(1, acohpsk.Nsymbrowpilot*M + nhfdelay);\r
+  \r
+  sync = 0; next_sync = 1;\r
+  sync_start = 1;\r
+  freq_offset_log = [];\r
+  sync_time_log = [];\r
+  ch_fdm_frame_log = [];\r
+  ch_symb_log = [];\r
+  tx_fdm_frame_log = [];\r
+\r
+  for f=1:frames\r
+\r
+    acohpsk.frame = f;\r
+\r
+    %\r
+    % Mod --------------------------------------------------------------------\r
+    %\r
+\r
+    tx_bits = tx_bits_coh(ptx_bits_coh:ptx_bits_coh+framesize-1);\r
+    ptx_bits_coh += framesize;\r
+    if ptx_bits_coh > length(tx_bits_coh)\r
+      ptx_bits_coh = 1;\r
+    end \r
+\r
+    [tx_symb tx_bits] = bits_to_qpsk_symbols(acohpsk, tx_bits, [], []);\r
+\r
+    tx_fdm_frame = [];\r
+    for r=1:acohpsk.Nsymbrowpilot\r
+      tx_onesymb = tx_symb(r,:);\r
+      [tx_baseband afdmdv] = tx_filter(afdmdv, tx_onesymb);\r
+      [tx_fdm afdmdv] = fdm_upconvert(afdmdv, tx_baseband);\r
+      tx_fdm_frame = [tx_fdm_frame tx_fdm];\r
+    end\r
+    tx_fdm_frame_log = [tx_fdm_frame_log tx_fdm_frame];\r
+\r
+    %\r
+    % Channel --------------------------------------------------------------------\r
+    %\r
+\r
+    ch_fdm_frame = zeros(1,acohpsk.Nsymbrowpilot*M);\r
+    for i=1:acohpsk.Nsymbrowpilot*M\r
+      foff_rect = exp(j*2*pi*foff/Fs);\r
+      foff += dfoff;\r
+      phase_ch *= foff_rect;\r
+      ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch;\r
+    end\r
+    phase_ch /= abs(phase_ch);\r
+\r
+    if fading_en\r
+      ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M);\r
+      ch_fdm_delay(nhfdelay+1:nhfdelay+acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;\r
+\r
+      for i=1:acohpsk.Nsymbrowpilot*M\r
+        ahf_model = hf_gain*(spread(hf_n)*ch_fdm_frame(i) + spread_2ms(hf_n)*ch_fdm_delay(i));\r
+        ch_fdm_frame(i) = ahf_model;\r
+        hf_n++;\r
+      end\r
+    end\r
+\r
+    % each carrier has power = 2, total power 2Nc, total symbol rate NcRs, noise BW B=Fs\r
+    % Es/No = (C/Rs)/(N/B), N = var = 2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No)\r
+\r
+    variance = 2*Fs/(acohpsk.Rs*EsNo);\r
+    uvnoise = sqrt(0.5)*(randn(1,acohpsk.Nsymbrowpilot*M) + j*randn(1,acohpsk.Nsymbrowpilot*M));\r
+    noise = sqrt(variance)*uvnoise;\r
+\r
+    ch_fdm_frame += noise;\r
+    ch_fdm_frame_log = [ch_fdm_frame_log ch_fdm_frame];\r
+\r
+    %\r
+    % Try to achieve sync --------------------------------------------------------------------\r
+    %\r
+\r
+    % store Nsw frames so we can rewind if we get a good candidate\r
+\r
+    ch_fdm_frame_buf(1:(Nsw-1)*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M);\r
+    ch_fdm_frame_buf((Nsw-1)*acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;\r
+\r
+    next_sync = sync;\r
+    sync = 0;\r
+    if sync == 0\r
+      acohpsk.f_est = Fcentre;\r
+    end\r
+\r
+    [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame_buf, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx);\r
+    [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, Nsw*acohpsk.Nsymbrowpilot, nin);\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
+    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
+    [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:), sync, next_sync);\r
+\r
+    % we've found a sync candidate\r
+\r
+    if (next_sync == 1)\r
+\r
+       % rewind and re-process last few frames with f_est\r
+\r
+       acohpsk.f_est -= acohpsk.f_fine_est;\r
+       printf("  [%d] trying sync and f_est: %f\n", f, acohpsk.f_est);\r
+       [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame_buf, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx);\r
+       [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, Nsw*acohpsk.Nsymbrowpilot, nin);\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
+       [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb((Nsw-1)*acohpsk.Nsymbrowpilot+1:Nsw*acohpsk.Nsymbrowpilot,:), sync, next_sync);\r
+       if abs(acohpsk.f_fine_est) > 2\r
+         printf("  [%d] Hmm %f is a bit big so back to coarse est ...\n", f, acohpsk.f_fine_est);\r
+         next_sync = 0;\r
+       end\r
+\r
+       % candidate checks out so log stats\r
+\r
+       if (next_sync == 1)\r
+         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
+         next_sync = 0; sync_start = f;\r
+       end\r
+    end\r
+\r
+    %printf("f: %d sync: %d next_sync: %d\n", f, sync, next_sync);\r
+    [sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync);\r
+\r
+  end\r
+\r
+  % ftx=fopen("coarse_tx.raw","wb"); fwrite(ftx, 1000*ch_fdm_frame_log, "short"); fclose(ftx);\r
+\r
+  sim_out.freq_offset_log = freq_offset_log;\r
+  sim_out.sync_time_log = sync_time_log;\r
+  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
+endfunction\r
+\r
+\r
+function freq_off_est_test_single\r
+  sim_in.frames    = 100;\r
+  sim_in.EsNodB    = 12;\r
+  sim_in.foff      = -20;\r
+  sim_in.dfoff     = 0;\r
+  sim_in.fading_en = 1;\r
+\r
+  sim_out = freq_off_est_test(sim_in);\r
+\r
+  figure(1)\r
+  subplot(211)\r
+  plot(sim_out.freq_offset_log)\r
+  subplot(212)\r
+  hist(sim_out.freq_offset_log)\r
+\r
+  figure(2)\r
+  subplot(211)\r
+  plot(sim_out.sync_time_log)\r
+  subplot(212)\r
+  hist(sim_out.sync_time_log)\r
+\r
+  figure(3)\r
+  subplot(211)\r
+  plot(real(sim_out.tx_fdm_frame_log(1:2*960)))\r
+  grid;\r
+  subplot(212)\r
+  plot(real(sim_out.ch_symb_log),'+')\r
+  grid;\r
+endfunction\r
+\r
+\r
+function [freq_off_log EsNodBSet] = freq_off_est_test_curves\r
+  EsNodBSet = [20 12 8];\r
+\r
+  sim_in.frames    = 100;\r
+  sim_in.foff      = -20;\r
+  sim_in.dfoff     = 0;\r
+  sim_in.fading_en = 1;\r
+  freq_off_log = 1E6*ones(sim_in.frames, length(EsNodBSet) );\r
+  sync_time_log = 1E6*ones(sim_in.frames, length(EsNodBSet) );\r
+\r
+  for i=1:length(EsNodBSet)\r
+    sim_in.EsNodB = EsNodBSet(i);\r
+    printf("%f\n", sim_in.EsNodB);\r
+\r
+    sim_out = freq_off_est_test(sim_in);\r
+    freq_off_log(1:length(sim_out.freq_offset_log),i) = sim_out.freq_offset_log;\r
+    sync_time_log(1:length(sim_out.sync_time_log),i) = sim_out.sync_time_log;\r
+  end\r
+\r
+  figure(1)\r
+  clf\r
+  hold on;\r
+  for i=1:length(EsNodBSet)\r
+    data = freq_off_log(find(freq_off_log(:,i) < 1E6),i);\r
+    s = std(data);\r
+    m = mean(data);\r
+    stdbar = [m-s; m+s];\r
+    plot(EsNodBSet(i), data, '+')\r
+    plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-')\r
+  end\r
+  hold off\r
+\r
+  axis([6 22 -25 25])\r
+  if sim_in.fading_en\r
+    title_str = sprintf('foff = %d Hz Fading', sim_in.foff);\r
+  else\r
+    title_str = sprintf('foff = %d Hz AWGN', sim_in.foff);\r
+  end\r
+  title(title_str);\r
+  xlabel('Es/No (dB)')\r
+  ylabel('freq offset error (Hz)');\r
\r
+  figure(2)\r
+  clf\r
+  hold on;\r
+  for i=1:length(EsNodBSet)\r
+    leg = sprintf("%d;%d dB;", i, EsNodBSet(i));\r
+    plot(freq_off_log(find(freq_off_log(:,i) < 1E6),i),leg)\r
+  end\r
+  hold off;\r
+  title(title_str);\r
+  xlabel('test');\r
+  ylabel('freq offset error (Hz)');\r
+  legend('boxoff');\r
+\r
+  figure(3)\r
+  clf\r
+  hold on;\r
+  for i=1:length(EsNodBSet)\r
+    data = sync_time_log(find(sync_time_log(:,i) < 1E6),i);\r
+    if length(data) \r
+      s = std(data);\r
+      m = mean(data);\r
+      stdbar = [m-s; m+s];\r
+      plot(EsNodBSet(i), data, '+')\r
+      plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-')\r
+    end\r
+  end \r
+  hold off;\r
+  axis([6 22 0 10])\r
+  ylabel('sync time (frames)')\r
+  xlabel('Es/No (dB)');\r
+  title(title_str);\r
+\r
+  figure(4)\r
+  clf\r
+  hold on;\r
+  for i=1:length(EsNodBSet)\r
+    leg = sprintf("%d;%d dB;", i, EsNodBSet(i));\r
+    plot(sync_time_log(find(sync_time_log(:,i) < 1E6),i),leg)\r
+  end\r
+  hold off;\r
+  title(title_str);\r
+  xlabel('test');\r
+  ylabel('sync time (frames)');\r
+  legend('boxoff');\r
+\r
+endfunction\r
+\r
+\r
+%freq_off_est_test_single;\r
+freq_off_est_test_curves;\r
+\r
+% 1. start with +/- 20Hz offset\r
+% 2. Measure frames to sync.  How to define sync?  Foff to withn 1 Hz. Sync state\r
+%    Need to see if we get false sync\r
+% 3. Try shortened filter\r
+% 4. Extend to parallel demods at +/- \r