cohpsk freq tracking unit test, integrated into tcohpsk and doing sensible things...

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

diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m
index f5ccd4332dc5c11caf9356fb2185f9b95edc40af..4d06ba80884b2a54ae1046eccd5996c9bcf5a8e3 100644 (file)
--- a/codec2-dev/octave/cohpsk.m
+++ b/codec2-dev/octave/cohpsk.m
@@ -589,17 +589,56 @@ function [next_sync cohpsk] = coarse_freq_offset_est(cohpsk, fdmdv, ch_fdm_frame
 endfunction
 
 
-function [ch_symb rx_timing rx_filt rx_baseband afdmdv] = rate_Fs_rx_processing(afdmdv, rx_fdm_frame_bb, nsymb, nin)
+function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, f_est, nsymb, nin, freq_track)
     M = afdmdv.M;
 
+    % figure(10); clf; hold on;
     for r=1:nsymb
+      % downconvert entire signal to nominal baseband 
+
+      [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame(1+(r-1)*M:r*M), -f_est, afdmdv.Fs, afdmdv.fbb_phase_rx);
+
       % downconvert each FDM carrier to Nc separate baseband signals
 
-      [rx_baseband afdmdv] = fdm_downconvert(afdmdv, rx_fdm_frame_bb(1+(r-1)*M:r*M), nin);
+      [rx_baseband afdmdv] = fdm_downconvert(afdmdv, rx_fdm_frame_bb, nin);
       [rx_filt afdmdv] = rx_filter(afdmdv, rx_baseband, nin);
       [rx_onesym rx_timing env afdmdv] = rx_est_timing(afdmdv, rx_filt, nin);     
 
       ch_symb(r,:) = rx_onesym;
+
+      % freq tracking, see test_ftrack.m for unit test.  Placed in
+      % this function as it needs to work on a symbol by symbol basis
+      % rather than frame by frame.  This means the control loop
+      % operates at a sample rate of Rs = 50Hz for say 1 Hz/s drift.
+
+      if freq_track
+        beta = 0.005;
+        g = 0.2;
+
+        % combine difference on phase from last symbol over Nc carriers
+
+        mod_strip = 0;
+        for c=1:afdmdv.Nc+1
+          adiff = rx_onesym(c) .* conj(afdmdv.prev_rx_symb(c));
+          afdmdv.prev_rx_symb(c) = rx_onesym(c);
+          amod_strip = adiff.^4;
+          amod_strip = abs(real(amod_strip)) + j*imag(amod_strip);
+          mod_strip += amod_strip;
+        end
+        %plot(mod_strip)
+
+        % 4th power strips QPSK modulation, by multiplying phase by 4
+        % Using the abs value of the real coord was found to help 
+        % non-linear issues when noise power was large
+
+        
+        % loop filter made up of 1st order IIR plus integrator.  Integerator
+        % was found to be reqd 
+        
+        afdmdv.filt = ((1-beta)*afdmdv.filt + beta*angle(mod_strip));
+        f_est += g*afdmdv.filt;
+
+      end
     end
 endfunction
 
@@ -683,7 +722,7 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne
     max_corr = 0;
     st = cohpsk.f_fine_est - 1;
     en = cohpsk.f_fine_est + 1;
-    for f_fine = st:0.25*en
+    for f_fine = st:0.25:en
         f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)';
         corr = 0; mag = 0;
         for c=1:Nc*Nd
@@ -764,10 +803,10 @@ function [sync cohpsk] = sync_state_machine(cohpsk, sync, next_sync)
     else
       cohpsk.sync_timer = 0;            
     end
-    %printf("  ratio: %f  sync timer: %d\n", cohpsk.ratio, cohpsk.sync_timer);
+    % printf("  ratio: %f  sync timer: %d\n", cohpsk.ratio, cohpsk.sync_timer);
 
-    if cohpsk.sync_timer == 5
-      printf("  lost sync ....\n");
+    if cohpsk.sync_timer == 10
+      printf("  [%d] lost sync ....\n", cohpsk.frame);
       next_sync = 0;
     end
   end

diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m
index 6fdc4f246b6de1421313d2c7434c7352372a4c6d..2a9cb7cd9e691a027a424090c81370a5a17a65de 100644 (file)
--- a/codec2-dev/octave/tcohpsk.m
+++ b/codec2-dev/octave/tcohpsk.m
@@ -22,9 +22,9 @@ randn('state',1);
 
 % test parameters ----------------------------------------------------------
 
-frames = 100;
-foff = 0;
-dfoff = 0;
+frames = 200;
+foff = -20;
+dfoff = 0/Fs;
 EsNodB = 12;
 fading_en = 1;
 hf_delay_ms = 2;
@@ -78,6 +78,9 @@ afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M;
 
 afdmdv.rx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter);
 
+afdmdv.filt = 0;
+afdmdv.prev_rx_symb = zeros(1,afdmdv.Nc+1);
+
 % COHPSK Init --------------------------------------------------------
 
 acohpsk = standard_init();
@@ -120,6 +123,8 @@ rx_fdm_frame_log = [];
 ct_symb_ff_log = [];
 rx_timing_log = [];
 ratio_log = [];
+foff_log = [];
+fest_log = [];
 
 % Channel modeling and BER measurement ----------------------------------------
 
@@ -177,7 +182,9 @@ for f=1:frames
     phase_ch *= foff_rect;
     ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch;
   end
+  foff_log = [foff_log foff];
   phase_ch /= abs(phase_ch);
+  % printf("foff: %f  ", foff);
 
   if fading_en
     ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M);
@@ -221,8 +228,8 @@ for f=1:frames
     % 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, Nsw*acohpsk.Nsymbrowpilot, nin);
+    %[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 f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame_buf, acohpsk.f_est, Nsw*acohpsk.Nsymbrowpilot, nin, 0);
     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
@@ -237,8 +244,8 @@ 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, Nsw*acohpsk.Nsymbrowpilot, nin);
+      %[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 f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame_buf, acohpsk.f_est, Nsw*acohpsk.Nsymbrowpilot, nin, 0);
       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
@@ -261,12 +268,17 @@ for f=1:frames
   % If in sync just do sample rate processing on latest frame
 
   if sync == 1
-    [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame, -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, acohpsk.Nsymbrowpilot, nin);
-    acohpsk.ct_symb_buf = update_ct_symb_buf(acohpsk.ct_symb_buf, ch_symb, acohpsk.Nct_sym_buf, acohpsk.Nsymbrowpilot);
+    %[rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame, -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, acohpsk.Nsymbrowpilot, nin);
+    [ch_symb rx_timing rx_filt rx_baseband afdmdv acohpsk.f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, acohpsk.f_est, acohpsk.Nsymbrowpilot, nin, 1);
+    [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb, sync, next_sync);
 
     acohpsk.ct_symb_ff_buf(1:2,:) = acohpsk.ct_symb_ff_buf(acohpsk.Nsymbrowpilot+1:acohpsk.Nsymbrowpilot+2,:);
     acohpsk.ct_symb_ff_buf(3:acohpsk.Nsymbrowpilot+2,:) = acohpsk.ct_symb_buf(acohpsk.ct+3:acohpsk.ct+acohpsk.Nsymbrowpilot+2,:);
+
+    ch_symb_log = [ch_symb_log; ch_symb];
+    rx_timing_log = [rx_timing_log rx_timing];
+    fest_log = [fest_log acohpsk.f_est];
   end
 
   % if we are in sync complete demodulation with symbol rate processing
@@ -373,12 +385,32 @@ else
 
   figure(1)
   plot(rx_symb_log*exp(j*pi/4),'+')
+  title('Scatter');
+
   figure(2)
   plot(rx_timing_log)
+  title('rx timing');
+
   figure(3)
   stem(nerr_log)
+  title('Bit Errors');
+
   figure(4)
   stem(ratio_log)
+  title('Sync ratio');
+
+  figure(5);
+  subplot(211)
+  plot(foff_log);
+  hold on;
+  plot(fest_log - Fcentre,'g');
+  hold off;
+  
+  title('freq offset error');
+
+  subplot(212)
+  plot(foff_log(1:length(fest_log)) - fest_log + Fcentre)
+  title('freq offset estimation error');
 
 end
 

diff --git a/codec2-dev/octave/test_ftrack.m b/codec2-dev/octave/test_ftrack.m
new file mode 100644 (file)
index 0000000..5b02e52
--- /dev/null
+++ b/codec2-dev/octave/test_ftrack.m
@@ -0,0 +1,168 @@
+% test_ftrack.m
+%
+% David Rowe May 2015
+%
+% Octave script that to test and develop symbol rate frequency offset
+% tracking
+%
+% [ ] add Es/No noise and measure BER
+% [ ] QPSK mod of random symbols
+% [ ] tracking 1 Hz/s freq drift with negl impl loss
+% [ ] simulate delay through filter, or implement R/N filter
+
+rand('state',1); 
+randn('state',1);
+graphics_toolkit ("gnuplot");
+
+Fs = 8000;
+Rs = 50;
+Nbits = 5000;
+Nsymb = Nbits/2;
+foff = 1;
+dfoff = -1/Fs;
+M = Fs/Rs;
+
+phase_ch = 1;
+prev_rx_symb = 1;
+
+EsNodB = 8;
+EsNo = 10^(EsNodB/10);    
+variance = 1/EsNo;
+
+Nsym = 5;
+
+beta = 0.005;
+g = 0.2;
+ftrack = 0;
+filt = 0;
+
+% Second order IIR filter coeffs, derived with pencil and paper
+
+a1 = (2-beta)/(1+g*beta);
+a2 = (1-beta)/(1+g*beta);
+
+% Which can be used to work out the polar coordinates of the pole
+
+gamma = sqrt(a2);           % radius (distance from origin)
+w = acos(a1/(2*gamma));     % angular frequency
+
+printf("2nd order loop w: %f   gamma: %f\n", w, gamma);
+
+function symbol = qpsk_mod(two_bits)
+    two_bits_decimal = sum(two_bits .* [2 1]); 
+    switch(two_bits_decimal)
+        case (0) symbol =  1;
+        case (1) symbol =  j;
+        case (2) symbol = -j;
+        case (3) symbol = -1;
+    endswitch
+endfunction
+
+% Gray coded QPSK demodulation function
+
+function two_bits = qpsk_demod(symbol)
+    if isscalar(symbol) == 0
+        printf("only works with scalars\n");
+        return;
+    end
+    bit0 = real(symbol*exp(j*pi/4)) < 0;
+    bit1 = imag(symbol*exp(j*pi/4)) < 0;
+    two_bits = [bit1 bit0];
+endfunction
+
+% modulator
+
+tx_bits = rand(1, Nbits) > 0.5;
+tx_symb = zeros(1, Nsymb);
+for i=1:Nsymb
+  tx_symb(i) = qpsk_mod(tx_bits(2*(i-1)+1:2*i));
+end
+
+% run symbol by symbol 
+
+rx_symb = zeros(1, Nsymb);
+rx_filt = zeros(1,Nsym);
+
+diff_log = [];
+foff_log = [];
+ftrack_log = [];
+mod_strip_log = [];
+
+for i=1:Nsymb
+
+  % channel
+
+  foff_rect = exp(j*2*pi*(foff-ftrack)*M/Fs);
+  foff += M*dfoff;
+  phase_ch *= foff_rect;
+  rx_symb(i) = tx_symb(i) * phase_ch;
+  noise = sqrt(variance*0.5)*(randn(1,1) + j*randn(1, 1));
+  rx_symb(i) += noise;
+
+  % simulate delay of root nyquist filter
+
+  rx_filt(1:Nsym-1) = rx_filt(2:Nsym);
+  rx_filt(Nsym) = rx_symb(i);
+
+  % freq tracking loop
+
+  diff = rx_filt(1) .* conj(prev_rx_symb);
+  prev_rx_symb = rx_filt(1);
+
+  % 4th power strips QPSK modulation, by multiplying phase by 4
+  % Using the abs value of the real coord was found to help 
+  % non-linear issues when noise power was large
+
+  mod_strip = diff.^4;
+  mod_strip = abs(real(mod_strip)) + j*imag(mod_strip);
+
+  % loop filter made up of 1st order IIR plus integrator.  Integerator
+  % was found to be reqd 
+  filt = ((1-beta)*filt + beta*angle(mod_strip));
+  ftrack += g*filt;
+
+  diff_log = [diff_log diff];
+  foff_log = [foff_log foff];
+  ftrack_log = [ftrack_log ftrack];
+  mod_strip_log = [mod_strip_log mod_strip];
+end
+
+% plot results
+
+figure(1)
+clf
+plot(angle(mod_strip_log))
+title('mod stripping phase');
+
+figure(2);
+clf
+plot(mod_strip_log,'+')
+axis([-2 2 -2 2])
+title('mod stripping scatter');
+
+figure(3)
+clf
+subplot(211)
+plot(real(mod_strip_log))
+title('mod stripping real and imag');
+subplot(212)
+plot(imag(mod_strip_log))
+
+figure(4)
+clf
+subplot(211)
+plot(foff_log,';foff;')
+hold on
+plot(ftrack_log,'g+;ftrack;')
+hold off;
+legend("boxoff");  
+ylabel('Freq Offset Hz');
+
+subplot(212)
+plot(foff_log-ftrack_log,';foff-ftrack;' )
+ylabel('Freq Offset Tracking Error Hz');
+xlabel('symbols')
+legend("boxoff");  
+
+figure(5)
+freqz(g*beta/(1+g*beta), [1 -2*gamma*cos(w) gamma*gamma])