From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Wed, 7 Mar 2012 23:35:45 +0000 (+0000)
Subject: first pass at pilot tone based freq offset estimation, seems to work OK
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=339252dc06da29badf4529f7216ab57454826b9c;p=freetel-svn-tracking.git

first pass at pilot tone based freq offset estimation, seems to work OK

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

diff --git a/codec2-dev/octave/README_fdmdv.txt b/codec2-dev/octave/README_fdmdv.txt
new file mode 100644
index 00000000..5d5c0e02
--- /dev/null
+++ b/codec2-dev/octave/README_fdmdv.txt
@@ -0,0 +1,5 @@
+
+
+Modeling sample clock errors using sox:
+
+sox -r 8000 -s -2 mod_dqpsk.raw -s -2 mod_dqpsk_8008hz.raw rate -h 8008
diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m
index ea06e3ff..7f5927bd 100644
--- a/codec2-dev/octave/fdmdv.m
+++ b/codec2-dev/octave/fdmdv.m
@@ -8,6 +8,11 @@
 % Version 2
 %
 
+% reqd to mak sure we get same random bits at mod and demod
+
+rand('state',1); 
+randn('state',1);
+
 % Constants
 
 global Fs = 8000;      % sample rate in Hz
@@ -80,7 +85,7 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
       msb = tx_bits_matrix(c,1); lsb = tx_bits_matrix(c,2);
 
       if ((msb == 0) && (lsb == 0))
-         tx_symbols(c) = prev_tx_symbols(c);
+	  tx_symbols(c) = prev_tx_symbols(c);
       endif  
       if ((msb == 0) && (lsb == 1))
          tx_symbols(c) = j*prev_tx_symbols(c);
@@ -94,7 +99,7 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation)
   end
   else
     % QPSK mapping
-    tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2); 
+    tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2);
   endif
 
 endfunction
@@ -145,7 +150,7 @@ function tx_fdm = fdm_upconvert(tx_filt)
 	tx_fdm(i) = tx_fdm(i) + tx_filt(c,i)*phase_tx(c);
       end
   end
-
+  
   % Nc/2 tones above centre freq  
 
   for c=Nc/2+1:Nc
@@ -155,7 +160,19 @@ function tx_fdm = fdm_upconvert(tx_filt)
       end
   end
 
-  tx_fdm = real(tx_fdm);
+  % add centre pilot tone at twice amplitude of other tones
+
+  c = Nc+1;
+  for i=1:M
+    phase_tx(c) = phase_tx(c) * freq(c);
+    tx_fdm(i) = tx_fdm(i) + 1/sqrt(2)*phase_tx(c);
+  end
+ 
+  % Scale such that total Carrier power C of real(tx_fdm) = Nc
+  % we return the complex (single sided) signal to make frequency
+  % shifting for the purpose of testing easier
+
+  tx_fdm = sqrt(2)*tx_fdm;
 endfunction
 
 
@@ -218,8 +235,62 @@ function rx_filt = rx_filter(rx_baseband)
 endfunction
 
 
-% Estimate optimum timing offset, and symbol receive symbols
+% Estimate frequency offset of FDM signal by peak picking pilot
+% tone at Fcentre
 
+function foff = rx_est_freq_offset(rx_fdm)
+  global M;
+  global Nc;
+  global Fs;
+  global freq;
+  global phase_rx;
+  global Npilotbaseband;
+  global Npilotlpf;
+  global Npilotcoeff;
+  global pilot_baseband;
+  global pilot_lpf;
+  global pilot_coeff;
+
+  % down convert latest M samples of pilot
+
+  pilot_baseband(1:Npilotbaseband-M) = pilot_baseband(M+1:Npilotbaseband);
+  c = Nc+1;
+  for i=1:M
+    phase_rx(c) = phase_rx(c) * freq(c);
+    pilot_baseband(Npilotbaseband-M+i) = rx_fdm(i) * phase_rx(c)';
+  end
+
+  % LPF cutoff 200Hz
+
+  pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf);
+  j = 1;
+  for i = Npilotlpf-M+1:Npilotlpf
+    pilot_lpf(i) = pilot_baseband(j:j+Npilotcoeff) * pilot_coeff';
+    j++;
+  end
+
+  % 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));
+
+  % peak pick and convert to Hz
+
+  [x ix] = max(S);
+  r = 2*200/Mpilotfft;     % maps FFT bin to frequency in Hz
+
+  if ix > Mpilotfft/2
+    foff = (ix - Mpilotfft - 1)*r;
+  else
+    foff = (ix - 1)*r;
+  endif
+
+endfunction
+
+
+% Estimate optimum timing offset, and symbol receive symbols
 
 function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
   global M;
@@ -250,7 +321,7 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband)
   % map phase to estimated optimum timing instant at rate M
   % the M/2 + 1 part was adjusted by experment, I know not why....
 
-  rx_timing = angle(x)*M/pi + M/2 + 1;
+  rx_timing = angle(x)*M/(2*pi) + M/4;
   if (rx_timing > M)
      rx_timing -= M;
   end
@@ -277,7 +348,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
 % all Nc carriers which should increase the SNR of phase estimate.
-% Maybe phase is coherent over a coupel of symbols in HF channel,not
+% 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.
 
@@ -369,10 +440,10 @@ function [sync bit_errors] = put_test_bits(rx_bits)
 
   % if less than a thresh we are aligned and in sync
 
-  ber =  bit_errors/Ntest_bits;
+  ber = bit_errors/Ntest_bits;
   
   sync = 0;
-  if (ber < 0.2)
+  if (ber < 0.1)
     sync = 1;
   endif
 endfunction
@@ -385,7 +456,7 @@ global rx_filter_memory = zeros(Nc, Nfilter);
 
 % phasors used for up and down converters
 
-global freq = zeros(Nc,1);;
+global freq = zeros(Nc+1,1);;
 for c=1:Nc/2
   carrier_freq = (-Nc/2 - 1 + c)*Fsep + Fcentre;
   freq(c) = exp(j*2*pi*carrier_freq/Fs);
@@ -395,10 +466,21 @@ for c=Nc/2+1:Nc
   freq(c) = exp(j*2*pi*carrier_freq/Fs);
 end
 
-global phase_tx = ones(Nc,1);
-global phase_rx = ones(Nc,1);
+freq(Nc+1) = exp(j*2*pi*Fcentre/Fs);
+
+global phase_tx = ones(Nc+1,1);
+global phase_rx = ones(Nc+1,1);
+
+% Freq offset estimator states
+
+global Npilotcoeff    = 30;                             % number of pilot LPF coeffs
+global pilot_coeff    = fir1(Npilotcoeff, 200/(Fs/2))'; % 200Hz LPF
+global Npilotbaseband = Npilotcoeff + 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
 
-% Timing estimateor states
+% Timing estimator states
 
 global rx_filter_mem_timing = zeros(Nc, Nt*P);
 global rx_baseband_mem_timing = zeros(Nc, Nfiltertiming);
diff --git a/codec2-dev/octave/fdmdv_demod.m b/codec2-dev/octave/fdmdv_demod.m
index 9cc707e6..7aff775b 100644
--- a/codec2-dev/octave/fdmdv_demod.m
+++ b/codec2-dev/octave/fdmdv_demod.m
@@ -11,6 +11,11 @@ function fdmdv_demod(rawfilename, nbits)
 
 fdmdv; % include modem code
 
+modulation = 'dqpsk';
+Foff_hz = 0;
+phase_offset = 1;
+freq_offset = exp(j*2*pi*Foff_hz/Fs);
+
 fin = fopen(rawfilename, "rb");
 rx_fdm = fread(fin, Inf, "short");
 gain = 1000;
@@ -18,21 +23,26 @@ rx_fdm /= gain;
 if (nargin == 1)
   frames = floor(length(rx_fdm)/M);
 else
-  frames = nbits/(Nc*Nb);
+    frames = nbits/(Nc*Nb);
 endif
 
+prev_rx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
+
 total_bit_errors = 0;
 total_bits = 0;
 
 rx_timing_log = [];
 rx_symbols_log = [];
 rx_phase_log = [];
-prev_rx_symbols = ones(Nc,1)*exp(j*pi/4);
-modulation = 'dqpsk';
+sync_log = [];
 
 % Main loop ----------------------------------------------------
 
 for i=1:frames
+  for n=1:M
+    phase_offset *= freq_offset;
+    rx_fdm((i-1)*M+1+n) *= phase_offset;
+  end
   rx_baseband = fdm_downconvert(rx_fdm((i-1)*M+1:i*M));
   rx_filt = rx_filter(rx_baseband);
 
@@ -43,6 +53,7 @@ for i=1:frames
   %rx_phase_log = [rx_phase_log rx_phase];
   %rx_symbols = rx_symbols*exp(j*rx_phase);
 
+
   if strcmp(modulation,'dqpsk')
     rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)];
   else
@@ -52,6 +63,7 @@ for i=1:frames
   prev_rx_symbols = rx_symbols;
 
   [sync bit_errors] = put_test_bits(rx_bits);
+  sync_log = [sync_log sync];
 
   if sync == 1
     total_bit_errors = total_bit_errors + bit_errors;
@@ -69,14 +81,18 @@ figure(1)
 clf;
 [n m] = size(rx_symbols_log);
 plot(real(rx_symbols_log(:,20:m)),imag(rx_symbols_log(:,20:m)),'+')
+%plot(rx_fdm);
 figure(2)
 clf;
-subplot(311)
+%subplot(211)
 plot(rx_timing_log)
-subplot(312)
-Nfft=Fs;
-S=fft(rx_fdm,Nfft);
-SdB=20*log10(abs(S));
-plot(SdB(1:Fs/4))
-subplot(313)
-%plot(rx_phase_log)
+%subplot(212)
+%plot(sync_log)
+figure(3)
+clf;
+for i=1:Nc
+  subplot(4,4,i);
+  mx = max(abs(rx_symbols_log(i,20:m)));
+  plot(real(rx_symbols_log(i,20:m)),imag(rx_symbols_log(i,20:m)),'+')
+  axis([-mx mx -mx mx]);
+end
diff --git a/codec2-dev/octave/fdmdv_mod.m b/codec2-dev/octave/fdmdv_mod.m
index 69a310c1..589047f6 100644
--- a/codec2-dev/octave/fdmdv_mod.m
+++ b/codec2-dev/octave/fdmdv_mod.m
@@ -14,14 +14,14 @@ fdmdv; % include modem code
 frames = floor(nbits/(Nc*Nb));
 tx_fdm = [];
 gain = 1000; % Scale up to 16 bit shorts
-prev_tx_symbols = ones(Nc,1)*exp(j*pi/4);
+prev_tx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
 
 for i=1:frames
   tx_bits = get_test_bits(Nc*Nb);
-  tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits,'dqpsk');
+  tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits,'qpsk');
   prev_tx_symbols = tx_symbols;
   tx_baseband = tx_filter(tx_symbols);
-  tx_fdm = [tx_fdm fdm_upconvert(tx_baseband)];
+tx_fdm = [tx_fdm real(fdm_upconvert(tx_baseband))];
 end
 
 tx_fdm *= gain;
diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m
index 562da5c4..96652504 100644
--- a/codec2-dev/octave/fdmdv_ut.m
+++ b/codec2-dev/octave/fdmdv_ut.m
@@ -10,14 +10,12 @@
 clear all;
 
 fdmdv;               % load modem code
-rand('state',1); 
-randn('state',1);
  
 % Simulation Parameters --------------------------------------
 
-frames = 50;
-EbNo_dB = 40;
-Foff_hz = 1;
+frames = 100;
+EbNo_dB = 7.3;
+Foff_hz = 124;
 modulation = 'dqpsk';
 
 % ------------------------------------------------------------
@@ -31,64 +29,103 @@ noise_pwr = 0;
 total_bit_errors = 0;
 total_bits = 0;
 rx_fdm_log = [];
+rx_baseband_log = [];
 rx_bits_offset = zeros(Nc*Nb*2);
-prev_tx_symbols = ones(Nc,1)*exp(j*pi/4);
-prev_rx_symbols = ones(Nc,1)*exp(j*pi/4);
+prev_tx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
+prev_rx_symbols = sqrt(2)*ones(Nc,1)*exp(j*pi/4);
+foff_log = [];
+
+Ndelay = M+20;
+rx_fdm_delay = zeros(Ndelay,1);
 
 % Eb/No calculations.  We need to work out Eb/No for each FDM carrier.
 % Total power is sum of power in all FDM carriers
 
-C = 1;  % power of each FDM carrier (energy/sample)
-N = 1;  % total noise power (energy/sample) of noise source before scaling 
-        % by Ngain
+C = 1; % power of each FDM carrier (energy/sample).  Total Carrier power should = Nc*C = Nc
+N = 1; % total noise power (energy/sample) of noise source across entire bandwidth
 
 % Eb  = Carrier power * symbol time / (bits/symbol)
-%     = C *(Rs/Fs) / 2
-Eb_dB = 10*log10(C) + 10*log10(Rs) - 10*log10(Fs) - 10*log10(2);
+%     = C *(1/Rs) / 2
+Eb_dB = 10*log10(C) - 10*log10(Rs) - 10*log10(2);
 
 No_dBHz = Eb_dB - EbNo_dB;
 
-% Noise power = Noise spectral density * bandwidth;
-N_dB = No_dBHz + 10*log10(Fs);
+% Noise power = Noise spectral density * bandwidth
+% Noise power = Noise spectral density * Fs/2 for real signals
+N_dB = No_dBHz + 10*log10(Fs/2);
 Ngain_dB = N_dB - 10*log10(N);
 Ngain = 10^(Ngain_dB/20);
 
-% C/No = Carrier Power/noise spectral denity
-%      = power per carrier*number of carriers / noise spectral denity
+% C/No = Carrier Power/noise spectral density
+%      = power per carrier*number of carriers / noise spectral density
 CNo_dB = 10*log10(C)  + 10*log10(Nc) - No_dBHz;
 
-% SNR in equivalent 2400 Hz SSB channel
+% SNR in equivalent 3000 Hz SSB channel
 
-B = 2400;
+B = 3000;
 SNR = CNo_dB - 10*log10(B);
 
 phase_offset = 1;
 freq_offset = exp(j*2*pi*Foff_hz/Fs);
+foff_phase = 1;
 
 % Main loop ----------------------------------------------------
 
 for i=1:frames
+
+  % -------------------
+  % Modulator
+  % -------------------
+
   tx_bits = get_test_bits(Nc*Nb);
-  %tx_bits = zeros(1,Nc*Nb);
-  %tx_bits = ones(1,Nc*Nb);
-  %tx_bits = [0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
-  %tx_bits = [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0];
   tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation);
   prev_tx_symbols = tx_symbols;
   tx_baseband = tx_filter(tx_symbols);
-  tx_fdm = fdm_upconvert(tx_baseband);
-  tx_pwr = 0.9*tx_pwr + 0.1*tx_fdm*tx_fdm'/(M);
+  tx_fdm = fdm_upconvert(tx_baseband) + 2*cos(2*pi*(0:M-1)*Fcentre/Fs);
+  tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M);
+
+  % -------------------
+  % Channel simulation
+  % -------------------
+
+  % frequency offset
 
-  noise = Ngain/sqrt(2)*[randn(1,M) + j*randn(1,M)];
-  noise_pwr = 0.9*noise_pwr + 0.1*noise*noise'/M;
   for i=1:M
     phase_offset *= freq_offset;
-    rx_fdm(i) = phase_offset*tx_fdm(i);
+    rx_fdm(i) = phase_offset*real(tx_fdm(i));
   end
+
+  % AWGN noise
+
+  noise = Ngain*randn(1,M);
+  noise_pwr = 0.9*noise_pwr + 0.1*noise*noise'/M;
   rx_fdm += noise;
   rx_fdm_log = [rx_fdm_log rx_fdm];
 
-  rx_baseband = fdm_downconvert(rx_fdm);
+  % delay
+
+  rx_fdm_delay(1:Ndelay-M) = rx_fdm_delay(M+1:Ndelay);
+  rx_fdm_delay(Ndelay-M+1:Ndelay) = rx_fdm;
+
+  % -------------------
+  % Demodulator
+  % -------------------
+
+  % frequency offset estimation and correction
+
+  foff = rx_est_freq_offset(rx_fdm);
+  foff_log = [ foff_log foff ];
+  foff_rect = exp(j*2*pi*foff/Fs);
+
+  for i=1:M
+    foff_phase *= foff_rect';
+    rx_fdm_delay(i) = rx_fdm_delay(i)*foff_phase;
+  end
+
+  % baseband processing
+
+  rx_baseband = fdm_downconvert(rx_fdm_delay(1:M));
+  rx_baseband_log = [rx_baseband_log rx_baseband];
   rx_filt = rx_filter(rx_baseband);
 
   [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband);
@@ -106,6 +143,8 @@ for i=1:frames
   rx_bits = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
   prev_rx_symbols = rx_symbols;
 
+  % count bit errors
+
   [sync bit_errors] = put_test_bits(rx_bits);
   if (sync == 1)
     total_bit_errors = total_bit_errors + bit_errors;
@@ -115,23 +154,28 @@ for i=1:frames
 end
 
 ber = total_bit_errors/total_bits;
-printf("Eb/No: %2.2f dB  %d bits  %d errors  Meas BER: %1.4f  Theor BER: %1.4f\n", EbNo_dB, 
-      total_bits, total_bit_errors, ber, 0.5*erfc(sqrt(10.^(EbNo_dB/10))));
+printf("Eb/No (meas): %2.2f (%2.2f) dB  %d bits  %d errors  QPSK BER (meas): %1.4f (%1.4f)\n", 
+       EbNo_dB, 10*log10(0.25*tx_pwr*Fs/(Rs*Nc*noise_pwr)),
+       total_bits, total_bit_errors, 0.5*erfc(sqrt(10.^(EbNo_dB/10))), ber );
 
 figure(1)
 clf;
 [n m] = size(rx_symbols_log);
 plot(real(rx_symbols_log(:,20:m)),imag(rx_symbols_log(:,20:m)),'+')
+
 figure(2)
 clf;
 subplot(211)
 plot(rx_timing_log)
 subplot(212)
-Nfft=Fs;
-S=fft(rx_fdm_log,Nfft);
-SdB=20*log10(abs(S));
+%Nfft=Fs;
+%S=fft(rx_fdm_log,Nfft);
+%SdB=20*log10(abs(S));
 %plot(-Fs/2+1:Fs/2,fftshift(SdB))
-plot(SdB(1:Fs/4))
+%plot(SdB(1:Fs/4))
+plot(foff_log)
+
+
 
 % TODO
 %   + handling sample slips, extra plus/minus samples
@@ -163,3 +207,11 @@ plot(SdB(1:Fs/4))
 %   Sensitivity to Fs
 %   Try BPSK
 %   second term of QPSK eqn
+
+% Faster sync:
+%
+% 1/ Maybe Ask Bill how we can sync in less than 10 symbols?  What to
+% put in filter memories?  If high SNR should sync very quickly
+% Maybe start feeding in symbols to middle of filter memory?  Or do timing
+% sync before rx filtering at start?
+