Octave and C agree with freq offset and freq drift and AWGN channel, yayyyyyyyy

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

diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m
index ba1c9503ca9e72876254ed7bef91684722953bef..3106997208bd3fdbefc02351cd7ea6781d3feb48 100644 (file)
--- a/codec2-dev/octave/cohpsk.m
+++ b/codec2-dev/octave/cohpsk.m
@@ -638,12 +638,14 @@ function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_proce
           mod_strip += amod_strip;
         end
         %plot(mod_strip)
-        
+        %printf("modstrip: %f %f\n", real(mod_strip), imag(mod_strip));
+
         % loop filter made up of 1st order IIR plus integrator.  Integerator
         % was found to be reqd 
         
         afdmdv.filt = (1-beta)*afdmdv.filt + beta*angle(mod_strip);
         f_est += g*afdmdv.filt;
+        %printf("filt: %f angle: %f\n", afdmdv.filt, angle(mod_strip));
 
       end
     end

diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m
index 84c8da3a0604c57ebe4806e71b20021197d02f5c..36a1e37095969766d42c4bd338342e3ae4d1bfc6 100644 (file)
--- a/codec2-dev/octave/tcohpsk.m
+++ b/codec2-dev/octave/tcohpsk.m
@@ -10,6 +10,19 @@
 % or (ii) can optionally be used to run an Octave version of the cohpsk
 %    modem to tune and develop it.
 
+%  TODO:
+%
+%  [ ] Test
+%      [X] AWGN channel
+%      [X] freq offset
+%      [X] fading channel
+%      [X] freq drift
+%      [ ] timing drift
+%  [X] tune perf/impl loss to get closer to ideal
+%      [X] linear interp of phase for better fading perf
+%  [X] freq offset/drift feedback loop 
+%  [ ] PAPR measurement and reduction
+
 graphics_toolkit ("gnuplot");
 more off;
 
@@ -30,8 +43,8 @@ test = 'compare to c';
 
 if strcmp(test, 'compare to c')
   frames = 35;
-  foff =  0;
-  dfoff = 0;
+  foff =  55.5;
+  dfoff = -0.5/Fs;
   EsNodB = 8;
   fading_en = 0;
   hf_delay_ms = 2;
@@ -111,7 +124,7 @@ afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M;
 afdmdv.rx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter);
 
 afdmdv.filt = 0;
-afdmdv.prev_rx_symb = zeros(1,afdmdv.Nc+1);
+afdmdv.prev_rx_symb = ones(1,afdmdv.Nc+1);
 
 % COHPSK Init --------------------------------------------------------
 
@@ -156,7 +169,7 @@ ct_symb_ff_log = [];
 rx_timing_log = [];
 ratio_log = [];
 foff_log = [];
-fest_log = [];
+f_est_log = [];
 
 % Channel modeling and BER measurement ----------------------------------------
 
@@ -323,7 +336,7 @@ for f=1:frames
   % If in sync just do sample rate processing on latest frame
 
   if sync == 1
-    [ch_symb rx_timing rx_filt rx_baseband afdmdv acohpsk.f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, acohpsk.f_est, acohpsk.Nsymbrowpilot, nin, 0);
+    [ch_symb rx_timing rx_filt rx_baseband afdmdv acohpsk.f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, acohpsk.f_est, acohpsk.Nsymbrowpilot, nin, 1);
     [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb, sync, next_sync);
 
     acohpsk.ct_symb_ff_buf(1:2,:) = acohpsk.ct_symb_ff_buf(acohpsk.Nsymbrowpilot+1:acohpsk.Nsymbrowpilot+2,:);
@@ -333,7 +346,8 @@ for f=1:frames
     rx_filt_log = [rx_filt_log rx_filt];
     ch_symb_log = [ch_symb_log; ch_symb];
     rx_timing_log = [rx_timing_log rx_timing];
-    fest_log = [fest_log acohpsk.f_est];
+    f_est_log = [f_est_log acohpsk.f_est];
+    %printf("%f\n", acohpsk.f_est);
   end
 
   % if we are in sync complete demodulation with symbol rate processing
@@ -388,6 +402,10 @@ printf("\nOctave EsNodB: %4.1f ber..: %4.3f Nerrs..: %d Tbits..: %d\n", EsNodB,
 
 if compare_with_c
 
+  % Output vectors from C port ---------------------------------------------------
+
+  load ../build_linux/unittest/tcohpsk_out.txt
+
   % Determine bit error rate
 
   sz = length(tx_bits_log_c);
@@ -396,10 +414,6 @@ if compare_with_c
   ber_c = Nerrs_c/Tbits_c;
   printf("C EsNodB.....: %4.1f ber_c: %4.3f Nerrs_c: %d Tbits_c: %d\n", EsNodB, ber_c, Nerrs_c, Tbits_c);
 
-  % Output vectors from C port ---------------------------------------------------
-
-  load ../build_linux/unittest/tcohpsk_out.txt
-
   stem_sig_and_error(1, 111, tx_bits_log_c, tx_bits_log - tx_bits_log_c, 'tx bits', [1 length(tx_bits) -1.5 1.5])
   stem_sig_and_error(2, 211, real(tx_symb_log_c), real(tx_symb_log - tx_symb_log_c), 'tx symb re', [1 length(tx_symb_log_c) -1.5 1.5])
   stem_sig_and_error(2, 212, imag(tx_symb_log_c), imag(tx_symb_log - tx_symb_log_c), 'tx symb im', [1 length(tx_symb_log_c) -1.5 1.5])
@@ -415,13 +429,15 @@ if compare_with_c
   [n m] = size(ch_symb_log);
   stem_sig_and_error(6, 211, real(ch_symb_log_c), real(ch_symb_log - ch_symb_log_c), 'ch symb re', [1 n -1.5 1.5])
   stem_sig_and_error(6, 212, imag(ch_symb_log_c), imag(ch_symb_log - ch_symb_log_c), 'ch symb im', [1 n -1.5 1.5])
-  %stem_sig_and_error(7, 211, real(ct_symb_ff_log_c), real(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff re', [1 n -1.5 1.5])
-  %stem_sig_and_error(7, 212, imag(ct_symb_ff_log_c), imag(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff im', [1 n -1.5 1.5])
+
+  [n m] = size(rx_symb_log);
   stem_sig_and_error(8, 211, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'Amp Est', [1 n -1.5 1.5])
   stem_sig_and_error(8, 212, rx_phi_log_c, rx_phi_log - rx_phi_log_c, 'Phase Est', [1 n -4 4])
   stem_sig_and_error(9, 211, real(rx_symb_log_c), real(rx_symb_log - rx_symb_log_c), 'rx symb re', [1 n -1.5 1.5])
   stem_sig_and_error(9, 212, imag(rx_symb_log_c), imag(rx_symb_log - rx_symb_log_c), 'rx symb im', [1 n -1.5 1.5])
-  stem_sig_and_error(10, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 n -1.5 1.5])
+
+  stem_sig_and_error(10, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 length(rx_bits_log) -1.5 1.5])
+  stem_sig_and_error(11, 111, f_est_log_c - Fcentre, f_est_log - f_est_log_c, 'f est', [1 length(f_est_log) foff-5 foff+5])
 
   check(tx_bits_log, tx_bits_log_c, 'tx_bits');
   check(tx_symb_log, tx_symb_log_c, 'tx_symb');
@@ -429,12 +445,13 @@ if compare_with_c
   check(ch_fdm_frame_log, ch_fdm_frame_log_c, 'ch_fdm_frame');
   %check(rx_fdm_frame_bb_log, rx_fdm_frame_bb_log_c, 'rx_fdm_frame_bb', 0.01);
 
-  check(ch_symb_log, ch_symb_log_c, 'ch_symb',0.01);
+  check(ch_symb_log, ch_symb_log_c, 'ch_symb',0.05);
   %check(ct_symb_ff_log, ct_symb_ff_log_c, 'ct_symb_ff',0.01);
   check(rx_amp_log, rx_amp_log_c, 'rx_amp_log',0.01);
-  check(rx_phi_log, rx_phi_log_c, 'rx_phi_log',0.01);
+  check(rx_phi_log, rx_phi_log_c, 'rx_phi_log',0.05);
   check(rx_symb_log, rx_symb_log_c, 'rx_symb',0.01);
   check(rx_bits_log, rx_bits_log_c, 'rx_bits');
+  check(f_est_log, f_est_log_c, 'f_est');
 
 
 else
@@ -466,13 +483,13 @@ else
   subplot(211)
   plot(foff_log,';freq offset;');
   hold on;
-  plot(fest_log - Fcentre,'g;freq offset est;');
+  plot(f_est_log - Fcentre,'g;freq offset est;');
   hold off;
   title('freq offset');
   legend("boxoff");  
 
   subplot(212)
-  plot(foff_log(1:length(fest_log)) - fest_log + Fcentre)
+  plot(foff_log(1:length(f_est_log)) - fest_log + Fcentre)
   title('freq offset estimation error');
 
 end

diff --git a/codec2-dev/src/cohpsk.c b/codec2-dev/src/cohpsk.c
index 21967d97d9e77e9e2f6af3384712105fa4013dd3..f3b37e49b588ebdf7e5c0eb9e53fa5dfa7688fae 100644 (file)
--- a/codec2-dev/src/cohpsk.c
+++ b/codec2-dev/src/cohpsk.c
@@ -5,20 +5,7 @@
   DATE CREATED: March 2015
                                                                              
   Functions that implement a coherent PSK FDM modem.
-                 
-  TODO:
-
-  [ ] Code to plot EB/No v BER curves to char perf for
-      [ ] AWGN channel
-      [ ] freq offset
-      [ ] fading channel
-      [ ] freq drift
-      [ ] timing drift
-  [ ] tune perf/impl loss to get closer to ideal
-      [X] linear interp of phase for better fading perf
-  [ ] freq offset/drift feedback loop 
-  [ ] smaller freq est block size to min ram req
-                                                      
+                                                                       
 \*---------------------------------------------------------------------------*/
 
 /*
@@ -564,6 +551,8 @@ void rate_Fs_rx_processing(struct COHPSK *coh, COMP ch_symb[][COHPSK_NC*ND], COM
 
             mod_strip.real = 0.0; mod_strip.imag = 0.0;
             for(c=0; c<fdmdv->Nc+1; c++) {
+                //printf("rx_onesym[%d] %f %f prev_rx_symbols[%d] %f %f\n", c, rx_onesym[c].real, rx_onesym[c].imag,
+                //       fdmdv->prev_rx_symbols[c].real, fdmdv->prev_rx_symbols[c].imag);
                 adiff = cmult(rx_onesym[c], cconj(fdmdv->prev_rx_symbols[c]));
                 fdmdv->prev_rx_symbols[c] = rx_onesym[c];
 
@@ -573,14 +562,16 @@ void rate_Fs_rx_processing(struct COHPSK *coh, COMP ch_symb[][COHPSK_NC*ND], COM
 
                 amod_strip = cmult(adiff, adiff);
                 amod_strip = cmult(amod_strip, amod_strip);
-                amod_strip.real = cabsolute(amod_strip);
+                amod_strip.real = fabsf(amod_strip.real);
                 mod_strip = cadd(mod_strip, amod_strip);
             }
-        
+            //printf("modstrip: %f %f\n", mod_strip.real, mod_strip.imag);
+
             /* loop filter made up of 1st order IIR plus integrator.  Integerator
                was found to be reqd  */
         
             fdmdv->filt = (1.0-beta)*fdmdv->filt + beta*atan2(mod_strip.imag, mod_strip.real);
+            //printf("filt: %f angle: %f\n", fdmdv->filt, atan2(mod_strip.imag, mod_strip.real));
             *f_est += g*fdmdv->filt;
         }    
 
@@ -663,7 +654,7 @@ void cohpsk_demod(struct COHPSK *coh, int rx_bits[], int *reliable_sync_bit, COM
         max_ratio = 0.0;
         for (coh->f_est = FDMDV_FCENTRE-40.0; coh->f_est <= FDMDV_FCENTRE+40.0; coh->f_est += 40.0) {
         
-            printf("  [%d] acohpsk.f_est: %f +/- 20\n", coh->frame, coh->f_est);
+            fprintf(stderr, "  [%d] acohpsk.f_est: %f +/- 20\n", coh->frame, coh->f_est);
 
             /* we are out of sync so reset f_est and process two frames to clean out memories */
 
@@ -728,7 +719,7 @@ void cohpsk_demod(struct COHPSK *coh, int rx_bits[], int *reliable_sync_bit, COM
     /* If in sync just do sample rate processing on latest frame */
 
     if (sync == 1) {
-        rate_Fs_rx_processing(coh, ch_symb, rx_fdm, &coh->f_est, NSYMROWPILOT, nin, 0);
+        rate_Fs_rx_processing(coh, ch_symb, rx_fdm, &coh->f_est, NSYMROWPILOT, nin, 1);
         frame_sync_fine_freq_est(coh, ch_symb, sync, &next_sync);
  
         for(r=0; r<2; r++)

diff --git a/codec2-dev/unittest/tcohpsk.c b/codec2-dev/unittest/tcohpsk.c
index c4990062ad7b35c75f3628d8ad96994bfad314be..71d1109ca69a554fa25a64b50000cd40dd647c5e 100644 (file)
--- a/codec2-dev/unittest/tcohpsk.c
+++ b/codec2-dev/unittest/tcohpsk.c
@@ -49,7 +49,8 @@
 #define FRAMESL     (SYNC_FRAMES*FRAMES)  /* worst case is every frame is out of sync                           */
 
 #define RS          50
-#define FOFF        0
+#define FOFF        55.5
+#define DFOFF       (-0.5/(float)FS)
 #define ESNODB      8
 
 extern float pilots_coh[][PILOTS_NC];
@@ -80,7 +81,8 @@ int main(int argc, char *argv[])
     FILE          *fout;
     int            f, r, c, log_r, log_data_r, noise_r, ff_log_r;
     int           *ptest_bits_coh, *ptest_bits_coh_end;
-    COMP           phase_ch;
+    double         foff;
+    COMP           foff_rect, phase_ch;
 
     struct FDMDV  *fdmdv;
     //COMP           rx_filt[COHPSK_NC*ND][P+1];
@@ -91,6 +93,8 @@ int main(int argc, char *argv[])
     //COMP           rx_onesym[COHPSK_NC*ND];
     //int            rx_baseband_log_col_index = 0;
     //COMP           rx_baseband_log[COHPSK_NC*ND][(M+M/P)*NSYMROWPILOT*FRAMES];
+    float            f_est_log[FRAMES];
+    int              f_est_samples;
 
     int            log_bits;
     float          EsNo, variance;
@@ -114,12 +118,13 @@ int main(int argc, char *argv[])
     
     /* init stuff */
 
-    log_r = log_data_r = noise_r = log_bits = ff_log_r = 0;
+    log_r = log_data_r = noise_r = log_bits = ff_log_r = f_est_samples = 0;
     ptest_bits_coh = (int*)test_bits_coh;
     ptest_bits_coh_end = (int*)test_bits_coh + sizeof(test_bits_coh)/sizeof(int);
     memcpy(tx_bits, test_bits_coh, sizeof(int)*COHPSK_BITS_PER_FRAME);
 
     phase_ch.real = 1.0; phase_ch.imag = 0.0; 
+    foff = FOFF;
      
     /*  each carrier has power = 2, total power 2Nc, total symbol rate
         NcRs, noise BW B=Fs Es/No = (C/Rs)/(N/B), N = var =
@@ -127,6 +132,7 @@ int main(int argc, char *argv[])
 
     EsNo = pow(10.0, ESNODB/10.0);
     variance = 2.0*FS/(RS*EsNo);
+    //fprintf(stderr, "doff: %e\n", DFOFF);
 
     /* Main Loop ---------------------------------------------------------------------*/
 
@@ -156,8 +162,17 @@ int main(int argc, char *argv[])
                                  Channel
        \*---------------------------------------------------------*/
 
-        fdmdv_freq_shift(ch_fdm_frame, tx_fdm_frame, FOFF, &phase_ch, NSYMROWPILOT*M);
+        //fdmdv_freq_shift(ch_fdm_frame, tx_fdm_frame, FOFF, &phase_ch, NSYMROWPILOT*M);
 
+        for(r=0; r<NSYMROWPILOT*M; r++) {
+            foff_rect.real = cos(2.0*M_PI*foff/FS); foff_rect.imag = sin(2.0*M_PI*foff/FS);
+            foff += DFOFF;
+            phase_ch = cmult(phase_ch, foff_rect);
+            ch_fdm_frame[r] = cmult(tx_fdm_frame[r], phase_ch);
+        }
+        phase_ch.real /= cabsolute(phase_ch);
+        phase_ch.imag /= cabsolute(phase_ch);
+        //fprintf(stderr, "%f\n", foff);
         for(r=0; r<NSYMROWPILOT*M; r++,noise_r++) {
             scaled_noise = fcmult(sqrt(variance), noise[noise_r]);
             ch_fdm_frame[r] = cadd(ch_fdm_frame[r], scaled_noise);
@@ -201,13 +216,14 @@ int main(int argc, char *argv[])
             }
             memcpy(&rx_bits_log[COHPSK_BITS_PER_FRAME*log_bits], rx_bits, sizeof(int)*COHPSK_BITS_PER_FRAME);
             log_bits++;
+            f_est_log[f_est_samples++] = coh->f_est;
         }
 
        assert(log_r <= NSYMROWPILOT*FRAMES);
        assert(noise_r <= NSYMROWPILOT*M*FRAMES);
        assert(log_data_r <= NSYMROW*FRAMES);
 
-        printf("\r[%d]", f+1);
+        printf("\r  [%d]", f+1);
     }
     printf("\n");
 
@@ -232,6 +248,7 @@ int main(int argc, char *argv[])
     octave_save_float(fout, "rx_phi_log_c", (float*)rx_phi_log, log_data_r, COHPSK_NC*ND, COHPSK_NC*ND);  
     octave_save_complex(fout, "rx_symb_log_c", (COMP*)rx_symb_log, log_data_r, COHPSK_NC*ND, COHPSK_NC*ND);  
     octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1, COHPSK_BITS_PER_FRAME*log_bits);
+    octave_save_float(fout, "f_est_log_c", &f_est_log[1], 1, f_est_samples-1, f_est_samples-1);  
     fclose(fout);
 
     cohpsk_destroy(coh);