moving freq est routines to functions, work OK, need to rename a few vars now
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 27 Mar 2015 01:25:54 +0000 (01:25 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 27 Mar 2015 01:25:54 +0000 (01:25 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2092 01035d8c-6547-0410-b346-abe4f91aad63

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

index 182737c0b484ac8423d87089674f0247b8928781..a8e990c9908976c49dd4fd12b93072451e2a005a 100644 (file)
@@ -63,17 +63,14 @@ function sim_in = symbol_rate_init(sim_in)
 
     sim_in.Nsymb         = Nsymb            = framesize/bps;
     sim_in.Nsymbrow      = Nsymbrow         = Nsymb/Nc;
-    sim_in.Npilotsframe  = Npilotsframe     = Nsymbrow/Ns;
-    sim_in.Nsymbrowpilot = Nsymbrowpilot    = Nsymbrow + Npilotsframe + 1;
+    sim_in.Npilotsframe  = Npilotsframe     = Nsymbrow/Ns+1;
+    sim_in.Nsymbrowpilot = Nsymbrowpilot    = Nsymbrow + Npilotsframe;
 
-    printf("Each frame is %d bits or %d symbols, transmitted as %d symbols by %d carriers.",
-           framesize, Nsymb, Nsymbrow, Nc);
-    printf("  There are %d pilot symbols in each carrier, seperated by %d data/parity symbols.",
-           Npilotsframe, Ns);
-    printf("  Including pilots, the frame is %d symbols long by %d carriers.\n\n", 
-           Nsymbrowpilot, Nc);
+    printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc);
+    printf("  There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns); 
+    printf("  Including pilots, the frame is %d symbols long by %d carriers.\n\n", Nsymbrowpilot, Nc);
 
-    assert(Npilotsframe == floor(Nsymbrow/Ns), "Npilotsframe must be an integer");
+    %assert(Npilotsframe == floor(Nsymbrow/Ns), "Npilotsframe must be an integer");
 
     sim_in.prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc*Nchip);
     sim_in.prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc*Nchip);
@@ -91,6 +88,11 @@ function sim_in = symbol_rate_init(sim_in)
       write_pilot_file(pilot, Nsymbrowpilot, Ns, Np, Nsymbrow, Npilotsframe, Nc);
     end
 
+    % we use first 2 pilots of next frame to help with frame sync and fine freq
+
+    sim_in.Nct_sym_buf = 2*Nsymbrowpilot + 2;
+    sim_in.ct_symb_buf = zeros(sim_in.Nct_sym_buf, Nc);
+
     % Init LDPC --------------------------------------------------------------------
 
     if ldpc_code
@@ -157,22 +159,10 @@ function [tx_symb tx_bits prev_sym_tx] = bits_to_qpsk_symbols(sim_in, tx_bits, c
         tx_symb(r,c) = qpsk_mod(tx_bits(2*(i-1)+1:2*i));
       end
     end
-    %tx_symb = zeros(Nsymbrow,Nc);
-
-    % insert pilots, one every Ns data symbols
-
-    tx_symb_pilot = zeros(Nsymbrowpilot, Nc);
-            
-    for p=1:Npilotsframe
-      tx_symb_pilot((p-1)*(Ns+1)+1,:)          = pilot(p,:);                 % row of pilots
-      %printf("%d %d %d %d\n", (p-1)*(Ns+1)+2, p*(Ns+1), (p-1)*Ns+1, p*Ns);
-      tx_symb_pilot((p-1)*(Ns+1)+2:p*(Ns+1),:) = tx_symb((p-1)*Ns+1:p*Ns,:); % payload symbols
-    end
-    tx_symb = tx_symb_pilot;
-
-    % Append extra col of pilots at the start
-
-    tx_symb = [ pilot(1,:);  tx_symb_pilot];
+    
+    % insert pilots at start of frame
+    
+    tx_symb = [pilot(1,:); pilot(2,:); tx_symb;];
 
     % Optionally copy to other carriers (spreading)
 
@@ -490,6 +480,8 @@ endfunction
 
 % Frequency offset estimation --------------------------------------------------
 
+% This function was used in initial Nov 2014 experiments
+
 function [f_max s_max] = freq_off_est(rx_fdm, tx_pilot, offset, n)
 
   Fs = 8000;
@@ -555,6 +547,100 @@ function [f_max s_max] = freq_off_est(rx_fdm, tx_pilot, offset, n)
 endfunction
 
 
+% Set of functions to implement latest and greatest freq offset
+% estimation, March 2015 ----------------------
+
+% returns an estimate of frequency offset, advances to next sync state
+
+function [next_sync cohpsk] = coarse_freq_offset_est(cohpsk, fdmdv, ch_fdm_frame, sync, next_sync)
+  Fcentre = fdmdv.Fcentre;
+  Nc      = fdmdv.Nc;
+  Fsep    = fdmdv.Fsep;
+  M       = fdmdv.M;
+  Fs      = fdmdv.Fs;
+  Ndft    = cohpsk.Ndft;
+
+  if sync == 0
+    f_start = Fcentre - ((Nc/2)+2)*Fsep; 
+    f_stop = Fcentre + ((Nc/2)+2)*Fsep;
+    T = abs(fft(ch_fdm_frame(1:6*M).* hanning(6*M)', Ndft)).^2;
+    sc = Ndft/Fs;
+    bin_start = floor(f_start*sc+0.5)+1;
+    bin_stop = floor(f_stop*sc+0.5)+1;
+    x = bin_start-1:bin_stop-1;
+    bin_est = x*T(bin_start:bin_stop)'/sum(T(bin_start:bin_stop));
+    cohpsk.f_est = bin_est/sc;
+    printf("coarse freq est: %f\n", cohpsk.f_est);
+    next_sync = 1;
+  end
+
+endfunction
+
+
+% returns index of start of frame and fine freq offset
+
+function [next_sync cohpsk] = frame_sync_fine_timing_est(cohpsk, ch_symb, sync, next_sync)
+  ct_symb_buf   = cohpsk.ct_symb_buf;
+  Nct_sym_buf   = cohpsk.Nct_sym_buf;
+  Rs            = cohpsk.Rs;
+  Nsymbrowpilot = cohpsk.Nsymbrowpilot;
+  Nc            = cohpsk.Nc;
+
+  % update memory in symbol buffer
+
+  for r=1:Nct_sym_buf-Nsymbrowpilot
+    ct_symb_buf(r,:) = ct_symb_buf(r+Nsymbrowpilot,:);
+  end
+  i = 1;
+  for r=Nct_sym_buf-Nsymbrowpilot+1:Nct_sym_buf
+    ct_symb_buf(r,:) = ch_symb(i,:);
+    i++;
+  end
+  cohpsk.ct_symb_buf = ct_symb_buf;
+
+  % sample pilots at start of this frame and start of next frame 
+
+  sampling_points = [1 2 7 8];
+  pilot2 = [ cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);];
+
+  if sync == 2
+
+    % sample correlation over 2D grid of time and fine freq points
+
+    max_corr = 0;
+    for f_fine=-20:1:20
+      f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)';
+      for t=0:cohpsk.Nsymbrowpilot-1
+        corr = 0; mag = 0;
+        for c=1:Nc
+          f_corr_vec = f_fine_rect .* ct_symb_buf(t+sampling_points,c);
+          for p=1:length(sampling_points)
+            corr += pilot2(p,c) * f_corr_vec(p);
+            mag  += abs(f_corr_vec(p));
+          end
+        end
+        %printf("  f: %f  t: %d corr: %f %f\n", f_fine, t, real(corr), imag(corr));
+        if corr > max_corr
+          max_corr = corr;
+          max_mag = mag;
+          cohpsk.ct = t;
+          cohpsk.f_fine_est = f_fine;
+        end
+      end
+    end
+
+    printf("  fine freq f: %f max_corr: %f max_mag: %f ct: %d\n", cohpsk.f_fine_est, abs(max_corr), max_mag, cohpsk.ct);
+    if max_corr/max_mag > 0.9
+      printf("in sync!\n");
+      next_sync = 4;
+    else
+      next_sync = 0;
+      printf("  back to coarse freq offset ets...\n");
+    end
+  end
+endfunction
+
+
 % Rate Rs BER tests ------------------------------------------------------------------------------
 
 function sim_out = ber_test(sim_in)
index 263978bde5c24e68c9e563fa566d1e933621ee1c..f73639dd894e64b8df97162c4aca3e82a8536996 100644 (file)
@@ -23,18 +23,18 @@ rand('state',1);
 randn('state',1);
 
 n = 2000;
-frames = 35;
-framesize = 160
-foff = 50;
+frames = 35*4;
+framesize = 32;
+foff = -80;
 
-EsNodB = 18;
+EsNodB = 8;
 EsNo = 10^(EsNodB/10);
 variance = 1/EsNo;
 
 load ../build_linux/unittest/tcohpsk_out.txt
 
 sim_in = standard_init();
-sim_in.framesize        = 160;
+sim_in.framesize        = framesize;
 sim_in.ldpc_code        = 0;
 sim_in.ldpc_code_rate   = 1;
 sim_in.Nc = Nc          = 4;
@@ -80,7 +80,10 @@ fdmdv.phase_tx = ones(fdmdv.Nc+1,1);
 freq_hz = Fsep*( -Nc/2 - 0.5 + (1:Nc) );
 fdmdv.freq_pol = 2*pi*freq_hz/Fs;
 fdmdv.freq = exp(j*fdmdv.freq_pol);
-Fcentre = 1500;
+fdmdv.Fcentre = 1500;
+
+cohpsk.Ndft = 1024;
+cohpsk.f_est = fdmdv.Fcentre;
 
 fdmdv.fbb_rect = exp(j*2*pi*Fcentre/Fs);
 fdmdv.fbb_phase_tx = 1;
@@ -209,27 +212,15 @@ for i=1:frames
 
   % Coarse Freq offset estimation
 
-  if sync == 0
-    f_start = Fcentre - ((Nc/2)+2)*Fsep; f_stop = Fcentre + ((Nc/2)+2)*Fsep;
-    T = abs(fft(ch_fdm_frame(1:8*M).* hanning(8*M)',Fs)).^2;
-    T  = T(f_start:f_stop);
-    x = f_start:f_stop;
-    f_est = x*T'/sum(T);
-    f_off_est = f_est - Fcentre;
-    f_err = f_off_est - foff;
-    if abs(f_err) > 5
-      f_err_fail++;
-    end
-    f_err_log = [f_err_log f_err];
-    printf("coarse freq est: %f offset: %f  f_err: %f\n", f_est, f_off_est, f_err);
-    fdmdv.fbb_rect_rx = exp(j*2*pi*(f_est)/Fs);
-    sync = 1;
-  end
+  next_sync = sync;
+  
+  [next_sync cohpsk] = coarse_freq_offset_est(cohpsk, fdmdv, ch_fdm_frame, sync, next_sync);
 
   nin = M;
 
   % shift frame down to complex baseband
-
+  fdmdv.fbb_rect_rx = exp(j*2*pi*cohpsk.f_est/fdmdv.Fs);
   rx_fdm_frame_bb = zeros(1, sim_in.Nsymbrowpilot*M);
   for r=1:sim_in.Nsymbrowpilot*M
     fbb_phase_rx = fbb_phase_rx*fdmdv.fbb_rect_rx';
@@ -239,6 +230,8 @@ for i=1:frames
   fbb_phase_rx /= mag;
   rx_fdm_log = [rx_fdm_log rx_fdm_frame_bb];
   
+  % sample rate demod processing
+
   ch_symb = zeros(sim_in.Nsymbrowpilot, Nc);
   for r=1:sim_in.Nsymbrowpilot
 
@@ -255,52 +248,21 @@ for i=1:frames
   ch_symb_log = [ch_symb_log; ch_symb];
 
   % coarse timing (frame sync) and initial fine freq est ---------------------------------------------
-
-  ct_symb_buf(1:sim_in.Nsymbrowpilot,:) = ct_symb_buf(sim_in.Nsymbrowpilot+1:2*sim_in.Nsymbrowpilot,:);
-  ct_symb_buf(sim_in.Nsymbrowpilot+1:2*sim_in.Nsymbrowpilot,:) = ch_symb;
-
-  sampling_points = [1 (2:sim_in.Ns+1:1+sim_in.Npilotsframe*sim_in.Ns+1)];
-  pilot2 = [ sim_in.pilot(1,:); sim_in.pilot];
+  
+  [next_sync sim_in] = frame_sync_fine_timing_est(sim_in, ch_symb, sync, next_sync);
 
   if sync == 1
-    max_corr = 0;
-    for f_fine=-15:1:15
-      f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)';
-      for t=0:sim_in.Nsymbrowpilot-1
-        corr = 0; mag = 0;
-        for c=1:Nc
-          f_corr_vec = f_fine_rect .* ct_symb_buf(t+sampling_points,c);
-          for p=1:sim_in.Npilotsframe+1
-            corr += pilot2(p,c) * f_corr_vec(p);
-            mag  += abs(f_corr_vec(p));
-          end
-        end
-        %printf("  f: %f  t: %d corr: %f %f\n", f_fine, t, real(corr), imag(corr));
-        if corr > max_corr
-          max_corr = corr;
-          max_mag = mag;
-          ct = t;
-          f_fine_est = f_fine;
-        end
-      end
-    end
-
-    printf("  fine freq f: %f max_corr: %f max_mag: %f ct: %d\n", f_fine_est, abs(max_corr), max_mag, ct);
-    if max_corr/max_mag > 0.8
-      sync = 2;
-    end
+    next_sync = 2;
   end
-  if (i==50)
-    figure(8)
-    f_fine_rect = exp(-j*f_fine_est*2*pi*sampling_points/Rs)';
-    plot(f_fine_rect,'+');
-    hold on;
-    plot(ct_symb_buf(ct+sampling_points,1),'b+');
-    hold off; 
+
+  printf("i: %d sync: %d next_sync: %d\n", i, sync, next_sync);
+  sync = next_sync;
+
+  if (i==10)
     xx
   end
 
+  if 0
   [rx_symb rx_bits rx_symb_linear amp_linear amp_ phi_ EsNo_ prev_sym_rx sim_in] = qpsk_symbols_to_bits(sim_in, ct_symb_buf(ct+1:ct+sim_in.Nsymbrowpilot,:), []);
   rx_symb_log = [rx_symb_log; rx_symb];
   rx_amp_log = [rx_amp_log; amp_];
@@ -317,6 +279,7 @@ for i=1:frames
   end
   prev_tx_bits2 = prev_tx_bits;
   prev_tx_bits = tx_bits;
+  end
 
 end