modified freq offset estimation to fix bug that occurs with 1 symbol delay. Improves...

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

diff --git a/codec2-dev/octave/README_fdmdv.txt b/codec2-dev/octave/README_fdmdv.txt
index ecf224457dc61e2e45fcdfb9abbe940ad2381b14..d6f415ced947e9ea7d5085a7e9425b6f09ad0a86 100644 (file)
--- a/codec2-dev/octave/README_fdmdv.txt
+++ b/codec2-dev/octave/README_fdmdv.txt
@@ -6,11 +6,23 @@ sox -r 8000 -s -2 mod_dqpsk.raw -s -2 mod_dqpsk_8008hz.raw rate -h 8008
 
 TODO
 
-[ ] Get file based mod and demod working again
-[ ] write file based ch simulator
+[X] Get file based mod and demod working again
+[ ] timing wraps around
+    + what is affect of bouncing back and forth over boundary?
+    + could mean we are going back and forth a symbol
+    + do we need to logic to lose or gain a frame?
+    + think so, add or remove samples, or a whole frame
 [ ] demod outputs ber (maybe after settling time)
 [ ] try to run from shell script
 [ ] run a few tests
 [ ] start coding in C and repeat tests
 [ ] arb bit stream input to Octave mod and demod
 
+Tests
+
+[ ] fdmdv_demod('mod_dqpsk_8008hz.raw',1400*60);
+[ ] fdmdv_demod('mod_dqpsk_7992hz.raw',1400*60);
+[ ] mod_dqpsk_awgn_4dB_8008hz.raw
+[ ] mod_dqpsk_good_4dB_8008hz.raw
+[ ] mod_dqpsk_moderate_4dB_8008hz.raw
+[ ] mod_dqpsk_moderate_4dB_7992hz.raw

diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m
index af8b8ee1e632927800044c98516094a6de6c1db6..35017d5d0928e80e6fd41c34f28380aa43235dfa 100644 (file)
--- a/codec2-dev/octave/fdmdv.m
+++ b/codec2-dev/octave/fdmdv.m
@@ -1,7 +1,7 @@
 % fdmdv.m
 %
 % Functions that implement a Frequency Divison Multiplexed Modem for
-% Digital Voice (FMDV)over HF channels.
+% Digital Voice (FDMDV) over HF channels.
 %
 % Copyright David Rowe 2012
 % This program is distributed under the terms of the GNU General Public License 
@@ -284,30 +284,45 @@ 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)
+function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev)
   global M;
-  global Nc;
-  global Fs;
-  global Rs;
-  global Fcentre;
-  global freq;
-  global freq_rx_pilot;
   global Npilotbaseband;
-  global Npilotlpf;
-  global Npilotcoeff;
-  global pilot_baseband;
-  global pilot_lpf;
-  global pilot_coeff;
+  global pilot_baseband1;
+  global pilot_baseband2;
+  global pilot_lpf1;
+  global pilot_lpf2;
 
   % down convert latest M samples of pilot by multiplying by
-  % ideal two-tone BPSK pilot signal
+  % ideal BPSK pilot signal we have generated locally
  
-  pilot_baseband(1:Npilotbaseband-M) = pilot_baseband(M+1:Npilotbaseband);
-  c = Nc+1;
+  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_baseband(Npilotbaseband-M+i) = rx_fdm(i)*conj(pilot(i)); 
+    pilot_baseband1(Npilotbaseband-M+i) = rx_fdm(i) * conj(pilot(i)); 
+    pilot_baseband2(Npilotbaseband-M+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);
+
+  if max1 > max2
+    foff = foff1;
+  else
+    foff = foff2;
+  end  
+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)
+  global M;
+  global Npilotlpf;
+  global Npilotcoeff;
+  global Fs;
+  global Mpilotfft;
+  global pilot_coeff;
+
   % LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset
 
   pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf);
@@ -320,20 +335,12 @@ function foff = rx_est_freq_offset(rx_fdm, pilot)
   % decimate to improve DFT resolution, window and DFT
 
   Mpilot = Fs/(2*200);  % calc decimation rate given new sample rate is twice LPF freq
-  Mpilotfft = 256;
   s = pilot_lpf(1:Mpilot:Npilotlpf) .* hanning(Npilotlpf/Mpilot)';
   S = abs(fft(s, Mpilotfft));
 
-  %figure(3)
-  %plot(real(pilot_baseband))
-  %plot(real(s))
-  %figure(4)
-  %plot(abs(fft(pilot_baseband)))
-  %plot(S)
-
   % peak pick and convert to Hz
 
-  [x ix] = max(S);
+  [imax ix] = max(S);
   r = 2*200/Mpilotfft;     % maps FFT bin to frequency in Hz
   
   if ix > Mpilotfft/2
@@ -552,6 +559,7 @@ phase_rx = ones(Nc+1,1);
 
 % Freq offset estimator constants
 
+global Mpilotfft      = 256;
 global Npilotcoeff    = 30;                             % number of pilot LPF coeffs
 global pilot_coeff    = fir1(Npilotcoeff, 200/(Fs/2))'; % 200Hz LPF
 global Npilotbaseband = Npilotcoeff + 4*M;              % number of pilot baseband samples reqd for pilot LPF
@@ -559,10 +567,14 @@ global Npilotlpf      = 4*M;                            % number of samples we D
 
 % Freq offset estimator states constants
 
-global pilot_baseband;
-pilot_baseband = zeros(1, Npilotbaseband);              % pilot baseband samples
-global pilot_lpf
-pilot_lpf = zeros(1, Npilotlpf);                        % LPF pilot samples
+global pilot_baseband1;
+global pilot_baseband2;
+pilot_baseband1 = zeros(1, Npilotbaseband);             % pilot baseband samples
+pilot_baseband2 = zeros(1, Npilotbaseband);             % pilot baseband samples
+global pilot_lpf1
+global pilot_lpf2
+pilot_lpf1 = zeros(1, Npilotlpf);                       % LPF pilot samples
+pilot_lpf2 = zeros(1, Npilotlpf);                       % LPF pilot samples
 
 % Timing estimator states
 

diff --git a/codec2-dev/octave/fdmdv_demod.m b/codec2-dev/octave/fdmdv_demod.m
index 29cd6dc6892b684f6c8437292c990c0fe510b34a..ed26e7a256c345d1daf0bf5121c3c411728d5ea1 100644 (file)
--- a/codec2-dev/octave/fdmdv_demod.m
+++ b/codec2-dev/octave/fdmdv_demod.m
@@ -33,6 +33,7 @@ function fdmdv_demod(rawfilename, nbits)
   pilot_freq = freq(Nc+1);
   pilot_phase = 1;
   pilot_filter_mem = zeros(1, Nfilter);
+  prev_pilot = zeros(M,1);
 
   % BER stats
 
@@ -40,6 +41,8 @@ function fdmdv_demod(rawfilename, nbits)
   total_bits = 0;
   bit_errors_log = [];
   sync_log = [];
+  test_frame_sync_log = [];
+  test_frame_sync_state = 0;
 
   rx_symbols_log = [];
   rx_timing_log = [];
@@ -53,8 +56,10 @@ function fdmdv_demod(rawfilename, nbits)
     % frequency offset estimation and correction
 
     [pilot pilot_rx_bit pilot_symbol pilot_filter_mem pilot_phase] = generate_pilot_fdm(pilot_rx_bit, pilot_symbol, pilot_filter_mem, pilot_phase, pilot_freq);
-    foff = rx_est_freq_offset(rx_fdm, pilot);
+    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
@@ -80,16 +85,38 @@ function fdmdv_demod(rawfilename, nbits)
     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
 
     [test_frame_sync bit_errors] = put_test_bits(rx_bits);
     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];
-    else
-      bit_errors_log = [bit_errors_log -1];
     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
 
   % ---------------------------------------------------------------------
@@ -134,11 +161,16 @@ function fdmdv_demod(rawfilename, nbits)
 
   figure(4)
   clf;
-  subplot(211)
+  subplot(311)
   stem(sync_log)
+  axis([0 frames 0 1.5]);
   title('BPSK Sync')
-  subplot(212)
+  subplot(312)
   stem(bit_errors_log);
-  title('Bit Errors for test data')
+  title('Bit Errors for test frames')
+  subplot(313)
+  plot(test_frame_sync_log);
+  axis([0 frames 0 1.5]);
+  title('Test Frame Sync')
 
 endfunction

diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m
index b9de56bd4e31cd157a9b8f54cf69d2ddf80e1600..4055739762d5e86e96d08586c3e37da9b2113f66 100644 (file)
--- 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 = 50;
-EbNo_dB = 73;
+EbNo_dB = 7.3;
 Foff_hz = 0;
 modulation = 'dqpsk';
 hpa_clip = 150;
@@ -25,8 +25,6 @@ rx_phase_log = 0;
 rx_timing_log = 0;
 tx_pwr = 0;
 noise_pwr = 0;
-total_bit_errors = 0;
-total_bits = 0;
 rx_fdm_log = [];
 rx_baseband_log = [];
 rx_bits_offset = zeros(Nc*Nb*2);
@@ -35,8 +33,15 @@ prev_rx_symbols = ones(Nc+1,1);
 foff_log = [];
 tx_baseband_log = [];
 tx_fdm_log = [];
-sync_log = [];
+
+% BER stats
+
+total_bit_errors = 0;
+total_bits = 0;
 bit_errors_log = [];
+sync_log = [];
+test_frame_sync_log = [];
+test_frame_sync_state = 0;
 
 % pilot states, used for copy of pilot at rx
 
@@ -45,6 +50,7 @@ pilot_symbol = sqrt(2);
 pilot_freq = freq(Nc+1);
 pilot_phase = 1;
 pilot_filter_mem = zeros(1, Nfilter);
+prev_pilot = zeros(M,1);
 
 % fixed delay simuation
 
@@ -148,7 +154,8 @@ for f=1:frames
   % frequency offset estimation and correction
 
   [pilot pilot_rx_bit pilot_symbol pilot_filter_mem pilot_phase] = generate_pilot_fdm(pilot_rx_bit, pilot_symbol, pilot_filter_mem, pilot_phase, pilot_freq);
-  foff = rx_est_freq_offset(rx_fdm_delay, pilot);
+  foff = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot);
+  prev_pilot = pilot;
   foff_log = [ foff_log foff ];
   %foff = 0;
   foff_rect = exp(j*2*pi*foff/Fs);
@@ -188,10 +195,30 @@ for f=1:frames
     total_bit_errors = total_bit_errors + bit_errors;
     total_bits = total_bits + Ntest_bits;
     bit_errors_log = [bit_errors_log bit_errors];
-  else
-    bit_errors_log = [bit_errors_log -1];
+  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
 
 % ---------------------------------------------------------------------
@@ -216,7 +243,7 @@ printf("Eb/No (meas): %2.2f (%2.2f) dB  %d bits  %d errors  BER: (%1.4f) PAPR: %
 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');
 
@@ -243,12 +270,17 @@ title('FDM Rx Spectrum');
 
 figure(4)
 clf;
-subplot(211)
+subplot(311)
 stem(sync_log)
+axis([0 frames 0 1.5]);
 title('BPSK Sync')
-subplot(212)
+subplot(312)
 stem(bit_errors_log);
-title('Bit Errors for test data')
+title('Bit Errors for test frames')
+subplot(313)
+plot(test_frame_sync_log);
+axis([0 frames 0 1.5]);
+title('Test Frame Sync')
 
 % TODO
 %   + handling sample slips, extra plus/minus samples