summaryrefslogtreecommitdiff
path: root/gr-digital/examples
diff options
context:
space:
mode:
Diffstat (limited to 'gr-digital/examples')
-rw-r--r--gr-digital/examples/CMakeLists.txt1
-rw-r--r--gr-digital/examples/Makefile.am3
-rwxr-xr-xgr-digital/examples/narrowband/digital_bert_rx.py2
-rwxr-xr-xgr-digital/examples/narrowband/digital_bert_tx.py7
-rwxr-xr-xgr-digital/examples/snr_estimators.py174
5 files changed, 184 insertions, 3 deletions
diff --git a/gr-digital/examples/CMakeLists.txt b/gr-digital/examples/CMakeLists.txt
index 2645557cc..7b94f745c 100644
--- a/gr-digital/examples/CMakeLists.txt
+++ b/gr-digital/examples/CMakeLists.txt
@@ -25,6 +25,7 @@ GR_PYTHON_INSTALL(PROGRAMS
example_timing.py
run_length.py
gen_whitener.py
+ snr_estimators.py
DESTINATION ${GR_PKG_DATA_DIR}/examples/digital
COMPONENT "digital_python"
)
diff --git a/gr-digital/examples/Makefile.am b/gr-digital/examples/Makefile.am
index 3841bbd0e..3e6b32e85 100644
--- a/gr-digital/examples/Makefile.am
+++ b/gr-digital/examples/Makefile.am
@@ -29,7 +29,8 @@ dist_basicdata_SCRIPTS = \
example_fll.py \
example_timing.py \
run_length.py \
- gen_whitener.py
+ gen_whitener.py \
+ snr_estimators.py
nbdatadir = $(dig_examples_prefix)/narrowband
dist_nbdata_SCRIPTS = \
diff --git a/gr-digital/examples/narrowband/digital_bert_rx.py b/gr-digital/examples/narrowband/digital_bert_rx.py
index 28331310d..dfed0eee7 100755
--- a/gr-digital/examples/narrowband/digital_bert_rx.py
+++ b/gr-digital/examples/narrowband/digital_bert_rx.py
@@ -113,9 +113,11 @@ class rx_psk_block(gr.top_block):
self._demodulator = self._demodulator_class(**demod_kwargs)
if(options.rx_freq is not None):
+ symbol_rate = options.bitrate / self._demodulator.bits_per_symbol()
self._source = uhd_receiver(options.args, options.bitrate,
options.samples_per_symbol,
options.rx_freq, options.rx_gain,
+ options.spec,
options.antenna, options.verbose)
options.samples_per_symbol = self._source._sps
diff --git a/gr-digital/examples/narrowband/digital_bert_tx.py b/gr-digital/examples/narrowband/digital_bert_tx.py
index 46f4f9097..f29e997af 100755
--- a/gr-digital/examples/narrowband/digital_bert_tx.py
+++ b/gr-digital/examples/narrowband/digital_bert_tx.py
@@ -67,9 +67,11 @@ class tx_psk_block(gr.top_block):
self._modulator = self._modulator_class(**mod_kwargs)
if(options.tx_freq is not None):
- self._sink = uhd_transmitter(options.args, options.bitrate,
+ symbol_rate = options.bitrate / self._modulator.bits_per_symbol()
+ self._sink = uhd_transmitter(options.args, symbol_rate,
options.samples_per_symbol,
options.tx_freq, options.tx_gain,
+ options.spec,
options.antenna, options.verbose)
options.samples_per_symbol = self._sink._sps
@@ -87,7 +89,8 @@ class tx_psk_block(gr.top_block):
verbose=options.verbose,
log=options.log)
- self.connect(self._transmitter, self._sink)
+ self.amp = gr.multiply_const_cc(options.amplitude)
+ self.connect(self._transmitter, self.amp, self._sink)
def get_options(mods):
diff --git a/gr-digital/examples/snr_estimators.py b/gr-digital/examples/snr_estimators.py
new file mode 100755
index 000000000..432abd455
--- /dev/null
+++ b/gr-digital/examples/snr_estimators.py
@@ -0,0 +1,174 @@
+#!/usr/bin/env python
+
+import sys
+
+try:
+ import scipy
+ from scipy import stats
+except ImportError:
+ print "Error: Program requires scipy (www.scipy.org)."
+ sys.exit(1)
+
+try:
+ import pylab
+except ImportError:
+ print "Error: Program requires Matplotlib (matplotlib.sourceforge.net)."
+ sys.exit(1)
+
+from gnuradio import gr, digital
+from optparse import OptionParser
+from gnuradio.eng_option import eng_option
+
+'''
+This example program uses Python and GNU Radio to calculate SNR of a
+noise BPSK signal to compare them.
+
+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):
+ 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)
+ mean = mean + delta_n
+ 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
+
+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
+ 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)
+
+ skw = y4*y4 / (y2*y2*y2);
+ snr_rat = y1*y1 / (y3 + skw*y1*y1)
+ 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))
+ return 10.0*scipy.log10(snr_rat), snr_rat
+
+def snr_est_svr(signal):
+ N = len(signal)
+ ssum = 0
+ msum = 0
+ for i in xrange(1, N):
+ ssum += (abs(signal[i])**2)*(abs(signal[i-1])**2)
+ msum += (abs(signal[i])**4)
+ 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)))
+ return 10.0*scipy.log10(snr_rat), snr_rat
+
+
+def main():
+ gr_estimators = {"simple": digital.SNR_EST_SIMPLE,
+ "skew": digital.SNR_EST_SKEW,
+ "m2m4": digital.SNR_EST_M2M4,
+ "svr": digital.SNR_EST_SVR}
+ py_estimators = {"simple": snr_est_simple,
+ "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]")
+ parser.add_option("", "--snr-min", type="float", default=-5,
+ help="Minimum SNR [default=%default]")
+ parser.add_option("", "--snr-max", type="float", default=20,
+ help="Maximum SNR [default=%default]")
+ parser.add_option("", "--snr-step", type="float", default=0.5,
+ help="SNR step amount [default=%default]")
+ parser.add_option("-t", "--type", type="choice",
+ choices=gr_estimators.keys(), default="simple",
+ help="Estimator type {0} [default=%default]".format(
+ gr_estimators.keys()))
+ (options, args) = parser.parse_args ()
+
+ N = options.nsamples
+ xx = scipy.random.randn(N)
+ xy = scipy.random.randn(N)
+ bits = 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
+
+ n_cpx = xx + 1j*xy
+
+ py_est = py_estimators[options.type]
+ gr_est = gr_estimators[options.type]
+
+ SNR_min = options.snr_min
+ SNR_max = options.snr_max
+ SNR_step = options.snr_step
+ 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)
+ yy = bits + n_cpx/scale
+ print "SNR: ", snr
+
+ Sknown = scipy.mean(yy**2)
+ Nknown = scipy.var(n_cpx/scale)/2
+ snr0 = Sknown/Nknown
+ snr0dB = 10.0*scipy.log10(snr0)
+ snr_known.append(snr0dB)
+
+ snrdB, snr = py_est(yy)
+ snr_python.append(snrdB)
+
+ gr_src = gr.vector_source_c(bits.tolist(), False)
+ gr_snr = digital.mpsk_snr_est_cc(gr_est, ntag, 0.001)
+ gr_chn = gr.channel_model(1.0/scale)
+ gr_snk = gr.null_sink(gr.sizeof_gr_complex)
+ tb = gr.top_block()
+ tb.connect(gr_src, gr_chn, gr_snr, gr_snk)
+ tb.run()
+
+ snr_gr.append(gr_snr.snr())
+
+ f1 = pylab.figure(1)
+ s1 = f1.add_subplot(1,1,1)
+ s1.plot(SNR_dB, snr_known, "k-o", linewidth=2, label="Known")
+ s1.plot(SNR_dB, snr_python, "b-o", linewidth=2, label="Python")
+ s1.plot(SNR_dB, snr_gr, "g-o", linewidth=2, label="GNU Radio")
+ s1.grid(True)
+ s1.set_title('SNR Estimators')
+ s1.set_xlabel('SNR (dB)')
+ s1.set_ylabel('Estimated SNR')
+ s1.legend()
+
+ pylab.show()
+
+
+if __name__ == "__main__":
+ main()
+