characterising freq offset estimation metrics
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sat, 7 Apr 2018 20:38:06 +0000 (20:38 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sat, 7 Apr 2018 20:38:06 +0000 (20:38 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@3447 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/ofdm_dev.m
codec2-dev/octave/ofdm_lib.m

index d39bbef81eba5e21c5c79fa9def6f4c5517d8eb0..40f17ac68d14fa2acdde4a8081d867ab81e6dfc9 100644 (file)
@@ -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')
index 5495a465831da8e51000e9a5c67cf48d1fc6b2f7..b41b93c35ad565874bb4dbd4d5f97a828a2f1c74 100644 (file)
@@ -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;