another pass at freq offset estimation, this time using DBPSK pilot

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

diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m
index 82d3bd374c2098a9bffe707f63611221eaaaeaee..d9a2a442060970a8fc32263b3cc4212493c2b170 100644 (file)
--- a/codec2-dev/octave/fdmdv.m
+++ b/codec2-dev/octave/fdmdv.m
@@ -28,7 +28,7 @@ global Fsep  = 75;     % Separation between carriers (Hz)
 global Fcentre = 1200; % Centre frequency, Nc/2 below this, N/c above (Hz)
 global Nt = 5;         % number of symbols we estimate timing over
 global P = 4;          % oversample factor used for rx symbol filtering
-global pilot_phase = pi/4;% current phase of symbol used to make pilot carrier
+global pilot_bit = 0;  % current phase of symbol used to make pilot carrier
 
 % Generate root raised cosine (Root Nyquist) filter ---------------
 
@@ -73,7 +73,7 @@ global gt_alpha5_root = real((ifft_GF_alpha5_root(1:Nfilter)));
 function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
   global Nc;
   global Nb;
-  global pilot_phase;
+  global pilot_bit;
 
   % re-arrange as (Nc,Nb) matrix
 
@@ -108,8 +108,16 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
   % +1 -1 +1 -1 BPSK sync carrier, once filtered becomes two spectral
   % lines at +/- Rs/2
  
-  pilot_phase *= -1;
-  tx_symbols(c+1) = pilot_phase;
+  if pilot_bit
+     tx_symbols(Nc+1) = -prev_tx_symbols(c+1);
+  else
+     tx_symbols(Nc+1) = prev_tx_symbols(c+1);
+  end
+  if pilot_bit 
+    pilot_bit = 0;
+  else
+    pilot_bit = 1;
+  end
 
 endfunction
 
@@ -141,7 +149,7 @@ endfunction
 % Construct FDM signal by frequency shifting each filtered symbol
 % stream
 
-function tx_fdm = fdm_upconvert(tx_filt)
+function [tx_fdm pilot] = fdm_upconvert(tx_filt)
   global Fs;
   global M;
   global Nc;
@@ -174,7 +182,8 @@ function tx_fdm = fdm_upconvert(tx_filt)
   c = Nc+1;
   for i=1:M
     phase_tx(c) = phase_tx(c) * freq(c);
-    tx_fdm(i) = tx_fdm(i) + 4*tx_filt(c,i)*phase_tx(c);
+    pilot(i) = sqrt(2)*2*tx_filt(c,i)*phase_tx(c);
+    tx_fdm(i) = tx_fdm(i) + pilot(i);
   end
  
   % Scale such that total Carrier power C of real(tx_fdm) = Nc.  This
@@ -215,6 +224,15 @@ function rx_baseband = fdm_downconvert(rx_fdm)
        rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
       end
   end
+
+  % Pilot
+
+  c = Nc+1;
+  for i=1:M
+    phase_rx(c) = phase_rx(c) * freq(c);
+    rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)';
+  end
+
 endfunction
 
 
@@ -229,7 +247,7 @@ function rx_filt = rx_filter(rx_baseband)
   global gt_alpha5_root;
   global Fsep;
 
-  rx_filt = zeros(Nc,P);
+  rx_filt = zeros(Nc+1,P);
 
   % 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
@@ -245,10 +263,10 @@ function rx_filt = rx_filter(rx_baseband)
 endfunction
 
 
-% Estimate frequency offset of FDM signal by peak picking pilot
-% tone at Fcentre
+% 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)
+function foff = rx_est_freq_offset(rx_fdm, pilot)
   global M;
   global Nc;
   global Fs;
@@ -265,20 +283,15 @@ function foff = rx_est_freq_offset(rx_fdm)
   global pilot_coeff;
 
   % down convert latest M samples of pilot by multiplying by
-  % ideal two-tone BPSK signal
+  % ideal two-tone BPSK pilot signal
  
   pilot_baseband(1:Npilotbaseband-M) = pilot_baseband(M+1:Npilotbaseband);
   c = Nc+1;
   for i=1:M
-    phase_rx_pilot(1) *= freq_rx_pilot(1);
-    phase_rx_pilot(2) *= freq_rx_pilot(2);
-    % phase_rx_pilot(1) *= freq(c);
-    pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i) * (phase_rx_pilot(1) + phase_rx_pilot(2)); 
-    %pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i) * phase_rx_pilot(1); 
-    %pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i); 
+    pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i)*conj(pilot(i)); 
   end
 
-  % LPF cutoff 200Hz
+  % LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset
 
   pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf);
   j = 1;
@@ -295,12 +308,11 @@ function foff = rx_est_freq_offset(rx_fdm)
   S = abs(fft(s, Mpilotfft));
 
   figure(3)
-  clf
-  subplot(211)
-  plot(real(pilot_baseband));
-  subplot(212)
-  plot(abs(fft(pilot_baseband)));
-  %plot(S);
+  %plot(real(pilot_baseband))
+  plot(real(s))
+  figure(4)
+  %plot(abs(fft(pilot_baseband)))
+  plot(S)
 
   % peak pick and convert to Hz
 
@@ -321,6 +333,7 @@ endfunction
 function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
   global M;
   global Nt;
+  global Nc;
   global rx_filter_mem_timing;
   global rx_baseband_mem_timing;
   global P;
@@ -331,7 +344,7 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
   % 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;
+  rx_filter_mem_timing(:,(Nt-1)*P+1:Nt*P) = rx_filt(1:Nc,:);
 
   % sum envelopes of all carriers
 
@@ -372,7 +385,7 @@ endfunction
 
 
 % Phase estimation function - probably won't work over a HF channel
-% Tries to operate over a sinle symbol but uses phase information from
+% 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
@@ -419,6 +432,7 @@ function rx_bits = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation)
       rx_bits(2*(c-1)+1) = msb;
       rx_bits(2*(c-1)+2) = lsb;
     end
+
   else
     % map (Nc,1) QPSK symbols back into an (1,Nc*Nb) array of bits
 
@@ -478,7 +492,7 @@ endfunction
 % Initialise ----------------------------------------------------
 
 global tx_filter_memory = zeros(Nc+1, Nfilter);
-global rx_filter_memory = zeros(Nc, Nfilter);
+global rx_filter_memory = zeros(Nc+1, Nfilter);
 
 % phasors used for up and down converters
 
@@ -501,17 +515,17 @@ global phase_rx = ones(Nc+1,1);
 
 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
+global Npilotbaseband = Npilotcoeff + 4*M;              % number of pilot baseband samples reqd for pilot LPF
 global Npilotlpf      = 4*M;                            % number of samples we DFT pilot over, pilot est window
 global pilot_baseband = zeros(1, Npilotbaseband);       % pilot baseband samples
 global pilot_lpf      = zeros(1, Npilotlpf);            % LPC pilot samples
 global phase_rx_pilot = [1 1];
-global freq_rx_pilot  = [ exp(-j*2*pi*(Fcentre-Rs/2)/Fs) exp(-j*2*pi*(Fcentre+Rs/2)/Fs) ];
+global freq_rx_pilot  = [ exp(-j*2*pi*(Fcentre-Rs/4)/Fs) exp(-j*2*pi*(Fcentre+Rs/4)/Fs) ];
 
 % Timing estimator states
 
 global rx_filter_mem_timing = zeros(Nc, Nt*P);
-global rx_baseband_mem_timing = zeros(Nc, Nfiltertiming);
+global rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming);
 
 % Test bit stream state variables
 

diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m
index 4e8588d32a11b85638aab1a0cdde9d3c14e96f34..7f0f421551c3d0049ffd14282b7946f946ba794b 100644 (file)
--- a/codec2-dev/octave/fdmdv_ut.m
+++ b/codec2-dev/octave/fdmdv_ut.m
@@ -15,13 +15,13 @@ fdmdv;               % load modem code
 
 frames = 50;
 EbNo_dB = 7.3;
-Foff_hz = 0;
+Foff_hz = -100;
 modulation = 'dqpsk';
 
 % ------------------------------------------------------------
 
 tx_filt = zeros(Nc,M);
-rx_symbols_log = zeros(Nc,1);
+rx_symbols_log = [];
 rx_phase_log = 0;
 rx_timing_log = 0;
 tx_pwr = 0;
@@ -31,8 +31,9 @@ total_bits = 0;
 rx_fdm_log = [];
 rx_baseband_log = [];
 rx_bits_offset = zeros(Nc*Nb*2);
-prev_tx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
-prev_rx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
+prev_tx_symbols = sqrt(2)*ones(Nc+1,1)*exp(j*pi/4);
+prev_tx_symbols(Nc+1) = 1;
+prev_rx_symbols = sqrt(2)*ones(Nc+1,1)*exp(j*pi/4);
 foff_log = [];
 tx_baseband_log = [];
 
@@ -84,7 +85,7 @@ for i=1:frames
   prev_tx_symbols = tx_symbols;
   tx_baseband = tx_filter(tx_symbols);
   tx_baseband_log = [tx_baseband_log tx_baseband];
-  tx_fdm = fdm_upconvert(tx_baseband);
+  [tx_fdm pilot] = fdm_upconvert(tx_baseband);
   tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M);
 
   % -------------------
@@ -94,8 +95,10 @@ for i=1:frames
   % frequency offset
 
   for i=1:M
-    Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs));
-    t++;
+    % Time varying freq offset
+    % Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs));
+    % t++;
+    Foff = Foff_hz;
     freq_offset = exp(j*2*pi*Foff/Fs);
     phase_offset *= freq_offset;
     rx_fdm(i) = phase_offset*real(tx_fdm(i));
@@ -119,7 +122,7 @@ for i=1:frames
 
   % frequency offset estimation and correction
 
-  foff = rx_est_freq_offset(rx_fdm);
+  foff = rx_est_freq_offset(rx_fdm, pilot);
   foff_log = [ foff_log foff ];
   %foff = 0;
   foff_rect = exp(j*2*pi*foff/Fs);
@@ -168,7 +171,7 @@ printf("Eb/No (meas): %2.2f (%2.2f) dB  %d bits  %d errors  QPSK BER (meas): %1.
 figure(1)
 clf;
 [n m] = size(rx_symbols_log);
-plot(real(rx_symbols_log(:,20:m)),imag(rx_symbols_log(:,20:m)),'+')
+plot(real(rx_symbols_log(1:Nc,20:m)),imag(rx_symbols_log(1:Nc,20:m)),'+')
 
 figure(2)
 clf;