added mean_amp estimation to C, tofdm passes

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

diff --git a/codec2-dev/octave/ofdm_ldpc_rx.m b/codec2-dev/octave/ofdm_ldpc_rx.m
index de0de06e4f7e6301a2cfc0a98521c28a08453f80..90e056a479378a5d3f704a64d5573cca47413f1e 100644 (file)
--- a/codec2-dev/octave/ofdm_ldpc_rx.m
+++ b/codec2-dev/octave/ofdm_ldpc_rx.m
@@ -145,6 +145,8 @@ function ofdm_ldpc_rx(filename, interleave_frames = 1, error_pattern_filename)
       rx_amp(1:Nsymbolsperinterleavedframe-Nsymbolsperframe) = rx_amp(Nsymbolsperframe+1:Nsymbolsperinterleavedframe);
       rx_amp(Nsymbolsperinterleavedframe-Nsymbolsperframe+1:Nsymbolsperinterleavedframe) = arx_amp(Nuwtxtsymbolsperframe+1:end);
            
+      mean_amp = states.mean_amp;
+      
       % de-interleave QPSK symbols and symbol amplitudes
 
       rx_np_de = gp_deinterleave(rx_np);
@@ -159,7 +161,6 @@ function ofdm_ldpc_rx(filename, interleave_frames = 1, error_pattern_filename)
 
       if strcmp(states.sync_state_interleaver,'search')
         st = 1; en = Ncodedbitsperframe/bps;
-        mean_amp = states.mean_amp;
         [rx_codeword parity_checks] = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_np_de(st:en)/mean_amp, min(EsNo,30), rx_amp_de(st:en)/mean_amp);
         Nerrs = code_param.data_bits_per_frame - max(parity_checks);
         %printf("Nerrs: %d\n", Nerrs);

diff --git a/codec2-dev/octave/tofdm.m b/codec2-dev/octave/tofdm.m
index 1dc8d3dd13d1b7bd6beee5a8bc589477a640425d..d468c75d40420b8889674940972e5be185100f4e 100644 (file)
--- a/codec2-dev/octave/tofdm.m
+++ b/codec2-dev/octave/tofdm.m
@@ -59,6 +59,7 @@ if cml_support
   ibits = payload_data_bits;
   codeword = LdpcEncode(ibits, code_param.H_rows, code_param.P_matrix);
   tx_bits(Nuwbits+Ntxtbits+1:end) = codeword;
+  states.mean_amp = 1;  % start this with something sensible otherwise LDPC decode fails
 else
   tx_bits(Nuwbits+Ntxtbits+1:end) = [payload_data_bits payload_data_bits];
 end
@@ -92,7 +93,7 @@ timing_est_log = timing_valid_log = timing_mx_log = [];
 coarse_foff_est_hz_log = []; sample_point_log = [];
 phase_est_pilot_log = []; rx_amp_log = [];
 rx_np_log = []; rx_bits_log = [];
-sig_var_log = noise_var_log = [];
+sig_var_log = noise_var_log = mean_amp_log = [];
 
 states.timing_en = 1;
 states.foff_est_en = 1;
@@ -138,11 +139,14 @@ for f=1:Nframes
   rx_bits_log = [rx_bits_log rx_bits];
   sig_var_log = [sig_var_log; states.sig_var];
   noise_var_log = [noise_var_log; states.noise_var];
+  mean_amp_log = [mean_amp_log; states.mean_amp];
   
   % Optional testing of LDPC functions
 
   if cml_support
-    symbol_likelihood = Demod2D(arx_np(Nuwtxtsymbolsperframe+1:end), S_matrix, EsNo, arx_amp(Nuwtxtsymbolsperframe+1:end));
+    mean_amp = states.mean_amp;
+    %mean_amp = 1;
+    symbol_likelihood = Demod2D(arx_np(Nuwtxtsymbolsperframe+1:end)/mean_amp, S_matrix, EsNo, arx_amp(Nuwtxtsymbolsperframe+1:end)/mean_amp);
     bit_likelihood = Somap(symbol_likelihood);
 
     [x_hat paritychecks] = MpDecode(-bit_likelihood(1:code_param.code_bits_per_frame), code_param.H_rows, code_param.H_cols, max_iterations, decoder_type, 1, 1);
@@ -239,4 +243,5 @@ if cml_support
 end
 check(sig_var_log, sig_var_log_c, 'sig_var_log');
 check(noise_var_log, noise_var_log_c, 'noise_var_log');
+check(mean_amp_log, mean_amp_log_c, 'mean_amp_log');
 

diff --git a/codec2-dev/src/mpdecode_core.c b/codec2-dev/src/mpdecode_core.c
index 8b425a3dc98da16a546aa89e5bb6dcf50e5cbecb..e7afa4a7ccf871f8737d4931befc3feff7bcaa57 100644 (file)
--- a/codec2-dev/src/mpdecode_core.c
+++ b/codec2-dev/src/mpdecode_core.c
@@ -778,6 +778,7 @@ void Demod2D(double  symbol_likelihood[],       /* output, M*number_symbols
              COMP    S_matrix[],                /* constellation of size M               */
              float   EsNo,
              float   fading[],                  /* real fading values, number_symbols    */
+             float   mean_amp,                  
              int     number_symbols)
 {
     int     M=QPSK_CONSTELLATION_SIZE;
@@ -788,10 +789,10 @@ void Demod2D(double  symbol_likelihood[],       /* output, M*number_symbols
   
     for (i=0;i<number_symbols;i++) {                /* go through each received symbol */
         for (j=0;j<M;j++) {                         /* each postulated symbol          */
-            tempsr = fading[i]*S_matrix[j].real;
-            tempsi = fading[i]*S_matrix[j].imag;
-            Er = r[i].real - tempsr;
-            Ei = r[i].imag - tempsi;
+            tempsr = fading[i]*S_matrix[j].real/mean_amp;
+            tempsi = fading[i]*S_matrix[j].imag/mean_amp;
+            Er = r[i].real/mean_amp - tempsr;
+            Ei = r[i].imag/mean_amp - tempsi;
             symbol_likelihood[i*M+j] = -EsNo*(Er*Er+Ei*Ei);
             //printf("symbol_likelihood[%d][%d] = %f\n", i,j,symbol_likelihood[i*M+j]);
         }
@@ -871,13 +872,13 @@ int extract_output(char out_char[], int DecodedBits[], int ParityCheckCount[], i
     return iter;
 }
 
-void symbols_to_llrs(double llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, int nsyms) {
+void symbols_to_llrs(double llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, float mean_amp, int nsyms) {
     int i;
     
     double symbol_likelihood[nsyms*QPSK_CONSTELLATION_SIZE];
     double bit_likelihood[nsyms*QPSK_BITS_PER_SYMBOL];
 
-    Demod2D(symbol_likelihood, rx_qpsk_symbols, S_matrix, EsNo, rx_amps,  nsyms);
+    Demod2D(symbol_likelihood, rx_qpsk_symbols, S_matrix, EsNo, rx_amps, mean_amp, nsyms);
     Somap(bit_likelihood, symbol_likelihood, nsyms);
     for(i=0; i<nsyms*QPSK_BITS_PER_SYMBOL; i++) {
         llr[i] = -bit_likelihood[i];

diff --git a/codec2-dev/src/mpdecode_core.h b/codec2-dev/src/mpdecode_core.h
index e0ba45a74a5a73ec6373ba77e6a63c32b9016185..50a2ad359564a76c03c7f0b9c7a7feb735a2f56a 100644 (file)
--- a/codec2-dev/src/mpdecode_core.h
+++ b/codec2-dev/src/mpdecode_core.h
@@ -36,9 +36,9 @@ int run_ldpc_decoder(struct LDPC *ldpc, char out_char[], double input[], int *pa
 
 void sd_to_llr(double llr[], double sd[], int n);
 
-void Demod2D(double symbol_likelihood[], COMP r[], COMP S_matrix[], float EsNo, float fading[], int number_symbols);
+void Demod2D(double symbol_likelihood[], COMP r[], COMP S_matrix[], float EsNo, float fading[], float mean_amp, int number_symbols);
 void Somap(double bit_likelihood[], double symbol_likelihood[], int number_symbols);
-void symbols_to_llrs(double llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, int nsyms);
+void symbols_to_llrs(double llr[], COMP rx_qpsk_symbols[], float rx_amps[], float EsNo, float mean_amp, int nsyms);
 
 struct v_node {
   int degree;

diff --git a/codec2-dev/src/ofdm.c b/codec2-dev/src/ofdm.c
index 7ee14b1acaf19fbf19ed836f429fbdccb64079a2..528175c8ffce31fe0ef8c7ac33a2d17fb2d52359 100644 (file)
--- a/codec2-dev/src/ofdm.c
+++ b/codec2-dev/src/ofdm.c
@@ -874,9 +874,9 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
         aphase_est_pilot_rect = aphase_est_pilot_rect + vector_sum(symbol, 3);
         aphase_est_pilot[i] = cargf(aphase_est_pilot_rect);
 
-        /* amplitude is estimated over 6 rows of pilots */
+        /* amplitude is estimated over 12 pilots */
 
-        aamp_est_pilot[i] = cabsf(aphase_est_pilot_rect / 6.0f);
+        aamp_est_pilot[i] = cabsf(aphase_est_pilot_rect / 12.0f);
     }
 
     /*
@@ -888,7 +888,8 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
     complex float rx_corr;
     int abit[2];
     int bit_index = 0;
-
+    float sum_amp = 0.0;
+    
     for (rr = 0; rr < OFDM_ROWSPERFRAME; rr++) {
         /*
          * Note the i has an index of 1 to 16, so we ignore carriers 0 and 17.
@@ -917,7 +918,8 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
              */
 
             ofdm->rx_amp[(rr * OFDM_NC) + (i - 1)] = aamp_est_pilot[i];
-
+            sum_amp += aamp_est_pilot[i];
+            
             /*
              * Note like amps in this implementation phase ests are the
              * same for each col, but we log them for each symbol anyway
@@ -939,6 +941,10 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
         }
     }
 
+    /* update mean amplitude estimate for LDPC decoder scaling */
+    
+    ofdm->mean_amp = 0.9*ofdm->mean_amp + 0.1*sum_amp/(OFDM_ROWSPERFRAME * OFDM_NC);
+    
     /* Adjust nin to take care of sample clock offset */
 
     ofdm->nin = OFDM_SAMPLESPERFRAME;
@@ -957,6 +963,37 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
             ofdm->sample_point += tshift;
         }
     }
+
+    /* estimate signal and noise power, see ofdm_lib.m, cohpsk.m for
+       more info */
+
+    float sig_var = 0.0;
+    complex float *rx_np = ofdm->rx_np;
+    for(i=0; i<OFDM_ROWSPERFRAME * OFDM_NC; i++) {
+        sig_var += crealf(rx_np[i])*crealf(rx_np[i]) + cimagf(rx_np[i])*cimagf(rx_np[i]);
+    }
+    sig_var /= (OFDM_ROWSPERFRAME * OFDM_NC);
+    float sig_rms = sqrtf(sig_var);
+
+    complex float s;
+    float sum_x = 0;
+    float sum_xx = 0;
+    int n = 0;
+    for (i=0; i<OFDM_ROWSPERFRAME * OFDM_NC; i++) {
+      s = rx_np[i];
+      if (fabs(crealf(s)) > sig_rms) {
+          sum_x  += cimagf(s);
+          sum_xx += cimagf(s) * cimagf(s);
+          n++;
+      }
+    }
+
+    float noise_var = 0;
+    if (n > 1) {
+      noise_var = (n*sum_xx - sum_x*sum_x)/(n*(n-1));
+    }
+    ofdm->sig_var = sig_var;
+    ofdm->noise_var = 2*noise_var;
 }
 
 

diff --git a/codec2-dev/src/ofdm_internal.h b/codec2-dev/src/ofdm_internal.h
index c4b30d5e6e4d8ae5eab1c12d833ef41d7a124899..5916e0845bff876e7e7841370b29fb3f0a037b92 100644 (file)
--- a/codec2-dev/src/ofdm_internal.h
+++ b/codec2-dev/src/ofdm_internal.h
@@ -124,7 +124,10 @@ struct OFDM {
     complex float rx_np[OFDM_ROWSPERFRAME * OFDM_NC];
     float rx_amp[OFDM_ROWSPERFRAME * OFDM_NC];
     float aphase_est_pilot_log[OFDM_ROWSPERFRAME * OFDM_NC];
-
+    float sig_var;
+    float noise_var;
+    float mean_amp;
+    
     /* modem sync state machine */
 
     int  tx_uw[OFDM_NUWBITS];

diff --git a/codec2-dev/unittest/tofdm.c b/codec2-dev/unittest/tofdm.c
index a471d7f14beab4e80be6594f7eb7bb04f1d0e2f4..a56cb72e788a677a46b3f497924ca6147fe1ed6e 100644 (file)
--- a/codec2-dev/unittest/tofdm.c
+++ b/codec2-dev/unittest/tofdm.c
@@ -165,6 +165,8 @@ int main(int argc, char *argv[])
     float          symbol_likelihood_log[ (CODED_BITSPERFRAME/OFDM_BPS) * (1<<OFDM_BPS) * NFRAMES];
     float          bit_likelihood_log[CODED_BITSPERFRAME * NFRAMES];        
     int            detected_data_log[CODED_BITSPERFRAME * NFRAMES];
+    float          sig_var_log[NFRAMES], noise_var_log[NFRAMES];        
+    float          mean_amp_log[NFRAMES];        
     
     FILE          *fout;
     int            f,i,j;
@@ -293,7 +295,11 @@ int main(int argc, char *argv[])
     int Nmaxsamperframe = ofdm_get_max_samples_per_frame();
     short rx_scaled[Nmaxsamperframe];
     #endif
-    
+
+    /* start this with something sensible otherwise LDPC decode fails in tofdm.m */
+
+    ofdm->mean_amp = 1.0;
+       
     for(f=0; f<NFRAMES; f++) {
         /* For initial testing, timing est is off, so nin is always
            fixed.  TODO: we need a constant for rxbuf_in[] size that
@@ -363,7 +369,7 @@ int main(int argc, char *argv[])
         }
         float *ldpc_codeword_symbol_amps = &ofdm->rx_amp[(OFDM_NUWBITS+OFDM_NTXTBITS)/OFDM_BPS];
                 
-        Demod2D(symbol_likelihood, ldpc_codeword_symbols, S_matrix, EsNo, ldpc_codeword_symbol_amps,  CODED_BITSPERFRAME/OFDM_BPS);
+        Demod2D(symbol_likelihood, ldpc_codeword_symbols, S_matrix, EsNo, ldpc_codeword_symbol_amps, ofdm->mean_amp, CODED_BITSPERFRAME/OFDM_BPS);
         Somap(bit_likelihood, symbol_likelihood, CODED_BITSPERFRAME/OFDM_BPS);
 
         int    iter;
@@ -424,10 +430,13 @@ int main(int argc, char *argv[])
 
         foff_hz_log[f] = ofdm->foff_est_hz;
         timing_est_log[f] = ofdm->timing_est + 1;     /* offset by 1 to match Octave */
-        timing_valid_log[f] = ofdm->timing_valid;     /* offset by 1 to match Octave */
-        timing_mx_log[f] = ofdm->timing_mx;           /* offset by 1 to match Octave */
+        timing_valid_log[f] = ofdm->timing_valid;     
+        timing_mx_log[f] = ofdm->timing_mx;           
         coarse_foff_est_hz_log[f] = ofdm->coarse_foff_est_hz;
         sample_point_log[f] = ofdm->sample_point + 1; /* offset by 1 to match Octave */
+        sig_var_log[f] = ofdm->sig_var;
+        noise_var_log[f] = ofdm->noise_var;
+        mean_amp_log[f] = ofdm->mean_amp;
 
         memcpy(&rx_bits_log[OFDM_BITSPERFRAME*f], rx_bits, sizeof(rx_bits));
 
@@ -469,6 +478,9 @@ int main(int argc, char *argv[])
     octave_save_float(fout, "symbol_likelihood_log_c", symbol_likelihood_log, (CODED_BITSPERFRAME/OFDM_BPS) * (1<<OFDM_BPS) * NFRAMES, 1, 1);
     octave_save_float(fout, "bit_likelihood_log_c", bit_likelihood_log, CODED_BITSPERFRAME * NFRAMES, 1, 1);
     octave_save_int(fout, "detected_data_log_c", detected_data_log, 1, CODED_BITSPERFRAME*NFRAMES);
+    octave_save_float(fout, "sig_var_log_c", sig_var_log, NFRAMES, 1, 1);
+    octave_save_float(fout, "noise_var_log_c", noise_var_log, NFRAMES, 1, 1);
+    octave_save_float(fout, "mean_amp_log_c", mean_amp_log, NFRAMES, 1, 1);
     fclose(fout);
 
     ofdm_destroy(ofdm);