diff options
Diffstat (limited to 'gr-digital')
25 files changed, 1547 insertions, 5 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() + diff --git a/gr-digital/grc/digital_dxpsk_demod.xml b/gr-digital/grc/digital_dxpsk_demod.xml index d5e742097..cfd474f68 100644 --- a/gr-digital/grc/digital_dxpsk_demod.xml +++ b/gr-digital/grc/digital_dxpsk_demod.xml @@ -33,6 +33,7 @@ <make>digital.$(type)_demod( samples_per_symbol=$samples_per_symbol, excess_bw=$excess_bw, + freq_bw=$freq_bw, phase_bw=$phase_bw, timing_bw=$timing_bw, gray_coded=$gray_coded, @@ -67,6 +68,12 @@ <type>real</type> </param> <param> + <name>FLL Bandwidth</name> + <key>freq_bw</key> + <value>6.28/100.0</value> + <type>real</type> + </param> + <param> <name>Phase Loop Bandwidth</name> <key>phase_bw</key> <value>6.28/100.0</value> diff --git a/gr-digital/include/CMakeLists.txt b/gr-digital/include/CMakeLists.txt index cf20bd1e7..81ed8d368 100644 --- a/gr-digital/include/CMakeLists.txt +++ b/gr-digital/include/CMakeLists.txt @@ -22,6 +22,7 @@ ######################################################################## install(FILES digital_api.h + digital_impl_mpsk_snr_est.h digital_binary_slicer_fb.h digital_clock_recovery_mm_cc.h digital_clock_recovery_mm_ff.h @@ -37,12 +38,14 @@ install(FILES digital_kurtotic_equalizer_cc.h digital_metric_type.h digital_mpsk_receiver_cc.h + digital_mpsk_snr_est_cc.h digital_ofdm_cyclic_prefixer.h digital_ofdm_frame_acquisition.h digital_ofdm_frame_sink.h digital_ofdm_insert_preamble.h digital_ofdm_mapper_bcv.h digital_ofdm_sampler.h + digital_probe_mpsk_snr_est_c.h digital_gmskmod_bc.h digital_cpmmod_bc.h DESTINATION ${GR_INCLUDE_DIR}/gnuradio diff --git a/gr-digital/include/Makefile.am b/gr-digital/include/Makefile.am index 8ce3a94e8..3cf186d3d 100644 --- a/gr-digital/include/Makefile.am +++ b/gr-digital/include/Makefile.am @@ -24,6 +24,7 @@ include $(top_srcdir)/Makefile.common # These headers get installed in ${prefix}/include/gnuradio grinclude_HEADERS = \ digital_api.h \ + digital_impl_mpsk_snr_est.h \ digital_binary_slicer_fb.h \ digital_clock_recovery_mm_cc.h \ digital_clock_recovery_mm_ff.h \ @@ -39,12 +40,14 @@ grinclude_HEADERS = \ digital_kurtotic_equalizer_cc.h \ digital_metric_type.h \ digital_mpsk_receiver_cc.h \ + digital_mpsk_snr_est_cc.h \ digital_ofdm_cyclic_prefixer.h \ digital_ofdm_frame_acquisition.h \ digital_ofdm_frame_sink.h \ digital_ofdm_insert_preamble.h \ digital_ofdm_mapper_bcv.h \ digital_ofdm_sampler.h \ + digital_probe_mpsk_snr_est_c.h \ digital_gmskmod_bc.h \ digital_cpmmod_bc.h diff --git a/gr-digital/include/digital_impl_mpsk_snr_est.h b/gr-digital/include/digital_impl_mpsk_snr_est.h new file mode 100644 index 000000000..df7dbadec --- /dev/null +++ b/gr-digital/include/digital_impl_mpsk_snr_est.h @@ -0,0 +1,279 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ +#ifndef INCLUDED_DIGITAL_IMPL_MPSK_SNR_EST_H +#define INCLUDED_DIGITAL_IMPL_MPSK_SNR_EST_H + +#include <digital_api.h> +#include <gr_sync_block.h> + +//! Enum for the type of SNR estimator to select +/*! \ingroup snr_blk + * \anchor ref_snr_est_types + * + * Below are some ROUGH estimates of what values of SNR each of these + * types of estimators is good for. In general, these offer a + * trade-off between accuracy and performance. + * + * \li SNR_EST_SIMPLE: Simple estimator (>= 7 dB) + * \li SNR_EST_SKEW: Skewness-base est (>= 5 dB) + * \li SNR_EST_M2M4: 2nd & 4th moment est (>= 1 dB) + * \li SNR_EST_SVR: SVR-based est (>= 0dB) +*/ +enum snr_est_type_t { + SNR_EST_SIMPLE = 0, // Simple estimator (>= 7 dB) + SNR_EST_SKEW, // Skewness-base est (>= 5 dB) + SNR_EST_M2M4, // 2nd & 4th moment est (>= 1 dB) + SNR_EST_SVR // SVR-based est (>= 0dB) +}; + +/*! \brief A parent class for SNR estimators, specifically for M-PSK + * signals in AWGN channels. + * \ingroup snr_blk + */ +class DIGITAL_API digital_impl_mpsk_snr_est +{ + protected: + double d_alpha, d_beta; + + public: + /*! Constructor + * + * Parameters: + * \param alpha: the update rate of internal running average + * calculations. + */ + digital_impl_mpsk_snr_est(double alpha); + virtual ~digital_impl_mpsk_snr_est(); + + //! Get the running-average coefficient + double alpha() const; + + //! Set the running-average coefficient + void set_alpha(double alpha); + + //! Update the current registers + virtual int update(int noutput_items, + const gr_complex *in); + + //! Use the register values to compute a new estimate + virtual double snr(); +}; + + +//! \brief SNR Estimator using simple mean/variance estimates. +/*! \ingroup snr_blk + * + * A very simple SNR estimator that just uses mean and variance + * estimates of an M-PSK constellation. This esimator is quick and + * cheap and accurate for high SNR (above 7 dB or so) but quickly + * starts to overestimate the SNR at low SNR. + */ +class DIGITAL_API digital_impl_mpsk_snr_est_simple : + public digital_impl_mpsk_snr_est +{ + private: + double d_y1, d_y2; + + public: + /*! Constructor + * + * Parameters: + * \param alpha: the update rate of internal running average + * calculations. + */ + digital_impl_mpsk_snr_est_simple(double alpha); + ~digital_impl_mpsk_snr_est_simple() {} + + int update(int noutput_items, + const gr_complex *in); + double snr(); +}; + + +//! \brief SNR Estimator using skewness correction. +/*! \ingroup snr_blk + * + * This is an estimator that came from a discussion between Tom + * Rondeau and fred harris with no known paper reference. The idea is + * that at low SNR, the variance estimations will be 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. + */ +class DIGITAL_API digital_impl_mpsk_snr_est_skew : + public digital_impl_mpsk_snr_est +{ + private: + double d_y1, d_y2, d_y3; + + public: + /*! Constructor + * + * Parameters: + * \param alpha: the update rate of internal running average + * calculations. + */ + digital_impl_mpsk_snr_est_skew(double alpha); + ~digital_impl_mpsk_snr_est_skew() {} + + int update(int noutput_items, + const gr_complex *in); + double snr(); +}; + + +//! \brief SNR Estimator using 2nd and 4th-order moments. +/*! \ingroup snr_blk + * + * An SNR estimator for M-PSK signals that uses 2nd (M2) and 4th (M4) + * order moments. This estimator uses knowledge of the kurtosis of + * the signal (k_a) and noise (k_w) to make its estimation. We use + * Beaulieu's approximations here to M-PSK signals and AWGN channels + * such that k_a=1 and k_w=2. These approximations significantly + * reduce the complexity of the calculations (and computations) + * required. + * + * Reference: + * D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR + * estimation techniques for the AWGN channel," IEEE + * Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000. + */ +class DIGITAL_API digital_impl_mpsk_snr_est_m2m4 : + public digital_impl_mpsk_snr_est +{ + private: + double d_y1, d_y2; + + public: + /*! Constructor + * + * Parameters: + * \param alpha: the update rate of internal running average + * calculations. + */ + digital_impl_mpsk_snr_est_m2m4(double alpha); + ~digital_impl_mpsk_snr_est_m2m4() {} + + int update(int noutput_items, + const gr_complex *in); + double snr(); +}; + + +//! \brief SNR Estimator using 2nd and 4th-order moments. +/*! \ingroup snr_blk + * + * An SNR estimator for M-PSK signals that uses 2nd (M2) and 4th (M4) + * order moments. This estimator uses knowledge of the kurtosis of + * the signal (k_a) and noise (k_w) to make its estimation. In this + * case, you can set your own estimations for k_a and k_w, the + * kurtosis of the signal and noise, to fit this estimation better to + * your signal and channel conditions. + * + * A word of warning: this estimator has not been fully tested or + * proved with any amount of rigor. The estimation for M4 in + * particular might be ignoring effectf of when k_a and k_w are + * different. Use this estimator with caution and a copy of the + * reference on hand. + * + * The digital_mpsk_snr_est_m2m4 assumes k_a and k_w to simplify the + * computations for M-PSK and AWGN channels. Use that estimator + * unless you have a way to guess or estimate these values here. + * + * 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. + * + * Reference used in derivation: + * D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR + * estimation techniques for the AWGN channel," IEEE + * Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000. + */ +class DIGITAL_API digital_impl_snr_est_m2m4 : + public digital_impl_mpsk_snr_est +{ + private: + double d_y1, d_y2; + double d_ka, d_kw; + + public: + /*! Constructor + * + * Parameters: + * \param alpha: the update rate of internal running average + * calculations. + * \param ka: estimate of the signal kurtosis (1 for PSK) + * \param kw: estimate of the channel noise kurtosis (2 for AWGN) + */ + digital_impl_snr_est_m2m4(double alpha, double ka, double kw); + ~digital_impl_snr_est_m2m4() {} + + int update(int noutput_items, + const gr_complex *in); + double snr(); +}; + + +//! \brief Signal-to-Variation Ratio SNR Estimator. +/*! \ingroup snr_blk + * + * This estimator actually comes from an SNR estimator for M-PSK + * signals in fading channels, but this implementation is + * specifically for AWGN channels. The math was simplified to assume + * a signal and noise kurtosis (k_a and k_w) for M-PSK signals in + * AWGN. These approximations significantly reduce the complexity of + * the calculations (and computations) required. + * + * Original paper: + * A. L. Brandao, L. B. Lopes, and D. C. McLernon, "In-service + * monitoring of multipath delay and cochannel interference for + * indoor mobile communication systems," Proc. IEEE + * Int. Conf. Communications, vol. 3, pp. 1458-1462, May 1994. + * + * Reference: + * D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR + * estimation techniques for the AWGN channel," IEEE + * Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000. + */ +class DIGITAL_API digital_impl_mpsk_snr_est_svr : + public digital_impl_mpsk_snr_est +{ + private: + double d_y1, d_y2; + + public: + /*! Constructor + * + * Parameters: + * \param alpha: the update rate of internal running average + * calculations. + */ + digital_impl_mpsk_snr_est_svr(double alpha); + ~digital_impl_mpsk_snr_est_svr() {} + + int update(int noutput_items, + const gr_complex *in); + double snr(); +}; + +#endif /* INCLUDED_DIGITAL_IMPL_MPSK_SNR_EST_H */ diff --git a/gr-digital/include/digital_mpsk_snr_est_cc.h b/gr-digital/include/digital_mpsk_snr_est_cc.h new file mode 100644 index 000000000..2cbd98bab --- /dev/null +++ b/gr-digital/include/digital_mpsk_snr_est_cc.h @@ -0,0 +1,115 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ +#ifndef INCLUDED_DIGITAL_MPSK_SNR_EST_CC_H +#define INCLUDED_DIGITAL_MPSK_SNR_EST_CC_H + +#include <digital_api.h> +#include <gr_sync_block.h> +#include <digital_impl_mpsk_snr_est.h> + +class digital_mpsk_snr_est_cc; +typedef boost::shared_ptr<digital_mpsk_snr_est_cc> digital_mpsk_snr_est_cc_sptr; + +DIGITAL_API digital_mpsk_snr_est_cc_sptr +digital_make_mpsk_snr_est_cc(snr_est_type_t type, + int tag_nsamples=10000, + double alpha=0.001); + +//! \brief A block for computing SNR of a signal. +/*! \ingroup snr_blk + * + * This block can be used to monitor and retrieve estimations of the + * signal SNR. It is designed to work in a flowgraph and passes all + * incoming data along to its output. + * + * The block is designed for use with M-PSK signals especially. The + * type of estimator is specified as the \p type parameter in the + * constructor. The estimators tend to trade off performance for + * accuracy, although experimentation should be done to figure out + * the right approach for a given implementation. Further, the + * current set of estimators are designed and proven theoretically + * under AWGN conditions; some amount of error should be assumed + * and/or estimated for real channel conditions. + */ +class DIGITAL_API digital_mpsk_snr_est_cc : public gr_sync_block +{ + private: + snr_est_type_t d_type; + int d_nsamples, d_count; + double d_alpha; + digital_impl_mpsk_snr_est *d_snr_est; + + //d_key is the tag name, 'snr', d_me is the block name + unique ID + pmt::pmt_t d_key, d_me; + + /*! Factory function returning shared pointer of this class + * + * Parameters: + * + * \param type: the type of estimator to use \ref ref_snr_est_types + * "snr_est_type_t" for details about the available types. + * \param tag_nsamples: after this many samples, a tag containing + * the SNR (key='snr') will be sent + * \param alpha: the update rate of internal running average + * calculations. + */ + friend DIGITAL_API digital_mpsk_snr_est_cc_sptr + digital_make_mpsk_snr_est_cc(snr_est_type_t type, + int tag_nsamples, + double alpha); + + // Private constructor + digital_mpsk_snr_est_cc(snr_est_type_t type, + int tag_nsamples, + double alpha); + +public: + + ~digital_mpsk_snr_est_cc(); + + int work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + + //! Return the estimated signal-to-noise ratio in decibels + double snr(); + + //! Return the type of estimator in use + snr_est_type_t type() const; + + //! Return how many samples between SNR tags + int tag_nsample() const; + + //! Get the running-average coefficient + double alpha() const; + + //! Set type of estimator to use + void set_type(snr_est_type_t t); + + //! Set the number of samples between SNR tags + void set_tag_nsample(int n); + + //! Set the running-average coefficient + void set_alpha(double alpha); +}; + +#endif /* INCLUDED_DIGITAL_MPSK_SNR_EST_CC_H */ diff --git a/gr-digital/include/digital_probe_mpsk_snr_est_c.h b/gr-digital/include/digital_probe_mpsk_snr_est_c.h new file mode 100644 index 000000000..a78e90412 --- /dev/null +++ b/gr-digital/include/digital_probe_mpsk_snr_est_c.h @@ -0,0 +1,113 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ +#ifndef INCLUDED_DIGITAL_PROBE_MPSK_SNR_EST_C_H +#define INCLUDED_DIGITAL_PROBE_MPSK_SNR_EST_C_H + +#include <digital_api.h> +#include <gr_sync_block.h> +#include <digital_impl_mpsk_snr_est.h> + +class digital_probe_mpsk_snr_est_c; +typedef boost::shared_ptr<digital_probe_mpsk_snr_est_c> digital_probe_mpsk_snr_est_c_sptr; + +DIGITAL_API digital_probe_mpsk_snr_est_c_sptr +digital_make_probe_mpsk_snr_est_c(snr_est_type_t type, + int msg_nsamples=10000, + double alpha=0.001); + +//! \brief A probe for computing SNR of a signal. +/*! \ingroup snr_blk + * + * This is a probe block (a sink) that can be used to monitor and + * retrieve estimations of the signal SNR. This probe is designed for + * use with M-PSK signals especially. The type of estimator is + * specified as the \p type parameter in the constructor. The + * estimators tend to trade off performance for accuracy, although + * experimentation should be done to figure out the right approach + * for a given implementation. Further, the current set of estimators + * are designed and proven theoretically under AWGN conditions; some + * amount of error should be assumed and/or estimated for real + * channel conditions. + */ +class DIGITAL_API digital_probe_mpsk_snr_est_c : public gr_sync_block +{ + private: + snr_est_type_t d_type; + int d_nsamples, d_count; + double d_alpha; + digital_impl_mpsk_snr_est *d_snr_est; + + //d_key is the message name, 'snr' + pmt::pmt_t d_key; + + /*! Factory function returning shared pointer of this class + * + * Parameters: + * + * \param type: the type of estimator to use \ref ref_snr_est_types + * "snr_est_type_t" for details about the available types. + * \param msg_nsamples: [not implemented yet] after this many + * samples, a message containing the SNR (key='snr') will be sent + * \param alpha: the update rate of internal running average + * calculations. + */ + friend DIGITAL_API digital_probe_mpsk_snr_est_c_sptr + digital_make_probe_mpsk_snr_est_c(snr_est_type_t type, + int msg_nsamples, + double alpha); + + //! Private constructor + digital_probe_mpsk_snr_est_c(snr_est_type_t type, + int msg_nsamples, + double alpha); + +public: + + ~digital_probe_mpsk_snr_est_c(); + + int work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); + + //! Return the estimated signal-to-noise ratio in decibels + double snr(); + + //! Return the type of estimator in use + snr_est_type_t type() const; + + //! Return how many samples between SNR messages + int msg_nsample() const; + + //! Get the running-average coefficient + double alpha() const; + + //! Set type of estimator to use + void set_type(snr_est_type_t t); + + //! Set the number of samples between SNR messages + void set_msg_nsample(int n); + + //! Set the running-average coefficient + void set_alpha(double alpha); +}; + +#endif /* INCLUDED_DIGITAL_PROBE_MPSK_SNR_EST_C_H */ diff --git a/gr-digital/lib/CMakeLists.txt b/gr-digital/lib/CMakeLists.txt index b90757111..779972ff3 100644 --- a/gr-digital/lib/CMakeLists.txt +++ b/gr-digital/lib/CMakeLists.txt @@ -32,6 +32,7 @@ link_directories(${Boost_LIBRARY_DIRS}) # Setup library ######################################################################## list(APPEND gr_digital_sources + digital_impl_mpsk_snr_est.cc digital_binary_slicer_fb.cc digital_clock_recovery_mm_cc.cc digital_clock_recovery_mm_ff.cc @@ -46,12 +47,14 @@ list(APPEND gr_digital_sources digital_lms_dd_equalizer_cc.cc digital_kurtotic_equalizer_cc.cc digital_mpsk_receiver_cc.cc + digital_mpsk_snr_est_cc.cc digital_ofdm_cyclic_prefixer.cc digital_ofdm_frame_acquisition.cc digital_ofdm_frame_sink.cc digital_ofdm_insert_preamble.cc digital_ofdm_mapper_bcv.cc digital_ofdm_sampler.cc + digital_probe_mpsk_snr_est_c.cc digital_gmskmod_bc.cc digital_cpmmod_bc.cc ) diff --git a/gr-digital/lib/Makefile.am b/gr-digital/lib/Makefile.am index 2860974ca..d5ad199e3 100644 --- a/gr-digital/lib/Makefile.am +++ b/gr-digital/lib/Makefile.am @@ -27,6 +27,7 @@ AM_CPPFLAGS = $(STD_DEFINES_AND_INCLUDES) $(PYTHON_CPPFLAGS) \ lib_LTLIBRARIES = libgnuradio-digital.la libgnuradio_digital_la_SOURCES = \ + digital_impl_mpsk_snr_est.cc \ digital_binary_slicer_fb.cc \ digital_clock_recovery_mm_cc.cc \ digital_clock_recovery_mm_ff.cc \ @@ -41,12 +42,14 @@ libgnuradio_digital_la_SOURCES = \ digital_lms_dd_equalizer_cc.cc \ digital_kurtotic_equalizer_cc.cc \ digital_mpsk_receiver_cc.cc \ + digital_mpsk_snr_est_cc.cc \ digital_ofdm_cyclic_prefixer.cc \ digital_ofdm_frame_acquisition.cc \ digital_ofdm_frame_sink.cc \ digital_ofdm_insert_preamble.cc \ digital_ofdm_mapper_bcv.cc \ digital_ofdm_sampler.cc \ + digital_probe_mpsk_snr_est_c.cc \ digital_gmskmod_bc.cc \ digital_cpmmod_bc.cc diff --git a/gr-digital/lib/digital_constellation.cc b/gr-digital/lib/digital_constellation.cc index 0c100f38e..d9a53c493 100644 --- a/gr-digital/lib/digital_constellation.cc +++ b/gr-digital/lib/digital_constellation.cc @@ -357,10 +357,8 @@ digital_constellation_psk::get_sector (const gr_complex *sample) float phase = arg(*sample); float width = M_TWOPI / n_sectors; int sector = floor(phase/width + 0.5); - unsigned int u_sector; if (sector < 0) sector += n_sectors; - u_sector = sector; return sector; } diff --git a/gr-digital/lib/digital_impl_mpsk_snr_est.cc b/gr-digital/lib/digital_impl_mpsk_snr_est.cc new file mode 100644 index 000000000..38177083f --- /dev/null +++ b/gr-digital/lib/digital_impl_mpsk_snr_est.cc @@ -0,0 +1,256 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <digital_impl_mpsk_snr_est.h> +#include <cstdio> + +digital_impl_mpsk_snr_est::digital_impl_mpsk_snr_est(double alpha) +{ + set_alpha(alpha); +} + +digital_impl_mpsk_snr_est::~digital_impl_mpsk_snr_est() +{} + +void +digital_impl_mpsk_snr_est::set_alpha(double alpha) +{ + d_alpha = alpha; + d_beta = 1.0-alpha; +} + +double +digital_impl_mpsk_snr_est::alpha() const +{ + return d_alpha; +} + +int +digital_impl_mpsk_snr_est::update(int noutput_items, + const gr_complex *in) +{ + throw std::runtime_error("digital_impl_mpsk_snr_est: Unimplemented"); +} + +double +digital_impl_mpsk_snr_est::snr() +{ + throw std::runtime_error("digital_impl_mpsk_snr_est: Unimplemented"); +} + + +/********************************************************************/ + + +digital_impl_mpsk_snr_est_simple::digital_impl_mpsk_snr_est_simple( + double alpha) : + digital_impl_mpsk_snr_est(alpha) +{ + d_y1 = 0; + d_y2 = 0; +} + +int +digital_impl_mpsk_snr_est_simple::update( + int noutput_items, + const gr_complex *in) +{ + for (int i = 0; i < noutput_items; i++){ + double y1 = abs(in[i]); + d_y1 = d_alpha*y1 + d_beta*d_y1; + + double y2 = real(in[i]*in[i]); + d_y2 = d_alpha*y2 + d_beta*d_y2; + } + return noutput_items; +} + +double +digital_impl_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); +} + + +/********************************************************************/ + + +digital_impl_mpsk_snr_est_skew::digital_impl_mpsk_snr_est_skew( + double alpha) : + digital_impl_mpsk_snr_est(alpha) +{ + d_y1 = 0; + d_y2 = 0; + d_y3 = 0; +} + + +int +digital_impl_mpsk_snr_est_skew::update( + int noutput_items, + const gr_complex *in) +{ + for (int i = 0; i < noutput_items; i++){ + double y1 = abs(in[i]); + d_y1 = d_alpha*y1 + d_beta*d_y1; + + double y2 = real(in[i]*in[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(in[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; + } + return noutput_items; +} + +double +digital_impl_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)); +} + + +/********************************************************************/ + + +digital_impl_mpsk_snr_est_m2m4::digital_impl_mpsk_snr_est_m2m4( + double alpha) : + digital_impl_mpsk_snr_est(alpha) +{ + d_y1 = 0; + d_y2 = 0; +} + +int +digital_impl_mpsk_snr_est_m2m4::update( + int noutput_items, + const gr_complex *in) +{ + for (int i = 0; i < noutput_items; i++){ + double y1 = abs(in[i])*abs(in[i]); + d_y1 = d_alpha*y1 + d_beta*d_y1; + + double y2 = abs(in[i])*abs(in[i])*abs(in[i])*abs(in[i]); + d_y2 = d_alpha*y2 + d_beta*d_y2; + } + return noutput_items; +} + +double +digital_impl_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))); +} + + +/********************************************************************/ + + +digital_impl_snr_est_m2m4::digital_impl_snr_est_m2m4( + double alpha, double ka, double kw) : + digital_impl_mpsk_snr_est(alpha) +{ + d_y1 = 0; + d_y2 = 0; + d_ka = ka; + d_kw = kw; +} + +int +digital_impl_snr_est_m2m4::update( + int noutput_items, + const gr_complex *in) +{ + for (int i = 0; i < noutput_items; i++) { + double y1 = abs(in[i])*abs(in[i]); + d_y1 = d_alpha*y1 + d_beta*d_y1; + + double y2 = abs(in[i])*abs(in[i])*abs(in[i])*abs(in[i]); + d_y2 = d_alpha*y2 + d_beta*d_y2; + } + return noutput_items; +} + +double +digital_impl_snr_est_m2m4::snr() +{ + double M2 = d_y1; + double M4 = d_y2; + double s = 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); +} + + +/********************************************************************/ + + +digital_impl_mpsk_snr_est_svr::digital_impl_mpsk_snr_est_svr( + double alpha) : + digital_impl_mpsk_snr_est(alpha) +{ + d_y1 = 0; + d_y2 = 0; +} + +int +digital_impl_mpsk_snr_est_svr::update( + int noutput_items, + const gr_complex *in) +{ + for (int i = 0; i < noutput_items; i++){ + double x = abs(in[i]); + double x1 = abs(in[i-1]); + 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; + } + return noutput_items; +} + +double +digital_impl_mpsk_snr_est_svr::snr() +{ + double x = d_y1 / (d_y2 - d_y1); + return 10.0*log10(2.*((x-1) + sqrt(x*(x-1)))); +} diff --git a/gr-digital/lib/digital_mpsk_snr_est_cc.cc b/gr-digital/lib/digital_mpsk_snr_est_cc.cc new file mode 100644 index 000000000..b5a60f0d3 --- /dev/null +++ b/gr-digital/lib/digital_mpsk_snr_est_cc.cc @@ -0,0 +1,186 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <digital_mpsk_snr_est_cc.h> +#include <gr_io_signature.h> +#include <cstdio> + +digital_mpsk_snr_est_cc_sptr +digital_make_mpsk_snr_est_cc(snr_est_type_t type, + int tag_nsamples, + double alpha) +{ + return gnuradio::get_initial_sptr(new digital_mpsk_snr_est_cc( + type, tag_nsamples, alpha)); +} + +digital_mpsk_snr_est_cc::digital_mpsk_snr_est_cc(snr_est_type_t type, + int tag_nsamples, + double alpha) + : gr_sync_block ("mpsk_snr_est_cc", + gr_make_io_signature(1, 1, sizeof(gr_complex)), + gr_make_io_signature(1, 1, sizeof(gr_complex))) +{ + d_snr_est = NULL; + + d_type = type; + d_nsamples = tag_nsamples; + d_count = 0; + set_alpha(alpha); + + set_type(type); + + // at least 1 estimator has to look back + set_history(2); + + std::stringstream str; + str << name() << unique_id(); + d_me = pmt::pmt_string_to_symbol(str.str()); + d_key = pmt::pmt_string_to_symbol("snr"); +} + +digital_mpsk_snr_est_cc::~digital_mpsk_snr_est_cc() +{ + if(d_snr_est) + delete d_snr_est; +} + +int +digital_mpsk_snr_est_cc::work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + // This is a pass-through block; copy input to output + memcpy(output_items[0], input_items[0], + noutput_items * sizeof(gr_complex)); + + const gr_complex *in = (const gr_complex*)input_items[0]; + + // Update, calculate, and issue an SNR tag every d_nsamples + int index = 0, x = 0; + int64_t nwritten = nitems_written(0); + while(index + (d_nsamples-d_count) <= noutput_items) { + x = d_nsamples - d_count; + nwritten += x; + + // Update the SNR estimate registers from the current input + d_snr_est->update(x, &in[index]); + + // Issue a tag with the SNR data + pmt::pmt_t pmt_snr = pmt::pmt_from_double(d_snr_est->snr()); + add_item_tag(0, // stream ID + nwritten, // tag's sample number + d_key, // snr key + pmt_snr, // SNR + d_me); // block src id + + index += x; + d_count = 0; + } + + // Keep track of remaining items and update estimators + x = noutput_items - index; + d_count += x; + d_snr_est->update(x, &in[index]); + + return noutput_items; +} + +double +digital_mpsk_snr_est_cc::snr() +{ + if(d_snr_est) + return d_snr_est->snr(); + else + throw std::runtime_error("digital_mpsk_snr_est_cc:: No SNR estimator defined.\n"); +} + +snr_est_type_t +digital_mpsk_snr_est_cc::type() const +{ + return d_type; +} + +int +digital_mpsk_snr_est_cc::tag_nsample() const +{ + return d_nsamples; +} + +double +digital_mpsk_snr_est_cc::alpha() const +{ + return d_alpha; +} + +void +digital_mpsk_snr_est_cc::set_type(snr_est_type_t t) +{ + d_type = t; + + if(d_snr_est) + delete d_snr_est; + + switch (d_type) { + case(SNR_EST_SIMPLE): + d_snr_est = new digital_impl_mpsk_snr_est_simple(d_alpha); + break; + case(SNR_EST_SKEW): + d_snr_est = new digital_impl_mpsk_snr_est_skew(d_alpha); + break; + case(SNR_EST_M2M4): + d_snr_est = new digital_impl_mpsk_snr_est_m2m4(d_alpha); + break; + case(SNR_EST_SVR): + d_snr_est = new digital_impl_mpsk_snr_est_svr(d_alpha); + break; + default: + throw std::invalid_argument("digital_mpsk_snr_est_cc: unknown type specified.\n"); + } +} + +void +digital_mpsk_snr_est_cc::set_tag_nsample(int n) +{ + if(n > 0) { + d_nsamples = n; + d_count = 0; // reset state + } + else + throw std::invalid_argument("digital_mpsk_snr_est_cc: tag_nsamples can't be <= 0\n"); +} + +void +digital_mpsk_snr_est_cc::set_alpha(double alpha) +{ + if((alpha >= 0) && (alpha <= 1.0)) { + d_alpha = alpha; + if(d_snr_est) + d_snr_est->set_alpha(d_alpha); + } + else + throw std::invalid_argument("digital_mpsk_snr_est_cc: alpha must be in [0,1]\n"); +} diff --git a/gr-digital/lib/digital_probe_mpsk_snr_est_c.cc b/gr-digital/lib/digital_probe_mpsk_snr_est_c.cc new file mode 100644 index 000000000..5cdfea96d --- /dev/null +++ b/gr-digital/lib/digital_probe_mpsk_snr_est_c.cc @@ -0,0 +1,152 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <digital_probe_mpsk_snr_est_c.h> +#include <gr_io_signature.h> +#include <cstdio> + +digital_probe_mpsk_snr_est_c_sptr +digital_make_probe_mpsk_snr_est_c(snr_est_type_t type, + int msg_nsamples, + double alpha) +{ + return gnuradio::get_initial_sptr( + new digital_probe_mpsk_snr_est_c(type, msg_nsamples, alpha)); +} + +digital_probe_mpsk_snr_est_c::digital_probe_mpsk_snr_est_c( + snr_est_type_t type, + int msg_nsamples, + double alpha) + : gr_sync_block ("probe_mpsk_snr_est_c", + gr_make_io_signature(1, 1, sizeof(gr_complex)), + gr_make_io_signature(0, 0, 0)) +{ + d_snr_est = NULL; + + d_type = type; + d_nsamples = msg_nsamples; + d_count = 0; + set_alpha(alpha); + + set_type(type); + + // at least 1 estimator has to look back + set_history(2); + + d_key = pmt::pmt_string_to_symbol("snr"); +} + +digital_probe_mpsk_snr_est_c::~digital_probe_mpsk_snr_est_c() +{ + if(d_snr_est) + delete d_snr_est; +} + +int +digital_probe_mpsk_snr_est_c::work(int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + const gr_complex *in = (const gr_complex*)input_items[0]; + return d_snr_est->update(noutput_items, in); +} + +double +digital_probe_mpsk_snr_est_c::snr() +{ + if(d_snr_est) + return d_snr_est->snr(); + else + throw std::runtime_error("digital_probe_mpsk_snr_est_c:: No SNR estimator defined.\n"); +} + +snr_est_type_t +digital_probe_mpsk_snr_est_c::type() const +{ + return d_type; +} + +int +digital_probe_mpsk_snr_est_c::msg_nsample() const +{ + return d_nsamples; +} + +double +digital_probe_mpsk_snr_est_c::alpha() const +{ + return d_alpha; +} + +void +digital_probe_mpsk_snr_est_c::set_type(snr_est_type_t t) +{ + d_type = t; + + if(d_snr_est) + delete d_snr_est; + + switch (d_type) { + case(SNR_EST_SIMPLE): + d_snr_est = new digital_impl_mpsk_snr_est_simple(d_alpha); + break; + case(SNR_EST_SKEW): + d_snr_est = new digital_impl_mpsk_snr_est_skew(d_alpha); + break; + case(SNR_EST_M2M4): + d_snr_est = new digital_impl_mpsk_snr_est_m2m4(d_alpha); + break; + case(SNR_EST_SVR): + d_snr_est = new digital_impl_mpsk_snr_est_svr(d_alpha); + break; + default: + throw std::invalid_argument("digital_probe_mpsk_snr_est_c: unknown type specified.\n"); + } +} + +void +digital_probe_mpsk_snr_est_c::set_msg_nsample(int n) +{ + if(n > 0) { + d_nsamples = n; + d_count = 0; // reset state + } + else + throw std::invalid_argument("digital_probe_mpsk_snr_est_c: msg_nsamples can't be <= 0\n"); +} + +void +digital_probe_mpsk_snr_est_c::set_alpha(double alpha) +{ + if((alpha >= 0) && (alpha <= 1.0)) { + d_alpha = alpha; + if(d_snr_est) + d_snr_est->set_alpha(d_alpha); + } + else + throw std::invalid_argument("digital_probe_mpsk_snr_est_c: alpha must be in [0,1]\n"); +} diff --git a/gr-digital/python/Makefile.am b/gr-digital/python/Makefile.am index ead6f7dfd..42bcc4dd2 100644 --- a/gr-digital/python/Makefile.am +++ b/gr-digital/python/Makefile.am @@ -45,6 +45,7 @@ noinst_PYTHON = \ qa_fll_band_edge.py \ qa_lms_equalizer.py \ qa_mpsk_receiver.py \ + qa_mpsk_snr_est.py \ qa_ofdm_insert_preamble.py digital_PYTHON = \ diff --git a/gr-digital/python/qa_fll_band_edge.py b/gr-digital/python/qa_fll_band_edge.py index 088eb2b68..7d89bc9ea 100755 --- a/gr-digital/python/qa_fll_band_edge.py +++ b/gr-digital/python/qa_fll_band_edge.py @@ -46,6 +46,7 @@ class test_fll_band_edge_cc(gr_unittest.TestCase): foffset = 0.2 / (2.0*math.pi) # Create a set of 1's and -1's, pulse shape and interpolate to sps + random.seed(0) data = [2.0*random.randint(0, 2) - 1.0 for i in xrange(200)] self.src = gr.vector_source_c(data, False) self.rrc = gr.interp_fir_filter_ccf(sps, rrc_taps) diff --git a/gr-digital/python/qa_mpsk_snr_est.py b/gr-digital/python/qa_mpsk_snr_est.py new file mode 100755 index 000000000..d392567bf --- /dev/null +++ b/gr-digital/python/qa_mpsk_snr_est.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python +# +# Copyright 2011 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. +# + +from gnuradio import gr, gr_unittest +import digital_swig as digital +import math, random + +def get_cplx(): + return complex(2*random.randint(0,1) - 1, 0) +def get_n_cplx(): + return complex(random.random()-0.5, random.random()-0.5) + +class test_mpsk_snr_est (gr_unittest.TestCase): + def setUp (self): + self.tb = gr.top_block () + + random.seed(0) # make repeatable + N = 10000 + self._noise = [get_n_cplx() for i in xrange(N)] + self._bits = [get_cplx() for i in xrange(N)] + + def tearDown (self): + self.tb = None + + def mpsk_snr_est_setup (self, op): + result = [] + for i in xrange(1,6): + src_data = [b+(i*n) for b,n in zip(self._bits, self._noise)] + + src = gr.vector_source_c (src_data) + dst = gr.null_sink (gr.sizeof_gr_complex) + + tb = gr.top_block () + tb.connect (src, op) + tb.connect (op, dst) + tb.run () # run the graph and wait for it to finish + + 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] + + N = 10000 + alpha = 0.001 + op = digital.mpsk_snr_est_cc (digital.SNR_EST_SIMPLE, N, alpha) + + actual_result = self.mpsk_snr_est_setup(op) + 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] + + N = 10000 + alpha = 0.001 + op = digital.mpsk_snr_est_cc (digital.SNR_EST_SKEW, N, alpha) + + actual_result = self.mpsk_snr_est_setup(op) + 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] + + N = 10000 + alpha = 0.001 + op = digital.mpsk_snr_est_cc (digital.SNR_EST_M2M4, N, alpha) + + actual_result = self.mpsk_snr_est_setup(op) + self.assertFloatTuplesAlmostEqual (expected_result, actual_result, 2) + + def test_mpsk_snr_est_svn (self): + expected_result = [10.90, 6.00, 4.76, 4.97, 5.49] + + N = 10000 + alpha = 0.001 + op = digital.mpsk_snr_est_cc (digital.SNR_EST_SVR, N, alpha) + + actual_result = self.mpsk_snr_est_setup(op) + 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] + + actual_result = [] + for i in xrange(1,6): + src_data = [b+(i*n) for b,n in zip(self._bits, self._noise)] + + src = gr.vector_source_c (src_data) + + N = 10000 + alpha = 0.001 + op = digital.probe_mpsk_snr_est_c (digital.SNR_EST_M2M4, N, alpha) + + tb = gr.top_block () + tb.connect (src, op) + tb.run () # run the graph and wait for it to finish + + actual_result.append(op.snr()) + self.assertFloatTuplesAlmostEqual (expected_result, actual_result, 2) + + +if __name__ == '__main__': + # Test various SNR estimators; we're not using a Gaussian + # 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") + diff --git a/gr-digital/swig/CMakeLists.txt b/gr-digital/swig/CMakeLists.txt index dd6097286..6f2c2251a 100644 --- a/gr-digital/swig/CMakeLists.txt +++ b/gr-digital/swig/CMakeLists.txt @@ -59,12 +59,14 @@ install( digital_lms_dd_equalizer_cc.i digital_kurtotic_equalizer_cc.i digital_mpsk_receiver_cc.i + digital_mpsk_snr_est_cc.i digital_ofdm_cyclic_prefixer.i digital_ofdm_frame_acquisition.i digital_ofdm_frame_sink.i digital_ofdm_insert_preamble.i digital_ofdm_mapper_bcv.i digital_ofdm_sampler.i + digital_probe_mpsk_snr_est_c.i digital_gmskmod_bc.i digital_cpmmod_bc.i DESTINATION ${GR_INCLUDE_DIR}/gnuradio/swig diff --git a/gr-digital/swig/Makefile.am b/gr-digital/swig/Makefile.am index 8591c8d74..97b47fafa 100644 --- a/gr-digital/swig/Makefile.am +++ b/gr-digital/swig/Makefile.am @@ -84,12 +84,14 @@ digital_swig_swiginclude_headers = \ digital_lms_dd_equalizer_cc.i \ digital_kurtotic_equalizer_cc.i \ digital_mpsk_receiver_cc.i \ + digital_mpsk_snr_est_cc.i \ digital_ofdm_cyclic_prefixer.i \ digital_ofdm_frame_acquisition.i \ digital_ofdm_frame_sink.i \ digital_ofdm_insert_preamble.i \ digital_ofdm_mapper_bcv.i \ digital_ofdm_sampler.i \ + digital_probe_mpsk_snr_est_c.i \ digital_gmskmod_bc.i \ digital_cpmmod_bc.i \ $(TOP_SWIG_DOC_IFILES) diff --git a/gr-digital/swig/digital_mpsk_snr_est_cc.i b/gr-digital/swig/digital_mpsk_snr_est_cc.i new file mode 100644 index 000000000..f0ca13f87 --- /dev/null +++ b/gr-digital/swig/digital_mpsk_snr_est_cc.i @@ -0,0 +1,45 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ + +GR_SWIG_BLOCK_MAGIC(digital,mpsk_snr_est_cc); + +digital_mpsk_snr_est_cc_sptr +digital_make_mpsk_snr_est_cc(snr_est_type_t type, + int tag_nsamples=10000, + double alpha=0.001); + +class digital_mpsk_snr_est_cc : public gr_sync_block +{ +private: + void digital_mpsk_snr_est_cc(snr_est_type_t type, + int tag_nsamples, + double alpha); + +public: + double snr(); + snr_est_type_t type() const; + int tag_nsample() const; + double alpha() const; + void set_type(snr_est_type_t t); + void set_tag_nsample(int n); + void set_alpha(double alpha); +}; diff --git a/gr-digital/swig/digital_probe_mpsk_snr_est_c.i b/gr-digital/swig/digital_probe_mpsk_snr_est_c.i new file mode 100644 index 000000000..93db4127a --- /dev/null +++ b/gr-digital/swig/digital_probe_mpsk_snr_est_c.i @@ -0,0 +1,45 @@ +/* -*- c++ -*- */ +/* + * Copyright 2011 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. + */ + +GR_SWIG_BLOCK_MAGIC(digital,probe_mpsk_snr_est_c); + +digital_probe_mpsk_snr_est_c_sptr +digital_make_probe_mpsk_snr_est_c(snr_est_type_t type, + int msg_nsamples=10000, + double alpha=0.001); + +class digital_probe_mpsk_snr_est_c : public gr_sync_block +{ +private: + void digital_probe_mpsk_snr_est_c(snr_est_type_t type, + int msg_nsamples, + double alpha); + +public: + double snr(); + snr_est_type_t type() const; + int msg_nsample() const; + double alpha() const; + void set_type(snr_est_type_t t); + void set_msg_nsample(int n); + void set_alpha(double alpha); +}; diff --git a/gr-digital/swig/digital_swig.i b/gr-digital/swig/digital_swig.i index 86b5cab13..a39ef9ab7 100644 --- a/gr-digital/swig/digital_swig.i +++ b/gr-digital/swig/digital_swig.i @@ -24,6 +24,15 @@ //load generated python docstrings %include "digital_swig_doc.i" +#if SWIGPYTHON +enum snr_est_type_t { + SNR_EST_SIMPLE = 0, // Simple estimator (>= 7 dB) + SNR_EST_SKEW, // Skewness-base est (>= 5 dB) + SNR_EST_M2M4, // 2nd & 4th moment est (>= 1 dB) + SNR_EST_SVR // SVR-based est (>= 0dB) +}; +#endif + %include <gri_control_loop.i> %{ @@ -41,12 +50,14 @@ #include "digital_kurtotic_equalizer_cc.h" #include "digital_lms_dd_equalizer_cc.h" #include "digital_mpsk_receiver_cc.h" +#include "digital_mpsk_snr_est_cc.h" #include "digital_ofdm_cyclic_prefixer.h" #include "digital_ofdm_frame_acquisition.h" #include "digital_ofdm_frame_sink.h" #include "digital_ofdm_insert_preamble.h" #include "digital_ofdm_mapper_bcv.h" #include "digital_ofdm_sampler.h" +#include "digital_probe_mpsk_snr_est_c.h" #include "digital_cpmmod_bc.h" #include "digital_gmskmod_bc.h" %} @@ -65,16 +76,26 @@ %include "digital_kurtotic_equalizer_cc.i" %include "digital_lms_dd_equalizer_cc.i" %include "digital_mpsk_receiver_cc.i" +%include "digital_mpsk_snr_est_cc.i" %include "digital_ofdm_cyclic_prefixer.i" %include "digital_ofdm_frame_acquisition.i" %include "digital_ofdm_frame_sink.i" %include "digital_ofdm_insert_preamble.i" %include "digital_ofdm_mapper_bcv.i" %include "digital_ofdm_sampler.i" +%include "digital_probe_mpsk_snr_est_c.i" %include "digital_cpmmod_bc.i" %include "digital_gmskmod_bc.i" #if SWIGGUILE + +enum snr_est_type_t { + SNR_EST_SIMPLE = 0, // Simple estimator (>= 7 dB) + SNR_EST_SKEW, // Skewness-base est (>= 5 dB) + SNR_EST_M2M4, // 2nd & 4th moment est (>= 1 dB) + SNR_EST_SVR // SVR-based est (>= 0dB) +}; + %scheme %{ (load-extension-global "libguile-gnuradio-digital_swig" "scm_init_gnuradio_digital_swig_module") %} |