git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2005,2006,2007 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 #ifdef HAVE_CONFIG_H
  24 #include "config.h"
  25 #endif
  26
  27 #include <gr_io_signature.h>
  28 #include <gr_prefs.h>
  29 #include <gr_mpsk_receiver_cc.h>
  30 #include <stdexcept>
  31 #include <gr_math.h>
  32 #include <gr_expj.h>
  33 #include <gri_mmse_fir_interpolator_cc.h>
  34
  35
  36 #define M_TWOPI (2*M_PI)
  37 #define VERBOSE_MM     0     // Used for debugging symbol timing loop
  38 #define VERBOSE_COSTAS 0     // Used for debugging phase and frequency tracking
  39
  40 // Public constructor
  41
  42 gr_mpsk_receiver_cc_sptr
  43 gr_make_mpsk_receiver_cc(unsigned int M, float theta,
  44                          float alpha, float beta,
  45                          float fmin, float fmax,
  46                          float mu, float gain_mu,
  47                          float omega, float gain_omega, float omega_rel)
  48 {
  49   return gr_mpsk_receiver_cc_sptr (new gr_mpsk_receiver_cc (M, theta,
  50                                                             alpha, beta,
  51                                                             fmin, fmax,
  52                                                             mu, gain_mu,
  53                                                             omega, gain_omega, omega_rel));
  54 }
  55
  56 gr_mpsk_receiver_cc::gr_mpsk_receiver_cc (unsigned int M, float theta,
  57                                           float alpha, float beta,
  58                                           float fmin, float fmax,
  59                                           float mu, float gain_mu,
  60                                           float omega, float gain_omega, float omega_rel)
  61   : gr_block ("mpsk_receiver_cc",
  62               gr_make_io_signature (1, 1, sizeof (gr_complex)),
  63               gr_make_io_signature (1, 1, sizeof (gr_complex))),
  64     d_M(M), d_theta(theta),
  65     d_alpha(alpha), d_beta(beta), d_freq(0), d_max_freq(fmax), d_min_freq(fmin), d_phase(0),
  66     d_current_const_point(0),
  67     d_mu(mu), d_gain_mu(gain_mu), d_gain_omega(gain_omega),
  68     d_omega_rel(omega_rel), d_max_omega(0), d_min_omega(0),
  69     d_p_2T(0), d_p_1T(0), d_p_0T(0), d_c_2T(0), d_c_1T(0), d_c_0T(0)
  70 {
  71   d_interp = new gri_mmse_fir_interpolator_cc();
  72   d_dl_idx = 0;
  73
  74   set_omega(omega);
  75
  76   if (omega <= 0.0)
  77     throw std::out_of_range ("clock rate must be > 0");
  78   if (gain_mu <  0  || gain_omega < 0)
  79     throw std::out_of_range ("Gains must be non-negative");
  80
  81   assert(d_interp->ntaps() <= DLLEN);
  82
  83   // zero double length delay line.
  84   for (unsigned int i = 0; i < 2 * DLLEN; i++)
  85     d_dl[i] = gr_complex(0.0,0.0);
  86
  87   // build the constellation vector from M
  88   make_constellation();
  89
  90   // Select a phase detector and a decision maker for the modulation order
  91   switch(d_M) {
  92   case 2:  // optimized algorithms for BPSK
  93     d_phase_error_detector = &gr_mpsk_receiver_cc::phase_error_detector_bpsk; //bpsk;
  94     d_decision = &gr_mpsk_receiver_cc::decision_bpsk;
  95     break;
  96
  97   case 4: // optimized algorithms for QPSK
  98     d_phase_error_detector = &gr_mpsk_receiver_cc::phase_error_detector_qpsk; //qpsk;
  99     d_decision = &gr_mpsk_receiver_cc::decision_qpsk;
 100     break;
 101
 102   default: // generic algorithms for any M (power of 2?) but not pretty
 103     d_phase_error_detector = &gr_mpsk_receiver_cc::phase_error_detector_generic;
 104     d_decision = &gr_mpsk_receiver_cc::decision_generic;
 105     break;
 106   }
 107 }
 108
 109 gr_mpsk_receiver_cc::~gr_mpsk_receiver_cc ()
 110 {
 111   delete d_interp;
 112 }
 113
 114 void
 115 gr_mpsk_receiver_cc::forecast(int noutput_items, gr_vector_int &ninput_items_required)
 116 {
 117   unsigned ninputs = ninput_items_required.size();
 118   for (unsigned i=0; i < ninputs; i++)
 119     ninput_items_required[i] = (int) ceil((noutput_items * d_omega) + d_interp->ntaps());
 120 }
 121
 122 // FIXME add these back in an test difference in performance
 123 float
 124 gr_mpsk_receiver_cc::phase_error_detector_qpsk(gr_complex sample) const
 125 {
 126   float phase_error = 0;
 127   if(fabsf(sample.real()) > fabsf(sample.imag())) {
 128     if(sample.real() > 0)
 129       phase_error = -sample.imag();
 130     else
 131       phase_error = sample.imag();
 132   }
 133   else {
 134     if(sample.imag() > 0)
 135       phase_error = sample.real();
 136     else
 137       phase_error = -sample.real();
 138   }
 139
 140   return phase_error;
 141 }
 142
 143 float
 144 gr_mpsk_receiver_cc::phase_error_detector_bpsk(gr_complex sample) const
 145 {
 146   return -(sample.real()*sample.imag());
 147 }
 148
 149 float gr_mpsk_receiver_cc::phase_error_detector_generic(gr_complex sample) const
 150 {
 151   //return gr_fast_atan2f(sample*conj(d_constellation[d_current_const_point]));
 152   return -arg(sample*conj(d_constellation[d_current_const_point]));
 153 }
 154
 155 unsigned int
 156 gr_mpsk_receiver_cc::decision_bpsk(gr_complex sample) const
 157 {
 158   return (gr_branchless_binary_slicer(sample.real()) ^ 1);
 159   //return gr_binary_slicer(sample.real()) ^ 1;
 160 }
 161
 162 unsigned int
 163 gr_mpsk_receiver_cc::decision_qpsk(gr_complex sample) const
 164 {
 165   unsigned int index;
 166
 167   //index = gr_branchless_quad_0deg_slicer(sample);
 168   index = gr_quad_0deg_slicer(sample);
 169   return index;
 170 }
 171
 172 unsigned int
 173 gr_mpsk_receiver_cc::decision_generic(gr_complex sample) const
 174 {
 175   unsigned int min_m = 0;
 176   float min_s = 65535;
 177
 178   // Develop all possible constellation points and find the one that minimizes
 179   // the Euclidean distance (error) with the sample
 180   for(unsigned int m=0; m < d_M; m++) {
 181     gr_complex diff = norm(d_constellation[m] - sample);
 182
 183     if(fabs(diff.real()) < min_s) {
 184       min_s = fabs(diff.real());
 185       min_m = m;
 186     }
 187   }
 188   // Return the index of the constellation point that minimizes the error
 189   return min_m;
 190 }
 191
 192
 193 void
 194 gr_mpsk_receiver_cc::make_constellation()
 195 {
 196   for(unsigned int m=0; m < d_M; m++) {
 197     d_constellation.push_back(gr_expj((M_TWOPI/d_M)*m));
 198   }
 199 }
 200
 201 void
 202 gr_mpsk_receiver_cc::mm_sampler(const gr_complex symbol)
 203 {
 204   gr_complex sample, nco;
 205
 206   d_mu--;             // skip a number of symbols between sampling
 207   d_phase += d_freq;  // increment the phase based on the frequency of the rotation
 208
 209   // Keep phase clamped and not walk to infinity
 210   while(d_phase > M_TWOPI)
 211     d_phase -= M_TWOPI;
 212   while(d_phase < -M_TWOPI)
 213     d_phase += M_TWOPI;
 214
 215   nco = gr_expj(d_phase+d_theta);   // get the NCO value for derotating the current sample
 216   sample = nco*symbol;      // get the downconverted symbol
 217
 218   // Fill up the delay line for the interpolator
 219   d_dl[d_dl_idx] = sample;
 220   d_dl[(d_dl_idx + DLLEN)] = sample;  // put this in the second half of the buffer for overflows
 221   d_dl_idx = (d_dl_idx+1) % DLLEN;    // Keep the delay line index in bounds
 222 }
 223
 224 void
 225 gr_mpsk_receiver_cc::mm_error_tracking(gr_complex sample)
 226 {
 227   gr_complex u, x, y;
 228   float mm_error = 0;
 229
 230   // Make sample timing corrections
 231
 232   // set the delayed samples
 233   d_p_2T = d_p_1T;
 234   d_p_1T = d_p_0T;
 235   d_p_0T = sample;
 236   d_c_2T = d_c_1T;
 237   d_c_1T = d_c_0T;
 238
 239   d_current_const_point = (*this.*d_decision)(d_p_0T);  // make a decision on the sample value
 240   d_c_0T = d_constellation[d_current_const_point];
 241
 242   x = (d_c_0T - d_c_2T) * conj(d_p_1T);
 243   y = (d_p_0T - d_p_2T) * conj(d_c_1T);
 244   u = y - x;
 245   mm_error = u.real();   // the error signal is in the real part
 246   mm_error = gr_branchless_clip(mm_error, 1.0); // limit mm_val
 247
 248   d_omega = d_omega + d_gain_omega * mm_error;  // update omega based on loop error
 249   d_omega = d_omega_mid + gr_branchless_clip(d_omega-d_omega_mid, d_omega_rel);   // make sure we don't walk away
 250
 251   d_mu += d_omega + d_gain_mu * mm_error;   // update mu based on loop error
 252
 253 #if VERBOSE_MM
 254   printf("mm: mu: %f   omega: %f  mm_error: %f  sample: %f+j%f  constellation: %f+j%f\n",
 255          d_mu, d_omega, mm_error, sample.real(), sample.imag(),
 256          d_constellation[d_current_const_point].real(), d_constellation[d_current_const_point].imag());
 257 #endif
 258 }
 259
 260
 261 void
 262 gr_mpsk_receiver_cc::phase_error_tracking(gr_complex sample)
 263 {
 264   float phase_error = 0;
 265
 266   // Make phase and frequency corrections based on sampled value
 267   phase_error = (*this.*d_phase_error_detector)(sample);
 268
 269   phase_error = gr_branchless_clip(phase_error, 1.0);
 270
 271   d_freq += d_beta*phase_error;             // adjust frequency based on error
 272   d_phase += d_freq + d_alpha*phase_error;  // adjust phase based on error
 273
 274   // Make sure we stay within +-2pi
 275   while(d_phase > M_TWOPI)
 276     d_phase -= M_TWOPI;
 277   while(d_phase < -M_TWOPI)
 278     d_phase += M_TWOPI;
 279
 280   // Limit the frequency range
 281   d_freq = gr_branchless_clip(d_freq, d_max_freq);
 282
 283 #if VERBOSE_COSTAS
 284   printf("cl: phase_error: %f  phase: %f  freq: %f  sample: %f+j%f  constellation: %f+j%f\n",
 285          phase_error, d_phase, d_freq, sample.real(), sample.imag(),
 286          d_constellation[d_current_const_point].real(), d_constellation[d_current_const_point].imag());
 287 #endif
 288 }
 289
 290 int
 291 gr_mpsk_receiver_cc::general_work (int noutput_items,
 292                                    gr_vector_int &ninput_items,
 293                                    gr_vector_const_void_star &input_items,
 294                                    gr_vector_void_star &output_items)
 295 {
 296   const gr_complex *in = (const gr_complex *) input_items[0];
 297   gr_complex *out = (gr_complex *) output_items[0];
 298
 299   int i=0, o=0;
 300
 301   while((o < noutput_items) && (i < ninput_items[0])) {
 302     while((d_mu > 1) && (i < ninput_items[0]))  {
 303       mm_sampler(in[i]);   // puts symbols into a buffer and adjusts d_mu
 304       i++;
 305     }
 306
 307     if(i < ninput_items[0]) {
 308       gr_complex interp_sample = d_interp->interpolate(&d_dl[d_dl_idx], d_mu);
 309
 310       mm_error_tracking(interp_sample);     // corrects M&M sample time
 311       phase_error_tracking(interp_sample);  // corrects phase and frequency offsets
 312
 313       out[o++] = interp_sample;
 314     }
 315   }
 316
 317   #if 0
 318   printf("ninput_items: %d   noutput_items: %d   consuming: %d   returning: %d\n",
 319          ninput_items[0], noutput_items, i, o);
 320   #endif
 321
 322   consume_each(i);
 323   return o;
 324 }