tracked down phase difference bug

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

diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m
index 4640a7fa035fbd7e8f9134f78ca9b362d70fba74..7fdabd4e5299a003fc67bda66813e1a5291f32a3 100644 (file)
--- a/codec2-dev/octave/fdmdv.m
+++ b/codec2-dev/octave/fdmdv.m
@@ -506,8 +506,11 @@ function [rx_bits sync_bit f_err phase_difference] = psk_to_bits(prev_rx_symbols
   assert((m == 4) || (m == 8));
 
   phase_difference = zeros(Nc+1,1);
-  phase_difference(1:Nc) = rx_symbols(1:Nc) .* conj(prev_rx_symbols(1:Nc)./(1E-6+abs(prev_rx_symbols(1:Nc))));
-  
+  for c=1:Nc 
+    norm = 1/(1E-6+abs(prev_rx_symbols(c)));  
+    phase_difference(c) = prev_rx_symbols(c) .* conj(prev_rx_symbols(c)) * norm;
+  end
+
   for c=1:Nc
 
     % determine index of constellation point received 0,1,...,m-1
@@ -542,7 +545,8 @@ function [rx_bits sync_bit f_err phase_difference] = psk_to_bits(prev_rx_symbols
 
   % Extract DBPSK encoded Sync bit
 
-  phase_difference(Nc+1,1) = rx_symbols(Nc+1) .* conj(prev_rx_symbols(Nc+1)./(1E-6+abs(prev_rx_symbols(Nc+1))));
+  norm = 1/(1E-6+abs(prev_rx_symbols(Nc+1)));
+  phase_difference(Nc+1) = rx_symbols(Nc+1) * conj(prev_rx_symbols(Nc+1)) * norm;
   if (real(phase_difference(Nc+1)) < 0)
     sync_bit = 1;
     f_err = imag(phase_difference(Nc+1));
@@ -552,9 +556,9 @@ function [rx_bits sync_bit f_err phase_difference] = psk_to_bits(prev_rx_symbols
   end
 
   % extra pi/4 rotation as we need for snr_update and scatter diagram
-
+  
   phase_difference *= exp(j*pi/4);
-
+  
 endfunction
 
 

diff --git a/codec2-dev/octave/tfdmdv.m b/codec2-dev/octave/tfdmdv.m
index 346ba0ec26760c108ccda93b3a5a8bc0cd29e4fb..86907f28009cf164b13ade0280d77a6524d03eee 100644 (file)
--- a/codec2-dev/octave/tfdmdv.m
+++ b/codec2-dev/octave/tfdmdv.m
@@ -48,6 +48,7 @@ rx_baseband_log = [];
 rx_filt_log = [];
 env_log = [];
 rx_timing_log = [];
+phase_difference_log = [];
 rx_symbols_log = [];
 rx_bits_log = []; 
 sync_bit_log = [];  
@@ -56,6 +57,11 @@ nin_log = [];
 sig_est_log = [];
 noise_est_log = [];
 
+% adjust this if the screen is getting a bit cluttered
+
+global no_plot_list;
+no_plot_list = [1 2 3 4 5 6 7 8 12];
+
 for f=1:frames
 
   % modulator
@@ -137,6 +143,7 @@ for f=1:frames
   nin_log = [nin_log nin];
 
   [rx_bits sync_bit foff_fine pd] = psk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
+  phase_difference_log = [phase_difference_log pd];
 
   [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
   sig_est_log = [sig_est_log sig_est];
@@ -162,6 +169,11 @@ load ../unittest/tfdmdv_out.txt
 % Helper functions to plot output of C verson and difference between Octave and C versions
 
 function stem_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec)
+  global no_plot_list;
+
+  if find(no_plot_list == plotnum)
+    return;
+  end
   figure(plotnum)
   subplot(subplotnum)
   stem(sig,'g;C version;');
@@ -175,6 +187,12 @@ function stem_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec)
 endfunction
 
 function plot_sig_and_error(plotnum, subplotnum, sig, error, titlestr, axisvec)
+  global no_plot_list;
+
+  if find(no_plot_list == plotnum)
+    return;
+  end
+
   figure(plotnum)
   subplot(subplotnum)
   plot(sig,'g;C version;');
@@ -243,22 +261,26 @@ plot_sig_and_error(11, 212, imag(rx_baseband_log(c,:)), imag(rx_baseband_log(c,:
 plot_sig_and_error(12, 211, real(rx_filt_log(c,:)), real(rx_filt_log(c,:) - rx_filt_log_c(c,:)), 'Rx filt real' )
 plot_sig_and_error(12, 212, imag(rx_filt_log(c,:)), imag(rx_filt_log(c,:) - rx_filt_log_c(c,:)), 'Rx filt imag' )
 
-st=1; en=3*Nt*P;
-plot_sig_and_error(13, 211, env_log(st:en), env_log(st:en) - env_log_c(st:en), 'env' )
-stem_sig_and_error(13, 212, real(rx_symbols_log(c,:)), real(rx_symbols_log(c,:) - rx_symbols_log_c(c,:)), 'rx symbols' )
-
-st=10*28;
-en = 12*28;
+st=1*28;
+en = 3*28;
 plot_sig_and_error(14, 211, rx_timing_log, rx_timing_log - rx_timing_log_c, 'Rx Timing' )
 stem_sig_and_error(14, 212, sync_bit_log_c, sync_bit_log - sync_bit_log_c, 'Sync bit', [1 n -1.5 1.5])
 
 stem_sig_and_error(15, 211, rx_bits_log_c(st:en), rx_bits_log(st:en) - rx_bits_log_c(st:en), 'RX bits', [1 en-st -1.5 1.5])
 stem_sig_and_error(15, 212, nin_log_c, nin_log - nin_log_c, 'nin')
 
-c = 1;
+c = 12;
 plot_sig_and_error(16, 211, sig_est_log(c,:), sig_est_log(c,:) - sig_est_log_c(c,:), 'sig est for SNR' )
 plot_sig_and_error(16, 212, noise_est_log(c,:), noise_est_log(c,:) - noise_est_log_c(c,:), 'noise est for SNR' )
 
+f=2;
+
+stem_sig_and_error(13, 211, real(rx_symbols_log(:,f)), real(rx_symbols_log(:,f) - rx_symbols_log_c(:,f)), 'rx symbols real' )
+stem_sig_and_error(13, 212, imag(rx_symbols_log(:,f)), imag(rx_symbols_log(:,f) - rx_symbols_log_c(:,f)), 'rx symbols imag' )
+
+stem_sig_and_error(17, 211, real(phase_difference_log(:,f)), real(phase_difference_log(:,f) - phase_difference_log_c(:,f)), 'phase difference real' )
+stem_sig_and_error(17, 212, imag(phase_difference_log(:,f)), imag(phase_difference_log(:,f) - phase_difference_log_c(:,f)), 'phase difference imag' )
+
 % ---------------------------------------------------------------------------------------
 % AUTOMATED CHECKS ------------------------------------------
 % ---------------------------------------------------------------------------------------

diff --git a/codec2-dev/src/fdmdv.c b/codec2-dev/src/fdmdv.c
index f912c5d3eab3a171f3a3ae6f922410e75b78c1c1..757da9b5440ffc1cce5603d62f73fb3debe27315 100644 (file)
--- a/codec2-dev/src/fdmdv.c
+++ b/codec2-dev/src/fdmdv.c
@@ -106,7 +106,7 @@ static COMP cadd(COMP a, COMP b)
 
 static float cabsolute(COMP a)
 {
-    return sqrt(pow(a.real, 2.0) + pow(a.imag, 2.0));
+    return sqrt(a.real*a.real + a.imag*a.imag);
 }
 
 /*---------------------------------------------------------------------------*\
@@ -1019,8 +1019,10 @@ float qpsk_to_bits(int rx_bits[], int *sync_bit, int Nc, COMP phase_difference[]
        leads to sensible scatter plots */
 
     for(c=0; c<Nc; c++) {
-        norm = 1.0/(cabsolute(prev_rx_symbols[c])+1E-6);
-       phase_difference[c] = cmult(cmult(rx_symbols[c], fcmult(norm,cconj(prev_rx_symbols[c]))), pi_on_4);
+        norm = 1.0/(1E-6 + cabsolute(prev_rx_symbols[c]));
+       prev_rx_symbols[c] = fcmult(norm, prev_rx_symbols[c]);
+        phase_difference[c] = cmult(prev_rx_symbols[c], cconj(prev_rx_symbols[c]));
+       phase_difference[c] = cmult(phase_difference[c],pi_on_4);        
     }
                                    
     /* map (Nc,1) DQPSK symbols back into an (1,Nc*Nb) array of bits */

diff --git a/codec2-dev/unittest/tfdmdv.c b/codec2-dev/unittest/tfdmdv.c
index 07eb896ff1ba2a4a3a660558538d0bedb67c7773..256ae55675429b99c6a53f1f1f102aaf7b442460 100644 (file)
--- a/codec2-dev/unittest/tfdmdv.c
+++ b/codec2-dev/unittest/tfdmdv.c
@@ -85,6 +85,7 @@ int main(int argc, char *argv[])
     float         env_log[NT*P*FRAMES];
     float         rx_timing_log[FRAMES];
     COMP          rx_symbols_log[FDMDV_NC+1][FRAMES];
+    COMP          phase_difference_log[FDMDV_NC+1][FRAMES];
     float         sig_est_log[FDMDV_NC+1][FRAMES];
     float         noise_est_log[FDMDV_NC+1][FRAMES];
     int           rx_bits_log[FDMDV_BITS_PER_FRAME*FRAMES];
@@ -106,7 +107,7 @@ int main(int argc, char *argv[])
     printf("sizeof FDMDV states: %d bytes\n", sizeof(struct FDMDV));
 
     for(f=0; f<FRAMES; f++) {
-
+        
        /* --------------------------------------------------------*\
                                  Modulator
        \*---------------------------------------------------------*/
@@ -170,6 +171,11 @@ int main(int argc, char *argv[])
        rx_filter(rx_filt, FDMDV_NC, rx_baseband, fdmdv->rx_filter_memory, nin);
        rx_timing = rx_est_timing(rx_symbols, FDMDV_NC, rx_filt, rx_baseband, fdmdv->rx_filter_mem_timing, env, fdmdv->rx_baseband_mem_timing, nin);     
        foff_fine = qpsk_to_bits(rx_bits, &sync_bit, FDMDV_NC, fdmdv->phase_difference, fdmdv->prev_rx_symbols, rx_symbols, 0);
+        //for(i=0; i<FDMDV_NC;i++)
+        //    printf("rx_symbols: %f %f prev_rx_symbols: %f %f phase_difference: %f %f\n", rx_symbols[i].real, rx_symbols[i].imag,
+        //          fdmdv->prev_rx_symbols[i].real, fdmdv->prev_rx_symbols[i].imag, fdmdv->phase_difference[i].real, fdmdv->phase_difference[i].imag);
+        //if (f==1)
+        //   exit(0);
        snr_update(fdmdv->sig_est, fdmdv->noise_est, FDMDV_NC, fdmdv->phase_difference);
        memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP)*(FDMDV_NC+1));
        
@@ -195,8 +201,6 @@ int main(int argc, char *argv[])
                tx_baseband_log[c][f*M+i] = tx_baseband[c][i]; 
        memcpy(&tx_fdm_log[M*f], tx_fdm, sizeof(COMP)*M);
 
-       /* freq offset estimation */
-
        memcpy(&pilot_baseband1_log[f*NPILOTBASEBAND], fdmdv->pilot_baseband1, sizeof(COMP)*NPILOTBASEBAND);
        memcpy(&pilot_baseband2_log[f*NPILOTBASEBAND], fdmdv->pilot_baseband2, sizeof(COMP)*NPILOTBASEBAND);
        memcpy(&pilot_lpf1_log[f*NPILOTLPF], fdmdv->pilot_lpf1, sizeof(COMP)*NPILOTLPF);
@@ -227,8 +231,10 @@ int main(int argc, char *argv[])
        memcpy(&env_log[NT*P*f], env, sizeof(float)*NT*P);
        rx_timing_log[f] = rx_timing;
        nin_log[f] = nin;
-       for(c=0; c<FDMDV_NC+1; c++)
+       for(c=0; c<FDMDV_NC+1; c++) {
            rx_symbols_log[c][f] = rx_symbols[c];
+           phase_difference_log[c][f] = fdmdv->phase_difference[c];
+        }
        
        /* qpsk_to_bits() */
 
@@ -271,6 +277,7 @@ int main(int argc, char *argv[])
     octave_save_float(fout, "env_log_c", env_log, 1, NT*P*FRAMES, NT*P*FRAMES);  
     octave_save_float(fout, "rx_timing_log_c", rx_timing_log, 1, FRAMES, FRAMES);  
     octave_save_complex(fout, "rx_symbols_log_c", (COMP*)rx_symbols_log, (FDMDV_NC+1), FRAMES, FRAMES);  
+    octave_save_complex(fout, "phase_difference_log_c", (COMP*)phase_difference_log, (FDMDV_NC+1), FRAMES, FRAMES);  
     octave_save_float(fout, "sig_est_log_c", (float*)sig_est_log, (FDMDV_NC+1), FRAMES, FRAMES);  
     octave_save_float(fout, "noise_est_log_c", (float*)noise_est_log, (FDMDV_NC+1), FRAMES, FRAMES);  
     octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1, FDMDV_BITS_PER_FRAME*FRAMES);