From ce90e01715e70c8392ffbb9ea4adc3103acd6cf9 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Sat, 19 Jan 2013 06:17:20 +0000
Subject: [PATCH] coherent demod simulations

git-svn-id: https://svn.code.sf.net/p/freetel/code@1156 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/octave/fdmdv.m           |  86 ++++++++--
 codec2-dev/octave/fdmdv_demod.m     |   2 +-
 codec2-dev/octave/fdmdv_demod_coh.m | 253 ++++++++++++++++++++++++++++
 codec2-dev/octave/fdmdv_ut.m        |   7 +-
 codec2-dev/octave/fdmdv_ut_coh.m    | 139 ++++-----------
 5 files changed, 366 insertions(+), 121 deletions(-)
 create mode 100644 codec2-dev/octave/fdmdv_demod_coh.m

diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m
index 432be81a..150a3cd5 100644
--- a/codec2-dev/octave/fdmdv.m
+++ b/codec2-dev/octave/fdmdv.m
@@ -30,7 +30,11 @@ global P = 4;          % oversample factor used for rx symbol filtering
 global Nfilter = Nsym*M;
 global Nfiltertiming = M+Nfilter+M;
 alpha = 0.5;
-global snr_coeff = 0.9 % SNR est averaging filter coeff
+global snr_coeff;
+snr_coeff = 0.9;       % SNR est averaging filter coeff
+global Nph;
+Nph = 9;               % number of symbols to estimate phase over
+                       % must be odd number as we take centre symbol
 
 % root raised cosine (Root Nyquist) filter 
 
@@ -382,19 +386,56 @@ function [rx_symbols rx_timing env] = rx_est_timing(rx_filt, rx_baseband, nin)
 endfunction
 
 
-% Phase estimation function - probably won't work over a HF channel
-% Tries to operate over a single symbol but uses phase information from
-% all Nc carriers which should increase the SNR of phase estimate.
-% Maybe phase is coherent over a couple of symbols in HF channel,not
-% sure but it's worth 3dB so worth experimenting or using coherent as
-% an option.
+% Experimental "feed forward" phase estimation function - estimates
+% phase over a windows of Nph (e.g. Nph = 9) symbols.  May not work
+% well on HF channels but lets see.  Has a phase ambiguity of m(pi/4)
+% where m=0,1,2 which needs to be corrected outside of this function
 
-function rx_phase = rx_est_phase(prev_rx_symbols, rx_symbols)
+function [phase_offsets ferr] = rx_est_phase(rx_symbols)
+  global rx_symbols_mem;
+  global prev_phase_offsets;
+  global phase_amb;
+  global Nph;
+  global Nc;
 
-  % modulation strip
+  % keep record of Nph symbols
 
-  rx_phase = angle(sum(rx_symbols .^ 4))/4;
+  rx_symbols_mem(:,1:Nph-1) = rx_symbols_mem(:,2:Nph);
+  rx_symbols_mem(:,Nph) = rx_symbols;
  
+  % estimate and correct phase offset based of modulation stripped samples
+
+  phase_offsets = zeros(Nc+1,1);
+  for c=1:Nc+1
+
+    % rotate QPSK constellation to a single point
+    mod_stripped = abs(rx_symbols_mem(c,:)) .* exp(j*4*angle(rx_symbols_mem(c,:)));
+    
+    % find average phase offset, which will be on -pi/4 .. pi/4
+    sum_real = sum(real(mod_stripped));
+    sum_imag = sum(imag(mod_stripped));
+    phase_offsets(c) = atan2(sum_imag, sum_real)/4;
+
+    % determine if phase has jumped from - -> +    
+    if (prev_phase_offsets(c) < -pi/8) && (phase_offsets(c) > pi/8)
+      phase_amb(c) -= pi/2;
+      if (phase_amb(c) < -pi)
+        phase_amb(c) += 2*pi;
+      end
+    end
+    
+    % determine if phase has jumped from + -> -    
+    if (prev_phase_offsets(c) > pi/8) && (phase_offsets(c) < -pi/8)
+      phase_amb(c) += pi/2;
+      if (phase_amb(c) > pi)
+        phase_amb(c) -= 2*pi;
+      end
+    end
+  end
+
+  ferr = mean(phase_offsets - prev_phase_offsets);
+  prev_phase_offsets = phase_offsets;
+
 endfunction
 
 
@@ -415,6 +456,7 @@ function [rx_bits sync_bit f_err phase_difference] = qpsk_to_bits(prev_rx_symbol
     % map (Nc,1) DQPSK symbols back into an (1,Nc*Nb) array of bits
 
     for c=1:Nc
+      msb = lsb = 0;
       d = phase_difference(c);
       if ((real(d) >= 0) && (imag(d) >= 0))
          msb = 0; lsb = 0;
@@ -467,13 +509,15 @@ function [sig_est noise_est] = snr_update(sig_est, noise_est, phase_difference)
     % vector of mags, one for each carrier.
 
     s = abs(phase_difference);
-
+    
     % signal mag estimate for each carrier is a smoothed version
     % of instantaneous magntitude, this gives us a vector of smoothed
     % mag estimates, one for each carrier.
-
+    
     sig_est = snr_coeff*sig_est + (1 - snr_coeff)*s;
 
+    %printf("s: %f sig_est: %f snr_coeff: %f\n", s(1), sig_est(1), snr_coeff);
+
     % noise mag estimate is distance of current symbol from average
     % location of that symbol.  We reflect all symbols into the first
     % quadrant for convenience.
@@ -486,6 +530,7 @@ function [sig_est noise_est] = snr_update(sig_est, noise_est, phase_difference)
     % noise power estimates, one for each carrier.
 
     noise_est = snr_coeff*noise_est + (1 - snr_coeff)*n;
+
 endfunction
 
 
@@ -525,6 +570,12 @@ function bits = get_test_bits(nbits)
  
   for i=1:nbits
     bits(i) = test_bits(current_test_bit++);
+    %if (mod(i,2) == 0)
+    %  bits(i) = 1;
+    %else
+    %  bits(i) = 0;
+    %end
+    
     if (current_test_bit > Ntest_bits)
       current_test_bit = 1;
     endif
@@ -536,9 +587,8 @@ endfunction
 % Accepts nbits from rx and attempts to sync with test_bits sequence.
 % if sync OK measures bit errors
 
-function [sync bit_errors] = put_test_bits(rx_bits)
+function [sync bit_errors] = put_test_bits(test_bits, rx_bits)
   global Ntest_bits;       % length of test sequence
-  global test_bits;
   global rx_test_bits_mem;
 
   % Append to our memory
@@ -933,3 +983,11 @@ current_test_bit = 1;
 global rx_test_bits_mem;
 rx_test_bits_mem = zeros(1,Ntest_bits);
 
+% Experimental phase estimator states ----------------------
+
+global rx_symbols_mem;
+rx_symbols_mem = zeros(Nc+1, Nph);
+global prev_phase_offsets;
+prev_phase_offsets = zeros(Nc+1, 1);
+global phase_amb;
+phase_amb = zeros(Nc+1, 1);
diff --git a/codec2-dev/octave/fdmdv_demod.m b/codec2-dev/octave/fdmdv_demod.m
index 489a616c..0e7a4e80 100644
--- a/codec2-dev/octave/fdmdv_demod.m
+++ b/codec2-dev/octave/fdmdv_demod.m
@@ -115,7 +115,7 @@ function fdmdv_demod(rawfilename, nbits, pngname)
 
     % count bit errors if we find a test frame
 
-    [test_frame_sync bit_errors] = put_test_bits(rx_bits);
+    [test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits);
     if (test_frame_sync == 1)
       total_bit_errors = total_bit_errors + bit_errors;
       total_bits = total_bits + Ntest_bits;
diff --git a/codec2-dev/octave/fdmdv_demod_coh.m b/codec2-dev/octave/fdmdv_demod_coh.m
new file mode 100644
index 00000000..94b08802
--- /dev/null
+++ b/codec2-dev/octave/fdmdv_demod_coh.m
@@ -0,0 +1,253 @@
+% fdmdv_demod_coh.m
+%
+% Demodulator function for FDMDV modem (Octave version).  Requires
+% 8kHz sample rate raw files as input.  This version uses experimental
+% psuedo coherent demodulation.
+%
+% Copyright David Rowe 2013
+% This program is distributed under the terms of the GNU General Public License 
+% Version 2
+%
+
+function fdmdv_demod_coh(rawfilename, nbits, pngname)
+
+  fdmdv; % include modem code
+
+  modulation = 'dqpsk';
+
+  fin = fopen(rawfilename, "rb");
+  gain = 1000;
+  frames = nbits/(Nc*Nb);
+
+  prev_rx_symbols = ones(Nc+1,1);
+  foff_phase = 1;
+
+  % BER stats
+
+  total_bit_errors = 0;
+  total_bits = 0;
+  bit_errors_log = [];
+  sync_log = [];
+  test_frame_sync_log = [];
+  test_frame_sync_state = 0;
+
+  % SNR states
+
+  sig_est = zeros(Nc+1,1);
+  noise_est = zeros(Nc+1,1);
+
+  % logs of various states for plotting
+
+  rx_symbols_log = [];
+  rx_timing_log = [];
+  foff_log = [];
+  rx_fdm_log = [];
+  snr_est_log = [];
+
+  % misc states
+
+  nin = M; % timing correction for sample rate differences
+  foff = 0;
+  track_log = [];
+  track = 0;
+  fest_state = 0;
+
+  % psuedo coherent demod states
+
+  rx_symbols_ph_log = [];
+  prev_rx_symbols_ph = ones(Nc+1,1);
+  rx_phase_offsets_log = [];
+  phase_amb_log = [];
+
+  % Main loop ----------------------------------------------------
+
+  for f=1:frames
+    
+    % obtain nin samples of the test input signal
+    
+    for i=1:nin
+      rx_fdm(i) = fread(fin, 1, "short")/gain;
+    end
+    
+    rx_fdm_log = [rx_fdm_log rx_fdm(1:nin)];
+
+    % frequency offset estimation and correction
+
+    [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, nin);
+    [foff_coarse S1 S2] = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin);
+    
+    if track == 0
+      foff  = foff_coarse;
+    end
+    foff_log = [ foff_log foff ];
+    foff_rect = exp(j*2*pi*foff/Fs);
+
+    for i=1:nin
+      foff_phase *= foff_rect';
+      rx_fdm(i) = rx_fdm(i)*foff_phase;
+    end
+
+    % baseband processing
+
+    rx_baseband = fdm_downconvert(rx_fdm, nin);
+    rx_filt = rx_filter(rx_baseband, nin);
+
+    [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin);    
+    rx_timing_log = [rx_timing_log rx_timing];
+
+    nin = M;
+    if rx_timing > 2*M/P
+       nin += M/P;
+    end
+    if rx_timing < 0;
+       nin -= M/P;
+    end
+
+    rx_symbols_log = [rx_symbols_log rx_symbols.*(conj(prev_rx_symbols)./abs(prev_rx_symbols))*exp(j*pi/4)];
+
+    % coherent phase offset estimation ------------------------------------
+
+    [rx_phase_offsets ferr] = rx_est_phase(rx_symbols);
+    rx_phase_offsets_log = [rx_phase_offsets_log rx_phase_offsets];
+    phase_amb_log = [phase_amb_log phase_amb];
+    rx_symbols_ph = rx_symbols_mem(:,floor(Nph/2)+1) .* exp(-j*(rx_phase_offsets + phase_amb));
+    rx_symbols_ph_log = [rx_symbols_ph_log rx_symbols_ph .* exp(j*pi/4)];
+    rx_symbols_ph = -1 + 2*(real(rx_symbols_ph .* exp(j*pi/4)) > 0) + j*(-1 + 2*(imag(rx_symbols_ph .* exp(j*pi/4)) > 0));
+
+    % Std differential (used for freq offset est and BPSK sync) and psuedo coherent detection -----------------------
+
+    [rx_bits_unused sync        f_err       pd       ] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+    [rx_bits        sync_unused ferr_unused pd_unused] = qpsk_to_bits(prev_rx_symbols_ph, rx_symbols_ph, 'dqpsk');
+ 
+    foff -= 0.5*f_err;
+    prev_rx_symbols = rx_symbols;
+    prev_rx_symbols_ph = rx_symbols_ph;
+    sync_log = [sync_log sync];
+
+    [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
+    snr_est = calc_snr(sig_est, noise_est);
+    snr_est_log = [snr_est_log snr_est];
+
+    % freq est state machine
+
+    [track fest_state] = freq_state(sync, fest_state);
+    track_log = [track_log track];
+
+    % count bit errors if we find a test frame
+
+    [test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits);
+    if (test_frame_sync == 1)
+      total_bit_errors = total_bit_errors + bit_errors;
+      total_bits = total_bits + Ntest_bits;
+      bit_errors_log = [bit_errors_log bit_errors/Ntest_bits];
+    else
+      bit_errors_log = [bit_errors_log 0];
+    end
+
+    % test frame sync state machine, just for more informative plots
+    
+    next_test_frame_sync_state = test_frame_sync_state;
+    if (test_frame_sync_state == 0)
+      if (test_frame_sync == 1)      
+        next_test_frame_sync_state = 1;
+	test_frame_count = 0;
+      end
+    end
+
+    if (test_frame_sync_state == 1)
+      % we only expect another test_frame_sync pulse every 4 symbols
+      test_frame_count++;
+      if (test_frame_count == 4)
+        test_frame_count = 0;
+        if ((test_frame_sync == 0))      
+          next_test_frame_sync_state = 0;
+        end
+      end
+    end
+    test_frame_sync_state = next_test_frame_sync_state;
+    test_frame_sync_log = [test_frame_sync_log test_frame_sync_state];
+
+  end
+  
+  % ---------------------------------------------------------------------
+  % Print Stats
+  % ---------------------------------------------------------------------
+
+  ber = total_bit_errors / total_bits;
+
+  printf("%d bits  %d errors  BER: %1.4f\n",total_bits, total_bit_errors, ber);
+
+  % ---------------------------------------------------------------------
+  % Plots
+  % ---------------------------------------------------------------------
+
+  xt = (1:frames)/Rs;
+  secs = frames/Rs;
+
+  figure(1)
+  clf;
+  [n m] = size(rx_symbols_log);
+  plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+')
+  axis([-2 2 -2 2]);
+  title('Scatter Diagram');
+
+  figure(2)
+  clf;
+  subplot(211)
+  plot(xt, rx_timing_log)
+  title('timing offset (samples)');
+  subplot(212)
+  plot(xt, foff_log, '-;freq offset;')
+  hold on;
+  plot(xt, track_log*75, 'r;course-fine;');
+  hold off;
+  title('Freq offset (Hz)');
+  grid
+
+  figure(3)
+  clf;
+  subplot(211)
+  [a b] = size(rx_fdm_log);
+  xt1 = (1:b)/Fs;
+  plot(xt1, rx_fdm_log);
+  title('Rx FDM Signal');
+  subplot(212)
+  spec(rx_fdm_log,8000);
+  title('FDM Rx Spectrogram');
+
+  figure(4)
+  clf;
+  subplot(311)
+  stem(xt, sync_log)
+  axis([0 secs 0 1.5]);
+  title('BPSK Sync')
+  subplot(312)
+  stem(xt, bit_errors_log);
+  title('Bit Errors for test frames')
+  subplot(313)
+  plot(xt, test_frame_sync_log);
+  axis([0 secs 0 1.5]);
+  title('Test Frame Sync')
+
+  figure(5)
+  clf;
+  plot(xt, snr_est_log);
+  title('SNR Estimates')
+ 
+  figure(6)
+  clf;
+  [n m] = size(rx_symbols_ph_log);
+  plot(real(rx_symbols_ph_log(1:Nc+1,15:m)),imag(rx_symbols_ph_log(1:Nc+1,15:m)),'+')
+  %plot(real(rx_symbols_ph_log(2,15:m)),imag(rx_symbols_ph_log(2,15:m)),'+')
+  axis([-2 2 -2 2]);
+  title('Scatter Diagram - after phase correction');
+
+  figure(7)
+  clf;
+  subplot(211)
+  plot(rx_phase_offsets_log(1,:))
+  subplot(212)
+  plot(phase_amb_log(1,:))
+  title('Rx Phase Offset Est')
+
+endfunction
diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m
index a78eedfe..147322df 100644
--- a/codec2-dev/octave/fdmdv_ut.m
+++ b/codec2-dev/octave/fdmdv_ut.m
@@ -12,7 +12,7 @@ fdmdv;               % load modem code
  
 % Simulation Parameters --------------------------------------
 
-frames = 25;
+frames = 100;
 EbNo_dB = 7.3;
 Foff_hz = 0;
 modulation = 'dqpsk';
@@ -96,6 +96,7 @@ fest_state = 0;
 track = 0;
 track_log = [];
 
+snr_log = [];
 
 % ---------------------------------------------------------------------
 % Main loop 
@@ -206,6 +207,7 @@ for f=1:frames
   % Update SNR est
 
   [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
+  snr_log = [snr_log calc_snr(sig_est, noise_est)];
 
   % count bit errors if we find a test frame
   % Allow 15 frames for filter memories to fill and time est to settle
@@ -316,3 +318,6 @@ plot(test_frame_sync_log);
 axis([0 frames 0 1.5]);
 title('Test Frame Sync')
 
+figure(5)
+clf
+plot(snr_log)
diff --git a/codec2-dev/octave/fdmdv_ut_coh.m b/codec2-dev/octave/fdmdv_ut_coh.m
index 7450b2fd..a1598878 100644
--- a/codec2-dev/octave/fdmdv_ut_coh.m
+++ b/codec2-dev/octave/fdmdv_ut_coh.m
@@ -19,9 +19,9 @@ fdmdv;               % load modem code
  
 % Simulation Parameters --------------------------------------
 
-frames = 100;
-EbNo_dB = 73;
-Foff_hz = 0;
+frames = 200;
+EbNo_dB = 7;
+Foff_hz = -100;
 hpa_clip = 150;
 
 % ------------------------------------------------------------
@@ -104,37 +104,10 @@ track_log = [];
 
 snr_log = [];
 
-% PSK demod ----
-
-% Set up matrix of tx symbols, one row for each carrier.  
-
-% The differential PSK encoder is driven by a sequence of test bits.
-% This test sequence is 4 frames (122 bits long), which means it
-% repeats every 4 frames.  However the DPSK encoder also has 4 states,
-% so it turns out that the transmitted sequence of PSK symbols repeats
-% itself every 16 frames.  Depending on the test data some carriers
-% repeat more often than that, but the sequence as a whole repeats
-% every 16 frames.
-
-% Each row in the matrix is the sequence of symbols for one carrier.
-% This repeats itself every 16 symbols.  We would like to demodulate
-% and measure the BER of this block of 16 symbols.  To do this we need
-% to correct any phase offset.
-
-prev_psk_test_symbols = ones(Nc+1,1);
-tx_symbols_mem = [];
-for i=1:6
-  psk_test_symbols = generate_psk_test_symbols(prev_psk_test_symbols);
-  prev_psk_test_symbols = psk_test_symbols(:,4);
-  tx_symbols_mem = [tx_symbols_mem psk_test_symbols .* exp(j*pi/4)];
-end
-
-% PSK BER counters
-
-psk_total_bits = 0;
-psk_total_bit_errors = 0;
-
 rx_symbols_ph_log = [];
+prev_rx_symbols_ph = ones(Nc+1,1);
+rx_phase_offsets_log = [];
+phase_amb_log = [];
 
 % ---------------------------------------------------------------------
 % Main loop 
@@ -206,7 +179,7 @@ for f=1:frames
     foff = foff_course;
   end
 
-  foff = 0; % disable for now
+  %foff = 0; % disable for now
 
   foff_log = [ foff_log foff ];
   foff_rect = exp(j*2*pi*foff/Fs);
@@ -223,13 +196,28 @@ for f=1:frames
   rx_filt = rx_filter(rx_baseband, M);
 
   [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, M);
+  rx_symbols_log = [rx_symbols_log rx_symbols.*(conj(prev_rx_symbols)./abs(prev_rx_symbols))*exp(j*pi/4)];
   rx_timing_log = [rx_timing_log rx_timing];
 
-  [rx_bits sync foff_fine pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
-  rx_symbols_log = [rx_symbols_log rx_symbols .* exp(j*pi/4)];
+  % coherent phase offset estimation ------------------------------------
+
+  [rx_phase_offsets ferr] = rx_est_phase(rx_symbols);
+  rx_phase_offsets_log = [rx_phase_offsets_log rx_phase_offsets];
+  phase_amb_log = [phase_amb_log phase_amb];
+  rx_symbols_ph = rx_symbols_mem(:,floor(Nph/2)+1) .* exp(-j*(rx_phase_offsets + phase_amb));
+  rx_symbols_ph_log = [rx_symbols_ph_log rx_symbols_ph .* exp(j*pi/4)];
+  rx_symbols_ph = -1 + 2*(real(rx_symbols_ph .* exp(j*pi/4)) > 0) + j*(-1 + 2*(imag(rx_symbols_ph .* exp(j*pi/4)) > 0));
+
+  % Std differential (used for freq offset est and BPSK sync) and psuedo coherent detection -----------------------
+
+  [rx_bits_unused sync        ferr        pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
+  [rx_bits        sync_unused ferr_unused pd] = qpsk_to_bits(prev_rx_symbols_ph, rx_symbols_ph, 'dqpsk');
+
+  %----------------------------------------------------------------------
 
   foff -= 0.5*ferr;
   prev_rx_symbols = rx_symbols;
+  prev_rx_symbols_ph = rx_symbols_ph;
   sync_log = [sync_log sync];
   
   % freq est state machine
@@ -242,48 +230,11 @@ for f=1:frames
   [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
   snr_log = [snr_log calc_snr(sig_est, noise_est)];
 
-  % coherent PSK demod --------------------------------------
-
-  % update memory of the last Nph rx symbols
-
-  rx_symbols_mem(:,1:Npom-1) = rx_symbols_mem(:,2:Npom);
-  rx_symbols_mem(:,Npom) = rx_symbols .* exp(j*pi/4);
-
-  % estimate and correct phase of each received symbol in 16 symbol frame
-
-  for c=1:Nc                  
-    for s=Nph2+1:Nph2+Nstates  % symbols in middle of matrix that we correct 
-      st = s - Nph2; en = s + Nph2;
-      corr = rx_symbols_mem(c,st:en) * tx_symbols_mem(c,st:en)';
-      rx_symbols_mem_ph(c,s) = rx_symbols_mem(c,s) .* exp(-j*angle(corr));
-    end
-  end
-  
-  % count bit errors in centre of window
-
-  psk_bit_errors = 0;
-  psk_bits = Nc*Nstates;
-  st = Nph2+1; en = Nph2+Nstates;
-  for c=1:Nc
-    psk_bit_errors += sum(real(rx_symbols_mem_ph(c,st:en)) .* real(tx_symbols_mem(c,st:en)) < 0);
-    psk_bit_errors += sum(imag(rx_symbols_mem_ph(c,st:en)) .* imag(tx_symbols_mem(c,st:en)) < 0);
-  end
-  printf("psk_bit_errors: %d corr: %f\n", psk_bit_errors, angle(corr));
-  ber = psk_bit_errors/psk_bits;
-  if ber < 0.1
-    psk_total_bit_errors += psk_bit_errors;
-    psk_total_bits += psk_bits;
-    % phase correction only valid when we are in sync
-    rx_symbols_ph_log = [rx_symbols_ph_log rx_symbols_mem_ph(:,st:en)];
-  end
-
-  % ---------------------------------------------------------
-
   % count bit errors if we find a test frame
 
   [test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits);
 
-  if test_frame_sync == 1
+  if (test_frame_sync == 1) && (f > 15)
     total_bit_errors = total_bit_errors + bit_errors;
     total_bits = total_bits + Ntest_bits;
     bit_errors_log = [bit_errors_log bit_errors];
@@ -334,10 +285,6 @@ printf("\nDPSK\n");
 printf("  bits......: %d\n", total_bits);
 printf("  errors....: %d\n", total_bit_errors);
 printf("  BER.......: %1.4f\n",  ber);
-printf("PSK\n");
-printf("  bits......: %d\n", psk_total_bits);
-printf("  errors....: %d\n", psk_total_bit_errors);
-printf("  BER.......: %1.4f\n",  psk_total_bit_errors/psk_total_bits);
 
 % ---------------------------------------------------------------------
 % Plots
@@ -365,18 +312,6 @@ title('Freq offset (Hz)');
 
 figure(3)
 clf;
-subplot(211)
-plot(real(tx_fdm_log));
-title('FDM Tx Signal');
-subplot(212)
-Nfft=Fs;
-S=fft(rx_fdm_log,Nfft);
-SdB=20*log10(abs(S));
-plot(SdB(1:Fs/4))
-title('FDM Rx Spectrum');
-
-figure(4)
-clf;
 subplot(311)
 stem(sync_log)
 axis([0 frames 0 1.5]);
@@ -389,24 +324,18 @@ plot(test_frame_sync_log);
 axis([0 frames 0 1.5]);
 title('Test Frame Sync')
 
-figure(5)
-clf;
-subplot(211)
-stem(real(rx_symbols_mem(2,:)))
-subplot(212)
-stem(imag(rx_symbols_mem(2,:)))
-
-figure(6)
-clf;
-subplot(211)
-stem(real(tx_symbols_mem(2,:)))
-subplot(212)
-stem(imag(tx_symbols_mem(2,:)))
-
-figure(7)
+figure(4)
 clf;
 [n m] = size(rx_symbols_ph_log);
 plot(real(rx_symbols_ph_log(1:Nc+1,15:m)),imag(rx_symbols_ph_log(1:Nc+1,15:m)),'+')
 %plot(real(rx_symbols_ph_log(2,15:m)),imag(rx_symbols_ph_log(2,15:m)),'+')
 axis([-3 3 -3 3]);
 title('Scatter Diagram - after phase correction');
+
+figure(5)
+clf;
+subplot(211)
+plot(rx_phase_offsets_log(1,:))
+subplot(212)
+plot(phase_amb_log(1,:))
+title('Rx Phase Offset Est')
-- 
2.25.1