4 files changed, 127 insertions, 86 deletions
diff --git a/gr-digital/examples/snr_estimators.py b/gr-digital/examples/snr_estimators.py
index 9eae7865e5..31efe6b83e 100755
--- a/gr-digital/examples/snr_estimators.py
+++ b/gr-digital/examples/snr_estimators.py
@@ -1,24 +1,24 @@
 #!/usr/bin/env python
 #
 # Copyright 2011-2013 Free Software Foundation, Inc.
-# 
+#
 # This file is part of GNU Radio
-# 
+#
 # GNU Radio is free software; you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
 # the Free Software Foundation; either version 3, or (at your option)
 # any later version.
-# 
+#
 # GNU Radio is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
-# 
+#
 # You should have received a copy of the GNU General Public License
 # along with GNU Radio; see the file COPYING.  If not, write to
 # the Free Software Foundation, Inc., 51 Franklin Street,
 # Boston, MA 02110-1301, USA.
-# 
+#
 
 import sys
 
@@ -34,7 +34,7 @@ try:
 except ImportError:
     print "Error: Program requires Matplotlib (matplotlib.sourceforge.net)."
     sys.exit(1)
-    
+
 from gnuradio import gr, digital, filter
 from gnuradio import blocks
 from gnuradio import channels
@@ -49,45 +49,45 @@ For an explination of the online algorithms, see:
 http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
 '''
 
-def online_skewness(data, alpha):
+def online_skewness(data):
     n = 0
     mean = 0
     M2 = 0
     M3 = 0
-    d_M3 = 0
 
     for n in xrange(len(data)):
         delta = data[n] - mean
         delta_n = delta / (n+1)
-        term1 = delta * delta_n * (n)
+        term1 = delta * delta_n * n
         mean = mean + delta_n
-        M3 = term1 * delta_n * (n - 1) - 3 * delta_n * M2
+        M3 = M3 + term1 * delta_n * (n - 1) - 3 * delta_n * M2
         M2 = M2 + term1
-        d_M3 = (0.001)*M3 + (1-0.001)*d_M3;
-        
-    return d_M3
+
+    return scipy.sqrt(len(data))*M3 / scipy.power(M2, 3.0/2.0);
 
 def snr_est_simple(signal):
-    y1 = scipy.mean(abs(signal))
-    y2 = scipy.real(scipy.mean(signal**2))
-    y3 = (y1*y1 - y2)
-    snr_rat = y1*y1/y3
+    s = scipy.mean(abs(signal)**2)
+    n = 2*scipy.var(abs(signal))
+    snr_rat = s/n
     return 10.0*scipy.log10(snr_rat), snr_rat
-    
+
 def snr_est_skew(signal):
     y1 = scipy.mean(abs(signal))
     y2 = scipy.mean(scipy.real(signal**2))
     y3 = (y1*y1 - y2)
-    y4 = online_skewness(abs(signal.real), 0.001)
+    y4 = online_skewness(signal.real)
+    #y4 = stats.skew(abs(signal.real))
 
     skw = y4*y4 / (y2*y2*y2);
-    snr_rat = y1*y1 / (y3 + skw*y1*y1)
+    s = y1*y1
+    n = 2*(y3 + skw*s)
+    snr_rat = s / n
     return 10.0*scipy.log10(snr_rat), snr_rat
 
 def snr_est_m2m4(signal):
     M2 = scipy.mean(abs(signal)**2)
     M4 = scipy.mean(abs(signal)**4)
-    snr_rat = 2*scipy.sqrt(2*M2*M2 - M4) / (M2 - scipy.sqrt(2*M2*M2 - M4))
+    snr_rat = scipy.sqrt(2*M2*M2 - M4) / (M2 - scipy.sqrt(2*M2*M2 - M4))
     return 10.0*scipy.log10(snr_rat), snr_rat
 
 def snr_est_svr(signal):
@@ -100,10 +100,10 @@ def snr_est_svr(signal):
     savg = (1.0/(float(N)-1.0))*ssum
     mavg = (1.0/(float(N)-1.0))*msum
     beta = savg / (mavg - savg)
-    
-    snr_rat = 2*((beta - 1) + scipy.sqrt(beta*(beta-1)))
+
+    snr_rat = ((beta - 1) + scipy.sqrt(beta*(beta-1)))
     return 10.0*scipy.log10(snr_rat), snr_rat
-    
+
 
 def main():
     gr_estimators = {"simple": digital.SNR_EST_SIMPLE,
@@ -114,8 +114,8 @@ def main():
                      "skew": snr_est_skew,
                      "m2m4": snr_est_m2m4,
                      "svr": snr_est_svr}
-    
-    
+
+
     parser = OptionParser(option_class=eng_option, conflict_handler="resolve")
     parser.add_option("-N", "--nsamples", type="int", default=10000,
                       help="Set the number of samples to process [default=%default]")
@@ -134,12 +134,14 @@ def main():
     N = options.nsamples
     xx = scipy.random.randn(N)
     xy = scipy.random.randn(N)
-    bits = 2*scipy.complex64(scipy.random.randint(0, 2, N)) - 1
+    bits =2*scipy.complex64(scipy.random.randint(0, 2, N)) - 1
+    #bits =(2*scipy.complex64(scipy.random.randint(0, 2, N)) - 1) + \
+    #    1j*(2*scipy.complex64(scipy.random.randint(0, 2, N)) - 1)
 
     snr_known = list()
     snr_python = list()
     snr_gr = list()
-    
+
     # when to issue an SNR tag; can be ignored in this example.
     ntag = 10000
 
@@ -154,17 +156,17 @@ def main():
     SNR_dB = scipy.arange(SNR_min, SNR_max+SNR_step, SNR_step)
     for snr in SNR_dB:
         SNR = 10.0**(snr/10.0)
-        scale = scipy.sqrt(SNR)
+        scale = scipy.sqrt(2*SNR)
         yy = bits + n_cpx/scale
         print "SNR: ", snr
 
         Sknown = scipy.mean(yy**2)
-        Nknown = scipy.var(n_cpx/scale)/2
+        Nknown = scipy.var(n_cpx/scale)
         snr0 = Sknown/Nknown
         snr0dB = 10.0*scipy.log10(snr0)
-        snr_known.append(snr0dB)
+        snr_known.append(float(snr0dB))
 
-        snrdB, snr = py_est(yy)        
+        snrdB, snr = py_est(yy)
         snr_python.append(snrdB)
 
         gr_src = blocks.vector_source_c(bits.tolist(), False)
@@ -188,9 +190,12 @@ def main():
     s1.set_ylabel('Estimated SNR')
     s1.legend()
 
+    f2 = pylab.figure(2)
+    s2 = f2.add_subplot(1,1,1)
+    s2.plot(yy.real, yy.imag, 'o')
+
     pylab.show()
 
 
 if __name__ == "__main__":
     main()
-    
diff --git a/gr-digital/include/gnuradio/digital/mpsk_snr_est.h b/gr-digital/include/gnuradio/digital/mpsk_snr_est.h
index 46df06164d..5398745383 100644
--- a/gr-digital/include/gnuradio/digital/mpsk_snr_est.h
+++ b/gr-digital/include/gnuradio/digital/mpsk_snr_est.h
@@ -1,19 +1,19 @@
 /* -*- c++ -*- */
 /*
  * Copyright 2011,2012 Free Software Foundation, Inc.
- * 
+ *
  * This file is part of GNU Radio
- * 
+ *
  * GNU Radio is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation; either version 3, or (at your option)
  * any later version.
- * 
+ *
  * GNU Radio is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
  * along with GNU Radio; see the file COPYING.  If not, write to
  * the Free Software Foundation, Inc., 51 Franklin Street,
@@ -58,6 +58,7 @@ namespace gr {
     {
     protected:
       double d_alpha, d_beta;
+      double d_signal, d_noise;
 
     public:
       /*! Constructor
@@ -81,6 +82,12 @@ namespace gr {
 
       //! Use the register values to compute a new estimate
       virtual double snr();
+
+      //! Returns the signal power estimate
+      virtual double signal();
+
+      //! Returns the noise power estimate
+      virtual double noise();
     };
 
 
@@ -97,7 +104,8 @@ namespace gr {
     {
     private:
       double d_y1, d_y2;
-  
+      double d_counter;
+
     public:
       /*! Constructor
        *
@@ -123,13 +131,16 @@ namespace gr {
      *  affected because of fold-over around the decision boundaries,
      *  which results in a skewness to the samples. We estimate the
      *  skewness and use this as a correcting term.
+     *
+     *  This algorithm only appears to work well for BPSK signals.
      */
     class DIGITAL_API mpsk_snr_est_skew :
       public mpsk_snr_est
     {
     private:
       double d_y1, d_y2, d_y3;
-  
+      double d_counter;
+
     public:
       /*! Constructor
        *
@@ -167,7 +178,7 @@ namespace gr {
     {
     private:
       double d_y1, d_y2;
-  
+
     public:
       /*! Constructor
        *
@@ -205,7 +216,7 @@ namespace gr {
      *  estimator unless you have a way to guess or estimate these
      *  values here.
      *
-     *  Original paper: 
+     *  Original paper:
      *  R. Matzner, "An SNR estimation algorithm for complex baseband
      *  signal using higher order statistics," Facta Universitatis
      *  (Nis), no. 6, pp. 41-52, 1993.
@@ -221,7 +232,7 @@ namespace gr {
     private:
       double d_y1, d_y2;
       double d_ka, d_kw;
-  
+
     public:
       /*! Constructor
        *
@@ -266,7 +277,7 @@ namespace gr {
     {
     private:
       double d_y1, d_y2;
-  
+
     public:
       /*! Constructor
        *
diff --git a/gr-digital/lib/mpsk_snr_est.cc b/gr-digital/lib/mpsk_snr_est.cc
index 6edb0274f5..4e7825b894 100644
--- a/gr-digital/lib/mpsk_snr_est.cc
+++ b/gr-digital/lib/mpsk_snr_est.cc
@@ -1,19 +1,19 @@
 /* -*- c++ -*- */
 /*
  * Copyright 2011,2012 Free Software Foundation, Inc.
- * 
+ *
  * This file is part of GNU Radio
- * 
+ *
  * GNU Radio is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation; either version 3, or (at your option)
  * any later version.
- * 
+ *
  * GNU Radio is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
  * along with GNU Radio; see the file COPYING.  If not, write to
  * the Free Software Foundation, Inc., 51 Franklin Street,
@@ -33,6 +33,8 @@ namespace gr {
 
     mpsk_snr_est::mpsk_snr_est(double alpha)
     {
+      d_signal = 0;
+      d_noise = 0;
       set_alpha(alpha);
     }
 
@@ -65,6 +67,18 @@ namespace gr {
       throw std::runtime_error("mpsk_snr_est: Unimplemented");
     }
 
+    double
+    mpsk_snr_est::signal()
+    {
+      return 10.0*log10(d_signal);
+    }
+
+    double
+    mpsk_snr_est::noise()
+    {
+      return 10.0*log10(d_noise);
+    }
+
 
     /*****************************************************************/
 
@@ -74,18 +88,27 @@ namespace gr {
     {
       d_y1 = 0;
       d_y2 = 0;
+      d_counter = 0;
     }
 
     int
     mpsk_snr_est_simple::update(int noutput_items,
 				const gr_complex *input)
     {
-      for(int i = 0; i < noutput_items; i++) {
-	double y1 = abs(input[i]);
-	d_y1 = d_alpha*y1 + d_beta*d_y1;
+      int i = 0;
+      if(d_counter == 0) {
+        d_y1 = abs(input[0]);
+        d_y2 = 0;
+        d_counter++;
+        i++;
+      }
 
-	double y2 = real(input[i]*input[i]);
-	d_y2 = d_alpha*y2 + d_beta*d_y2;
+      for(; i < noutput_items; i++) {
+        double x = abs(input[i]);
+        double m = d_y1 + (x - d_y1)/d_counter;
+        d_y2 = d_y2 + (x - d_y1)*(x - m);
+        d_y1 = m;
+        d_counter++;
       }
       return noutput_items;
     }
@@ -93,9 +116,10 @@ namespace gr {
     double
     mpsk_snr_est_simple::snr()
     {
-      double y1_2 = d_y1*d_y1;
-      double y3 = y1_2 - d_y2 + 1e-20;
-      return 10.0*log10(y1_2/y3);
+      d_signal = d_y1;
+      d_noise = 2 * d_y2 / (d_counter - 1);
+
+      return 10.0*log10(d_signal/d_noise);
     }
 
 
@@ -107,7 +131,7 @@ namespace gr {
     {
       d_y1 = 0;
       d_y2 = 0;
-      d_y3 = 0;
+      d_counter = 1;
     }
 
     int
@@ -118,16 +142,15 @@ namespace gr {
 	double y1 = abs(input[i]);
 	d_y1 = d_alpha*y1 + d_beta*d_y1;
 
-	double y2 = real(input[i]*input[i]);
+        double y2 = real(input[i]*input[i]);
 	d_y2 = d_alpha*y2 + d_beta*d_y2;
 
 	// online algorithm for calculating skewness
 	// See:
 	// http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
-	double d = abs(input[i]) - d_y1;
-	double d_i = d / (i+1);
-	double y3 = (d*d_i*i)*d_i*(i-1) - 3.0*d_i*d_y2;
-	d_y3 = d_alpha*y3 + d_beta*d_y3;
+        double d = abs(input[i]) - d_y1;
+        double d_i = d / (i+1);
+        d_y3 = (d*d_i*i)*d_i*(i-1) - 3.0*d_i*d_y2;
       }
       return noutput_items;
     }
@@ -136,9 +159,11 @@ namespace gr {
     mpsk_snr_est_skew::snr()
     {
       double y3 = d_y3*d_y3 / (d_y2*d_y2*d_y2);
-      double y1_2 = d_y1*d_y1;
-      double x = y1_2 - d_y2;
-      return 10.0*log10(y1_2 / (x + y3*y1_2));
+      d_signal = d_y1*d_y1;
+      double x = d_signal - d_y2;
+      d_noise = 2*(x + y3*d_signal);
+
+      return 10.0*log10(d_signal/d_noise);
     }
 
 
@@ -170,8 +195,9 @@ namespace gr {
     mpsk_snr_est_m2m4::snr()
     {
       double y1_2 = d_y1*d_y1;
-      return 10.0*log10(2.0*sqrt(2*y1_2 - d_y2) / 
-			(d_y1 - sqrt(2*y1_2 - d_y2)));
+      d_signal = sqrt(2*y1_2 - d_y2);
+      d_noise = d_y1 - sqrt(2*y1_2 - d_y2);
+      return 10.0*log10(d_signal / d_noise);
     }
 
 
@@ -206,12 +232,12 @@ namespace gr {
     {
       double M2 = d_y1;
       double M4 = d_y2;
-      double s = M2*(d_kw - 2) + 
+      d_signal = M2*(d_kw - 2) +
 	sqrt((4.0-d_ka*d_kw)*M2*M2 + M4*(d_ka+d_kw-4.0)) /
 	(d_ka + d_kw - 4.0);
-      double n = M2 - s;
-  
-      return 10.0*log10(s / n);
+      d_noise = M2 - d_signal;
+
+      return 10.0*log10(d_signal / d_noise);
     }
 
 
@@ -234,7 +260,7 @@ namespace gr {
 	double x1 = abs(input[i]);
 	double y1 = (x*x)*(x1*x1);
 	d_y1 = d_alpha*y1 + d_beta*d_y1;
-    
+
 	double y2 = x*x*x*x;
 	d_y2 = d_alpha*y2 + d_beta*d_y2;
       }
@@ -245,7 +271,7 @@ namespace gr {
     mpsk_snr_est_svr::snr()
     {
       double x = d_y1 / (d_y2 - d_y1);
-      return 10.0*log10(2.*((x-1) + sqrt(x*(x-1))));
+      return 10.0*log10(x - 1 + sqrt(x*(x-1)));
     }
 
   } /* namespace digital */
diff --git a/gr-digital/python/digital/qa_mpsk_snr_est.py b/gr-digital/python/digital/qa_mpsk_snr_est.py
index 032edf1c73..97d31c7686 100755
--- a/gr-digital/python/digital/qa_mpsk_snr_est.py
+++ b/gr-digital/python/digital/qa_mpsk_snr_est.py
@@ -1,24 +1,24 @@
 #!/usr/bin/env python
 #
 # Copyright 2011-2013 Free Software Foundation, Inc.
-# 
+#
 # This file is part of GNU Radio
-# 
+#
 # GNU Radio is free software; you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
 # the Free Software Foundation; either version 3, or (at your option)
 # any later version.
-# 
+#
 # GNU Radio is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
-# 
+#
 # You should have received a copy of the GNU General Public License
 # along with GNU Radio; see the file COPYING.  If not, write to
 # the Free Software Foundation, Inc., 51 Franklin Street,
 # Boston, MA 02110-1301, USA.
-# 
+#
 import random
 
 from gnuradio import gr, gr_unittest, digital, blocks
@@ -44,7 +44,7 @@ class test_mpsk_snr_est(gr_unittest.TestCase):
         result = []
         for i in xrange(1,6):
             src_data = [b+(i*n) for b,n in zip(self._bits, self._noise)]
-            
+
             src = blocks.vector_source_c(src_data)
             dst = blocks.null_sink(gr.sizeof_gr_complex)
 
@@ -55,9 +55,9 @@ class test_mpsk_snr_est(gr_unittest.TestCase):
 
             result.append(op.snr())
         return result
-            
+
     def test_mpsk_snr_est_simple(self):
-	expected_result = [11.48, 5.91, 3.30, 2.08, 1.46]
+	expected_result = [8.20, 4.99, 3.23, 2.01, 1.03]
 
         N = 10000
         alpha = 0.001
@@ -67,7 +67,7 @@ class test_mpsk_snr_est(gr_unittest.TestCase):
         self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 2)
 
     def test_mpsk_snr_est_skew(self):
-	expected_result = [11.48, 5.91, 3.30, 2.08, 1.46]
+	expected_result = [8.31, 1.83, -1.68, -3.56, -4.68]
 
         N = 10000
         alpha = 0.001
@@ -77,7 +77,7 @@ class test_mpsk_snr_est(gr_unittest.TestCase):
         self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 2)
 
     def test_mpsk_snr_est_m2m4(self):
-	expected_result = [11.02, 6.20, 4.98, 5.16, 5.66]
+	expected_result = [8.01, 3.19, 1.97, 2.15, 2.65]
 
         N = 10000
         alpha = 0.001
@@ -87,7 +87,7 @@ class test_mpsk_snr_est(gr_unittest.TestCase):
         self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 2)
 
     def test_mpsk_snr_est_svn(self):
-	expected_result = [10.92, 6.02, 4.78, 4.98, 5.51]
+	expected_result = [7.91, 3.01, 1.77, 1.97, 2.49]
 
         N = 10000
         alpha = 0.001
@@ -97,12 +97,12 @@ class test_mpsk_snr_est(gr_unittest.TestCase):
         self.assertFloatTuplesAlmostEqual(expected_result, actual_result, 2)
 
     def test_probe_mpsk_snr_est_m2m4(self):
-	expected_result = [11.02, 6.20, 4.98, 5.16, 5.66]
+	expected_result = [8.01, 3.19, 1.97, 2.15, 2.65]
 
         actual_result = []
         for i in xrange(1,6):
             src_data = [b+(i*n) for b,n in zip(self._bits, self._noise)]
-            
+
             src = blocks.vector_source_c(src_data)
 
             N = 10000
@@ -121,4 +121,3 @@ if __name__ == '__main__':
     # noise source, so these estimates have no real meaning;
     # just a sanity check.
     gr_unittest.run(test_mpsk_snr_est, "test_mpsk_snr_est.xml")
-