working two state freq est/tracking, sample clock correction logic, tested over paths...
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 13 Apr 2012 01:20:40 +0000 (01:20 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 13 Apr 2012 01:20:40 +0000 (01:20 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@366 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/README_fdmdv.txt
codec2-dev/octave/fdmdv.m
codec2-dev/octave/fdmdv_demod.m
codec2-dev/octave/fdmdv_ut.m

index 57d29c5e9bd0d8fe06ec40fcef0afacbb4747762..a22b65a5cf61f4ce88cc4e2962ced4a3b2e4d094 100644 (file)
@@ -18,6 +18,7 @@ TODO
     + does modem fall over?
 [ ] try non-flat channel, e.g. 3dB difference between hi and low tones
     + make sure all estimators keep working
+[ ] test rx level sensitivity, i.e. 0 to 20dB attenuation
 [ ] try to run from shell script
 [ ] run a few tests
 [ ] start coding in C and repeat tests
index 444fef8e04cbc5acbd959c2e29a66b3ce8f78702..7f65353c83563e1d7556b06e998c7fc429dc113c 100644 (file)
@@ -82,8 +82,6 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
 
   if (strcmp(modulation,'dqpsk')) 
  
-
-  if 1
     % map to (Nc,1) DQPSK symbols
 
     for c=1:Nc
@@ -102,21 +100,6 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
          tx_symbols(c) = -j*prev_tx_symbols(c);
       endif 
     end
-  else
-    % map to pi/4 DQPSK from spra341 Eq (6) & (7)
-
-    for c=1:Nc
-
-      msb = tx_bits_matrix(c,1); lsb = tx_bits_matrix(c,2);
-      a = 2*msb - 1;
-      b = 2*lsb - 1;
-      p = prev_tx_symbols(c);
-      inphase    = (real(p)*a - imag(p)*b)*0.707;
-      quadrature = (imag(p)*a + real(p)*b)*0.707;
-      tx_symbols(c) = inphase + j*quadrature;
-    end
-  end
-
   else
     % QPSK mapping
     tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2);
@@ -215,7 +198,7 @@ endfunction
 
 % Frequency shift each modem carrier down to Nc baseband signals
 
-function rx_baseband = fdm_downconvert(rx_fdm)
+function rx_baseband = fdm_downconvert(rx_fdm, nin)
   global Fs;
   global M;
   global Nc;
@@ -223,12 +206,12 @@ function rx_baseband = fdm_downconvert(rx_fdm)
   global phase_rx;
   global freq;
 
-  rx_baseband = zeros(1,M);
+  rx_baseband = zeros(1,nin);
 
   % Nc/2 tones below centre freq
   
   for c=1:Nc/2
-      for i=1:M
+      for i=1:nin
         phase_rx(c) = phase_rx(c) * freq(c);
        rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
       end
@@ -237,7 +220,7 @@ function rx_baseband = fdm_downconvert(rx_fdm)
   % Nc/2 tones above centre freq  
 
   for c=Nc/2+1:Nc
-      for i=1:M
+      for i=1:nin
         phase_rx(c) = phase_rx(c) * freq(c);
        rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
       end
@@ -246,7 +229,7 @@ function rx_baseband = fdm_downconvert(rx_fdm)
   % Pilot
 
   c = Nc+1;
-  for i=1:M
+  for i=1:nin
     phase_rx(c) = phase_rx(c) * freq(c);
     rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
   end
@@ -256,7 +239,7 @@ endfunction
 
 % Receive filter each baseband signal at oversample rate P
 
-function rx_filt = rx_filter(rx_baseband)
+function rx_filt = rx_filter(rx_baseband, nin)
   global Nc;
   global M;
   global P;
@@ -265,14 +248,14 @@ function rx_filt = rx_filter(rx_baseband)
   global gt_alpha5_root;
   global Fsep;
 
-  rx_filt = zeros(Nc+1,P);
+  rx_filt = zeros(Nc+1,nin*P/M);
 
   % rx filter each symbol, generate P filtered output samples for each symbol.
   % Note we keep memory at rate M, it's just the filter output at rate P
 
   N=M/P;
   j=1;
-  for i=1:N:M
+  for i=1:N:nin
     rx_filter_memory(:,Nfilter-N+1:Nfilter) = rx_baseband(:,i:i-1+N);
     rx_filt(:,j) = rx_filter_memory * gt_alpha5_root';
     rx_filter_memory(:,1:Nfilter-N) = rx_filter_memory(:,1+N:Nfilter);
@@ -284,7 +267,7 @@ endfunction
 % Estimate frequency offset of FDM signal using BPSK pilot.  This is quite
 % sensitive to pilot tone level wrt other carriers
 
-function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev)
+function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev, nin)
   global M;
   global Npilotbaseband;
   global pilot_baseband1;
@@ -295,15 +278,15 @@ function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev)
   % down convert latest M samples of pilot by multiplying by
   % ideal BPSK pilot signal we have generated locally
  
-  pilot_baseband1(1:Npilotbaseband-M) = pilot_baseband1(M+1:Npilotbaseband);
-  pilot_baseband2(1:Npilotbaseband-M) = pilot_baseband2(M+1:Npilotbaseband);
-  for i=1:M
-    pilot_baseband1(Npilotbaseband-M+i) = rx_fdm(i) * conj(pilot(i)); 
-    pilot_baseband2(Npilotbaseband-M+i) = rx_fdm(i) * conj(pilot_prev(i)); 
+  pilot_baseband1(1:Npilotbaseband-nin) = pilot_baseband1(nin+1:Npilotbaseband);
+  pilot_baseband2(1:Npilotbaseband-nin) = pilot_baseband2(nin+1:Npilotbaseband);
+  for i=1:nin
+    pilot_baseband1(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot(i)); 
+    pilot_baseband2(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot_prev(i)); 
   end
 
-  [foff1 max1 pilot_lpf1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1);
-  [foff2 max2 pilot_lpf2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2);
+  [foff1 max1 pilot_lpf1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1, nin);
+  [foff2 max2 pilot_lpf2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2, nin);
 
   if max1 > max2
     foff = foff1;
@@ -315,7 +298,7 @@ endfunction
 
 % LPF and peak pick part of freq est, put in a function as we call it twice
 
-function  [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf)
+function  [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin)
   global M;
   global Npilotlpf;
   global Npilotcoeff;
@@ -325,9 +308,9 @@ function  [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf)
 
   % LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset
 
-  pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf);
+  pilot_lpf(1:Npilotlpf-nin) = pilot_lpf(nin+1:Npilotlpf);
   j = 1;
-  for i = Npilotlpf-M+1:Npilotlpf
+  for i = Npilotlpf-nin+1:Npilotlpf
     pilot_lpf(i) = pilot_baseband(j:j+Npilotcoeff) * pilot_coeff';
     j++;
   end
@@ -354,7 +337,7 @@ endfunction
 
 % Estimate optimum timing offset, and symbol receive symbols
 
-function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
+function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin)
   global M;
   global Nt;
   global Nc;
@@ -365,10 +348,22 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
   global Nfiltertiming;
   global gt_alpha5_root;
 
+  % nin  adjust 
+  % --------------------------------
+  % 120  -1 (one less rate P sample)
+  % 160   0 (nominal)
+  % 200   1 (one more rate P sample)
+
+  adjust = P - nin*P/M;
+
   % update buffer of Nt rate P filtered symbols
 
-  rx_filter_mem_timing(:,1:(Nt-1)*P) = rx_filter_mem_timing(:,P+1:Nt*P);
-  rx_filter_mem_timing(:,(Nt-1)*P+1:Nt*P) = rx_filt(1:Nc,:);
+  rx_filter_mem_timing(:,1:(Nt-1)*P+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P);
+  %size((Nt-1)*P+1+adjust:Nt*P)
+  %size(rx_filt)
+  %adjust
+  %nin
+  rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(1:Nc,:);
 
   % sum envelopes of all carriers
 
@@ -396,8 +391,8 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
   % baseband signal at rate M this enables us to resample the filtered
   % rx symbol with M sample precision once we have rx_timing
 
-  rx_baseband_mem_timing(:,1:Nfiltertiming-M) = rx_baseband_mem_timing(:,M+1:Nfiltertiming);
-  rx_baseband_mem_timing(:,Nfiltertiming-M+1:Nfiltertiming) = rx_baseband;
+  rx_baseband_mem_timing(:,1:Nfiltertiming-nin) = rx_baseband_mem_timing(:,nin+1:Nfiltertiming);
+  rx_baseband_mem_timing(:,Nfiltertiming-nin+1:Nfiltertiming) = rx_baseband;
 
   % sample right in the middle of the timing estimator window, by filtering
   % at rate M
@@ -426,7 +421,7 @@ endfunction
 
 % convert symbols back to an array of bits
 
-function [rx_bits sync_bit] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
+function [rx_bits sync_bit f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
   global Nc;
   global Nb;
   global Nb;
@@ -462,8 +457,10 @@ function [rx_bits sync_bit] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulati
     phase_difference(Nc+1) = rx_symbols(Nc+1) .* conj(prev_rx_symbols(Nc+1));
     if (real(phase_difference(Nc+1)) < 0)
       sync_bit = 0;
+      f_err = imag(phase_difference(Nc+1));
     else
       sync_bit = 1;
+      f_err = -imag(phase_difference(Nc+1));
     end
 
   else
@@ -709,3 +706,82 @@ function [buf_out t nin] = resample(buf_in, t, ratio, nout)
   t -= delta;
 
 endfunction
+
+
+% freq offset state machine.  Moves between acquire and track states based
+% on BPSK pilot sequence.  Freq offset estimator occasionally makes mistakes
+% when used continuously.  So we use it unit we have acquired the BPSK pilot,
+% then switch to a more robust tracking algorithm.  If we lose sync we switch
+% back to acquire mode for fast-requisition.
+
+function [track state] = freq_state(sync_bit, state)
+
+  % acquire state, look for 6 symbol 010101 sequence from sync bit
+
+  next_state = state;
+  if state == 0
+    if sync_bit == 0
+      next_state = 1;
+    end        
+  end
+  if state == 1
+    if sync_bit == 1
+      next_state = 2;
+    else 
+      next_state = 0;
+    end        
+  end
+  if state == 2
+    if sync_bit == 0
+      next_state = 3;
+    else 
+      next_state = 0;
+    end        
+  end
+  if state == 3
+    if sync_bit == 1
+      next_state = 4;
+    else 
+      next_state = 0;
+    end        
+  end
+  if state == 4
+    if sync_bit == 0
+      next_state = 5;
+    else 
+      next_state = 0;
+    end        
+  end
+  if state == 5
+    if sync_bit == 1
+      next_state = 6;
+    else 
+      next_state = 0;
+    end        
+  end
+
+  % states 6 and above are track mode, make sure we keep getting 0101 sync bit sequence
+
+  if state == 6
+    if sync_bit == 0
+      next_state = 7;
+    else 
+      next_state = 0;
+    end        
+  end
+  if state == 7
+    if sync_bit == 1
+      next_state = 6;
+    else 
+      next_state = 0;
+    end        
+  end
+
+  state = next_state;
+  if state >= 6
+    track = 1;
+  else
+    track = 0;
+  end
+endfunction
+
index 0d52c5cf54b67bcdc64daa6617b43f744729be29..a5baa49c8d4c533779e88dc8affa4d163536f3db 100644 (file)
@@ -1,6 +1,6 @@
 % fdmdv_demod.m
 %
-% Demodulator function for FDMDV modem.  Requires 48kHz sample rate raw files
+% Demodulator function for FDMDV modem.  Requires 8kHz sample rate raw files
 % as input
 %
 % Copyright David Rowe 2012
@@ -25,7 +25,7 @@ function fdmdv_demod(rawfilename, nbits)
   
   pilot_lut = generate_pilot_lut;
   pilot_lut_index = 1;
-  prev_pilot = zeros(1,M);
+  prev_pilot_lut_index = 3*M+1;
 
   % BER stats
 
@@ -39,101 +39,93 @@ function fdmdv_demod(rawfilename, nbits)
   rx_symbols_log = [];
   rx_timing_log = [];
   foff_log = [];
+  rx_fdm_log = [];
 
-  % resampler states
+  % misc states
 
-  t = 3;
-  ratio = 1.000;
-  F=6;
-  MF=M*F;
-  nin = MF;
-  nin_size = MF+6;
-  buf_in = zeros(1,nin_size);
-  rx_fdm_buf = [];
-  ratio_log = [];
+  nin = M; % timing correction for sample rate differences
+  foff = 0;
+  track_log = [];
+  track = 1;
+  fest_state = 0;
 
   % Main loop ----------------------------------------------------
 
   for f=1:frames
-    % update buf_in memory
-
-    m = nin_size - nin;
-    for i=1:m
-      buf_in(i) = buf_in(i+nin);  
-    end
-    
     % obtain nin samples of the test input signal
-
-    for i=m+1:nin_size
-      buf_in(i) = fread(fin, 1, "short")/gain; 
+    
+    for i=1:nin
+      rx_fdm(i) = fread(fin, 1, "short")/gain; 
     end
 
-    % resample at 48kHz and decimate to 8kHz
-
-    [rx_fdm_mf t nin] = resample(buf_in, t, 1.0, MF);
-    rx_fdm = rx_fdm_mf(1:F:MF);
-
-    %rx_fdm = buf_in(m+1:m+n);
-
-    %for i=1:M
-    %  rx_fdm(i) = fread(fin, 1, "short")/gain; 
-    %end
-    rx_fdm_buf = [rx_fdm_buf rx_fdm];
+    rx_fdm_log = [rx_fdm_log rx_fdm(1:nin)];
 
     % frequency offset estimation and correction
 
-    for i=1:M
+    for i=1:nin
       pilot(i) = pilot_lut(pilot_lut_index);
       pilot_lut_index++;
       if pilot_lut_index > 4*M
         pilot_lut_index = 1;
       end
+      prev_pilot(i) = pilot_lut(prev_pilot_lut_index);
+      prev_pilot_lut_index++;
+      if prev_pilot_lut_index > 4*M
+        prev_pilot_lut_index = 1;
+      end
+    end
+    foff_coarse = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin);
+    if track == 0
+      foff  = foff_coarse;
     end
-    foff = rx_est_freq_offset(rx_fdm, pilot, prev_pilot);
-    prev_pilot = pilot;
     foff_log = [ foff_log foff ];
-    foff = 0;
     foff_rect = exp(j*2*pi*foff/Fs);
 
-    for i=1:M
+    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);
-    rx_filt = rx_filter(rx_baseband);
+    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);
+    [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin);
     rx_timing_log = [rx_timing_log rx_timing];
-    beta = 1E-7;
-    ratio += beta*rx_timing;
-    if (ratio > 1.002)
-       ratio = 1.002;
+    nin = M;
+    if rx_timing > 2*M/P
+       nin += M/P;
     end
-    if (ratio < 0.998)
-       ratio = 0.998;
+    if rx_timing < 0;
+       nin -= M/P;
     end
-    ratio_log = [ratio_log ratio];
 
     if strcmp(modulation,'dqpsk')
       rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)];
     else
       rx_symbols_log = [rx_symbols_log rx_symbols];
     endif
-    [rx_bits sync] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+    [rx_bits sync f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+    foff -= 0.5*f_err;
     prev_rx_symbols = rx_symbols;
     sync_log = [sync_log sync];
 
+    % 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
     % Allow 15 frames for filter memories to fill and time est to settle
 
     [test_frame_sync bit_errors] = put_test_bits(rx_bits);
-    if (test_frame_sync == 1 && f > 15)
+    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];
+      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
@@ -173,50 +165,56 @@ function fdmdv_demod(rawfilename, nbits)
   % 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,10:m)),imag(rx_symbols_log(1:Nc+1,10:m)),'+')
+  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(rx_timing_log(15:m))
+  plot(xt, rx_timing_log)
   title('timing offset (samples)');
   subplot(212)
-  plot(foff_log)
+  plot(xt, foff_log)
   title('Freq offset (Hz)');
   grid
 
   figure(3)
   clf;
   subplot(211)
-  %plot(rx_fdm_buf);
-  %title('FDM Rx Signal');
-  plot(ratio_log-1);
-  title('Sampling Clock error (ppm)');
+  [a b] = size(rx_fdm_log);
+  xt1 = (1:b)/Fs;
+  plot(xt1, rx_fdm_log);
+  title('Rx FDM Signal');
   subplot(212)
   Nfft=Fs;
-  S=fft(rx_fdm_buf,Nfft);
+  S=fft(rx_fdm_log,Nfft);
   SdB=20*log10(abs(S));
   plot(SdB(1:Fs/4))
-  title('FDM Tx Spectrum');
+  title('FDM Rx Spectrum');
 
   figure(4)
   clf;
   subplot(311)
-  stem(sync_log)
-  axis([0 frames 0 1.5]);
+  stem(xt, sync_log)
+  axis([0 secs 0 1.5]);
   title('BPSK Sync')
   subplot(312)
-  stem(bit_errors_log);
+  stem(xt, bit_errors_log);
   title('Bit Errors for test frames')
   subplot(313)
-  plot(test_frame_sync_log);
-  axis([0 frames 0 1.5]);
+  plot(xt, test_frame_sync_log);
+  axis([0 secs 0 1.5]);
   title('Test Frame Sync')
 
-  mean(ratio_log)
+  figure(5)
+  clf;
+  plot(xt, track_log);
+  axis([0 secs 0 1.5]);
 endfunction
index 65cb0d54dbdb93f9e282c5b6ceb2f3efdf6c5220..ad6df497bf9fc6756915c5397d6c85abffec81ca 100644 (file)
@@ -11,9 +11,9 @@ fdmdv;               % load modem code
  
 % Simulation Parameters --------------------------------------
 
-frames = 50;
+frames = 100;
 EbNo_dB = 7.3;
-Foff_hz = 0;
+Foff_hz = 10;
 modulation = 'dqpsk';
 hpa_clip = 150;
 
@@ -30,6 +30,8 @@ rx_baseband_log = [];
 rx_bits_offset = zeros(Nc*Nb*2);
 prev_tx_symbols = ones(Nc+1,1);
 prev_rx_symbols = ones(Nc+1,1);
+f_err = 0;
+foff = 0;
 foff_log = [];
 tx_baseband_log = [];
 tx_fdm_log = [];
@@ -47,6 +49,7 @@ test_frame_sync_state = 0;
 
 pilot_lut = generate_pilot_lut;
 pilot_lut_index = 1;
+prev_pilot_lut_index = 3*M+1;
 
 % fixed delay simuation
 
@@ -111,11 +114,12 @@ for f=1:frames
 
   % frequency offset
 
+  %Foff_hz += 1/Rs;
+  Foff = Foff_hz;
   for i=1:M
     % Time varying freq offset
-    % Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs));
-    % t++;
-    Foff = Foff_hz;
+    %Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs));
+    %t++;
     freq_offset = exp(j*2*pi*Foff/Fs);
     phase_offset *= freq_offset;
     tx_fdm(i) = phase_offset*tx_fdm(i);
@@ -155,9 +159,13 @@ for f=1:frames
     if pilot_lut_index > 4*M
       pilot_lut_index = 1;
     end
+    prev_pilot(i) = pilot_lut(prev_pilot_lut_index);
+    prev_pilot_lut_index++;
+    if prev_pilot_lut_index > 4*M
+      prev_pilot_lut_index = 1;
+    end
   end
-  foff = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot);
-  prev_pilot = pilot;
+  %foff = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M);
   foff_log = [ foff_log foff ];
   %foff = 0;
   foff_rect = exp(j*2*pi*foff/Fs);
@@ -169,11 +177,11 @@ for f=1:frames
 
   % baseband processing
 
-  rx_baseband = fdm_downconvert(rx_fdm_delay(1:M));
+  rx_baseband = fdm_downconvert(rx_fdm_delay(1:M), M);
   rx_baseband_log = [rx_baseband_log rx_baseband];
-  rx_filt = rx_filter(rx_baseband);
+  rx_filt = rx_filter(rx_baseband, M);
 
-  [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband);
+  [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, M);
   rx_timing_log = [rx_timing_log rx_timing];
 
   %rx_phase = rx_est_phase(rx_symbols);
@@ -185,10 +193,11 @@ for f=1:frames
   else
     rx_symbols_log = [rx_symbols_log rx_symbols];
   endif
-  [rx_bits sync] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+  [rx_bits sync f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+  foff -= 0.5*f_err;
   prev_rx_symbols = rx_symbols;
   sync_log = [sync_log sync];
-
+  
   % count bit errors if we find a test frame
   % Allow 15 frames for filter memories to fill and time est to settle
 
@@ -246,16 +255,16 @@ 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]);
+axis([-3 3 -3 3]);
 title('Scatter Diagram');
 
 figure(2)
 clf;
 subplot(211)
-plot(rx_timing_log)
+plot(rx_timing_log(15:m))
 title('timing offset (samples)');
 subplot(212)
-plot(foff_log)
+plot(foff_log(15:m))
 title('Freq offset (Hz)');
 
 figure(3)