[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Commit-gnuradio] [gnuradio] 01/01: digital: fixes a bug in the mpsk SNR
|
From: |
git |
|
Subject: |
[Commit-gnuradio] [gnuradio] 01/01: digital: fixes a bug in the mpsk SNR estimators where they were reporting 3 dB too high. |
|
Date: |
Fri, 29 Aug 2014 20:50:22 +0000 (UTC) |
This is an automated email from the git hooks/post-receive script.
trondeau pushed a commit to branch maint
in repository gnuradio.
commit d56ad00edb24060d6cfe8d0f31f2dbd456ca7f8a
Author: Tom Rondeau <address@hidden>
Date: Fri Aug 29 15:51:26 2014 -0400
digital: fixes a bug in the mpsk SNR estimators where they were reporting 3
dB too high.
---
gr-digital/examples/snr_estimators.py | 71 +++++++++---------
gr-digital/include/gnuradio/digital/mpsk_snr_est.h | 31 +++++---
gr-digital/lib/mpsk_snr_est.cc | 84 ++++++++++++++--------
gr-digital/python/digital/qa_mpsk_snr_est.py | 27 ++++---
4 files changed, 127 insertions(+), 86 deletions(-)
diff --git a/gr-digital/examples/snr_estimators.py
b/gr-digital/examples/snr_estimators.py
index 9eae786..31efe6b 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 46df061..5398745 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 6edb027..4e7825b 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 032edf1..97d31c7 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")
-