[Commit-gnuradio] [gnuradio] 13/24: digital: Dealing with scaling issues with corr_est.

[git]

This is an automated email from the git hooks/post-receive script.

jcorgan pushed a commit to branch packet2
in repository gnuradio.

commit 5d1f389b571b8ad5674fbe328c1f86b4e05b63e2
Author: Tom Rondeau <address@hidden>
Date:   Tue Apr 12 16:25:19 2016 -0400

    digital: Dealing with scaling issues with corr_est.
    
    The correlation estimator detector assumes the input is scaled to near
    +/-1. This is not a guaranteed condition and AGCs can be slow to
    handle the bursts well.
    
    Instead, this plays with a threshold based on the second moment (more
    or less) of the recieved samples. Anything too far outside this range
    is detected as a (potential) peak.
    
    The second thing that this commit does is add a new "amp_est"
    tag. This is the reciprocal of the amplitude of the sample when the
    correlation peak is detected. It is meant so that a downstream block
    can use this tag to rescale all samples by the value of this tag.
---
 gr-digital/include/gnuradio/digital/corr_est_cc.h |  8 +--
 gr-digital/lib/corr_est_cc_impl.cc                | 73 +++++++++++++----------
 gr-digital/lib/corr_est_cc_impl.h                 |  3 +
 3 files changed, 48 insertions(+), 36 deletions(-)

diff --git a/gr-digital/include/gnuradio/digital/corr_est_cc.h 
b/gr-digital/include/gnuradio/digital/corr_est_cc.h
index e8211cf..237a499 100644
--- a/gr-digital/include/gnuradio/digital/corr_est_cc.h
+++ b/gr-digital/include/gnuradio/digital/corr_est_cc.h
@@ -42,6 +42,7 @@ namespace gr {
      * \li tag 'phase_est': estimate of phase offset
      * \li tag 'time_est': estimate of symbol timing offset
      * \li tag 'corr_est': the correlation value of the estimates
+     * \li tag 'amp_est': 1 over the estimated amplitude
      * \li tag 'corr_start': the start sample of the correlation and the value
      *
      * \li Optional 2nd output stream providing the advanced correlator output
@@ -51,12 +52,11 @@ namespace gr {
      * offset estimate. These estimates are passed downstream as
      * stream tags for use by follow-on synchronization blocks.
      *
-     * The sync word is provided as a set of symbols along with a
-     * baseband matched filter which we use to create the filtered and
-     * upsampled symbols that we will receive over-the-air.
+     * The sync word is provided as a set of symbols after being
+     * filtered by a baseband matched filter.
      *
      * The phase_est tag can be used by downstream blocks to adjust
-     * their phase estimatopm/correction loops, and is currently
+     * their phase estimator/correction loops, and is currently
      * implemented by the gr::digital::costas_loop_cc block.
      *
      * The time_est tag can be used to adjust the sampling timing
diff --git a/gr-digital/lib/corr_est_cc_impl.cc 
b/gr-digital/lib/corr_est_cc_impl.cc
index 25e6b7f..dfcd966 100644
--- a/gr-digital/lib/corr_est_cc_impl.cc
+++ b/gr-digital/lib/corr_est_cc_impl.cc
@@ -55,6 +55,17 @@ namespace gr {
     {
       d_sps = sps;
 
+      // In order to easily support the optional second output,
+      // don't deal with an unbounded max number of output items.
+      // For the common case of not using the optional second output,
+      // this ensures we optimally call the volk routines.
+      const size_t nitems = 24*1024;
+      set_max_noutput_items(nitems);
+      d_corr = (gr_complex *)
+               volk_malloc(sizeof(gr_complex)*nitems, volk_get_alignment());
+      d_corr_mag = (float *)
+                   volk_malloc(sizeof(float)*nitems, volk_get_alignment());
+
       // Create time-reversed conjugate of symbols
       d_symbols = symbols;
       for(size_t i=0; i < d_symbols.size(); i++) {
@@ -96,16 +107,7 @@ namespace gr {
       //  volk_get_alignment() / sizeof(gr_complex);
       //set_alignment(std::max(1,alignment_multiple));
 
-      // In order to easily support the optional second output,
-      // don't deal with an unbounded max number of output items.
-      // For the common case of not using the optional second output,
-      // this ensures we optimally call the volk routines.
-      const size_t nitems = 24*1024;
-      set_max_noutput_items(nitems);
-      d_corr = (gr_complex *)
-               volk_malloc(sizeof(gr_complex)*nitems, volk_get_alignment());
-      d_corr_mag = (float *)
-                   volk_malloc(sizeof(float)*nitems, volk_get_alignment());
+      d_scale = 1.0f;
     }
 
     corr_est_cc_impl::~corr_est_cc_impl()
@@ -193,17 +195,7 @@ namespace gr {
     corr_est_cc_impl::_set_threshold(float threshold)
     {
       d_stashed_threshold = threshold;
-
-      // Compute a correlation threshold.
-      // Compute the value of the discrete autocorrelation of the matched
-      // filter with offset 0 (aka the autocorrelation peak).
-      float corr = 0;
-      for(size_t i = 0; i < d_symbols.size(); i++)
-        corr += abs(d_symbols[i]*conj(d_symbols[i]));
-
-      // Half the user-specified threshold to account for summing two
-      // consecutive correlation values in the work function.
-      d_thresh = (threshold + (1-threshold)/2.0f)*corr*corr;
+      d_pfa = -logf(1.0f-threshold);
     }
 
     void
@@ -231,10 +223,6 @@ namespace gr {
       // Our correlation filter length
       unsigned int hist_len = history() - 1;
 
-      // Delay the output by our correlation filter length so we can
-      // tag backwards in time
-      memcpy(out, &in[0], sizeof(gr_complex)*noutput_items);
-
       // Calculate the correlation of the non-delayed input with the
       // known symbols.
       d_filter->filter(noutput_items, &in[hist_len], corr);
@@ -242,24 +230,31 @@ namespace gr {
       // Find the magnitude squared of the correlation
       volk_32fc_magnitude_squared_32f(&d_corr_mag[0], corr, noutput_items);
 
+      float detection = 0;
+      for(int i = 0; i < noutput_items; i++) {
+        detection += d_corr_mag[i];
+      }
+      detection /= static_cast<float>(noutput_items);
+      detection *= d_pfa;
+
       int isps = (int)(d_sps + 0.5f);
       int i = 0;
       while(i < noutput_items) {
-        // Look for the correlator output to cross the threshold Sum
-        // power over two consecutive symbols in case we're offset in
-        // time. If off by 1/2 a symbol, the peak of any one point is
-        // much lower.
+        // Look for the correlator output to cross the threshold.
+        // Sum power over two consecutive symbols in case we're offset
+        // in time. If off by 1/2 a symbol, the peak of any one point
+        // is much lower.
         float corr_mag = d_corr_mag[i] + d_corr_mag[i+1];
-        if (corr_mag <= d_thresh) {
+        if(corr_mag <= 4*detection) {
           i++;
           continue;
         }
 
         // Go to (just past) the current correlator output peak
         while ((i < (noutput_items-1)) &&
-               (d_corr_mag[i] < d_corr_mag[i+1]))
+               (d_corr_mag[i] < d_corr_mag[i+1])) {
           i++;
-
+        }
         // Delaying the primary signal output by the matched filter
         // length using history(), means that the the peak output of
         // the matched filter aligns with the start of the desired
@@ -294,6 +289,12 @@ namespace gr {
         center = (center - 2.0); // adjust for bias in center of mass 
calculation
 #endif
 
+        // Estimated scaling factor for the input stream to normalize
+        // the output to +/-1.
+        uint32_t maxi;
+        volk_32fc_index_max_16u(&maxi, (gr_complex*)in, noutput_items);
+        d_scale = 1 / std::abs(in[maxi]);
+
         // Calculate the phase offset of the incoming signal.
         //
         // The analytic cross-correlation is:
@@ -322,6 +323,8 @@ namespace gr {
         // N.B. the appropriate d_corr_mag[] index is "i", not "index".
         add_item_tag(0, nitems_written(0) + index, pmt::intern("corr_est"),
                      pmt::from_double(d_corr_mag[i]), d_src_id);
+        add_item_tag(0, nitems_written(0) + index, pmt::intern("amp_est"),
+                     pmt::from_double(d_scale), d_src_id);
 
         if (output_items.size() > 1) {
           // N.B. these debug tags are not offset to avoid walking off out buf
@@ -331,6 +334,8 @@ namespace gr {
                        pmt::from_double(center), d_src_id);
           add_item_tag(1, nitems_written(0) + i, pmt::intern("corr_est"),
                        pmt::from_double(d_corr_mag[i]), d_src_id);
+          add_item_tag(1, nitems_written(0) + i, pmt::intern("amp_est"),
+                       pmt::from_double(d_scale), d_src_id);
         }
 
         // Skip ahead to the next potential symbol peak
@@ -343,6 +348,10 @@ namespace gr {
       //               pmt::intern("ce_eow"), pmt::from_uint64(noutput_items),
       //               d_src_id);
 
+      // Delay the output by our correlation filter length so we can
+      // tag backwards in time
+      memcpy(out, &in[0], sizeof(gr_complex)*noutput_items);
+
       return noutput_items;
     }
 
diff --git a/gr-digital/lib/corr_est_cc_impl.h 
b/gr-digital/lib/corr_est_cc_impl.h
index 6e8dd17..e0d4c8e 100644
--- a/gr-digital/lib/corr_est_cc_impl.h
+++ b/gr-digital/lib/corr_est_cc_impl.h
@@ -44,6 +44,9 @@ namespace gr {
       gr_complex *d_corr;
       float *d_corr_mag;
 
+      float d_scale;
+      float d_pfa; // probability of false alarm
+
       void _set_mark_delay(unsigned int mark_delay);
       void _set_threshold(float threshold);