[Commit-gnuradio] r4618 - gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general

[trondeau]

Author: trondeau
Date: 2007-02-24 14:06:57 -0700 (Sat, 24 Feb 2007)
New Revision: 4618

Modified:
   
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.cc
   
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.h
Log:
corrections to the mpsk receiver and some cleanup to make it more readable.


Modified: 
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.cc
===================================================================
--- 
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.cc
     2007-02-24 03:15:09 UTC (rev 4617)
+++ 
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.cc
     2007-02-24 21:06:57 UTC (rev 4618)
@@ -1,6 +1,6 @@
 /* -*- c++ -*- */
 /*
- * Copyright 2005,2006 Free Software Foundation, Inc.
+ * Copyright 2005,2006,2007 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -34,8 +34,9 @@
 
 
 #define M_TWOPI (2*M_PI)
+#define VERBOSE_MM     0     // Used for debugging symbol timing loop
+#define VERBOSE_COSTAS 0     // Used for debugging phase and frequency tracking
 
-
 // Public constructor
 
 gr_mpsk_receiver_cc_sptr 
@@ -62,6 +63,7 @@
              gr_make_io_signature (1, 1, sizeof (gr_complex))),
     d_M(M), d_theta(theta), 
     d_alpha(alpha), d_beta(beta), d_freq(0), d_max_freq(fmax), 
d_min_freq(fmin), d_phase(0),
+    d_current_const_point(0),
     d_mu(mu), d_gain_mu(gain_mu), d_gain_omega(gain_omega), 
     d_omega_rel(omega_rel), d_max_omega(0), d_min_omega(0),
     d_p_2T(0), d_p_1T(0), d_p_0T(0), d_c_2T(0), d_c_1T(0), d_c_0T(0)
@@ -82,17 +84,19 @@
   for (unsigned int i = 0; i < 2 * DLLEN; i++)
     d_dl[i] = gr_complex(0.0,0.0);
 
-
+  // build the constellation vector from M
+  make_constellation();
+  
   // Select a phase detector and a decision maker for the modulation order
   switch(d_M) {
   case 2:  // optimized algorithms for BPSK
-    d_phase_error_detector = &gr_mpsk_receiver_cc::phase_error_detector_bpsk;
-    d_decision = &gr_mpsk_receiver_cc::decision_bpsk;
+    d_phase_error_detector = 
&gr_mpsk_receiver_cc::phase_error_detector_generic; //bpsk;
+    d_decision = &gr_mpsk_receiver_cc::decision_generic; //bpsk;
     break;
 
   case 4: // optimized algorithms for QPSK
-    d_phase_error_detector = &gr_mpsk_receiver_cc::phase_error_detector_qpsk;
-    d_decision = &gr_mpsk_receiver_cc::decision_qpsk;
+    d_phase_error_detector = 
&gr_mpsk_receiver_cc::phase_error_detector_generic; //qpsk;
+    d_decision = &gr_mpsk_receiver_cc::decision_generic; //qpsk;
     break;
 
   default: // generic algorithms for any M (power of 2?) but not pretty
@@ -100,6 +104,8 @@
     d_decision = &gr_mpsk_receiver_cc::decision_generic;
     break;
   }
+
+  set_history(3);                      // ensure 2 extra input sample is 
available
 }
 
 gr_mpsk_receiver_cc::~gr_mpsk_receiver_cc ()
@@ -112,14 +118,17 @@
 {
   unsigned ninputs = ninput_items_required.size();
   for (unsigned i=0; i < ninputs; i++)
-    ninput_items_required[i] = (int)(d_omega);
+    ninput_items_required[i] = (int) ceil((noutput_items * d_omega) + 
d_interp->ntaps());
+  //ninput_items_required[i] = (int)(d_omega);
+
 }
 
 float
 gr_mpsk_receiver_cc::phase_error_detector_qpsk(gr_complex sample) const
 {
-  return ((sample.real()>0 ? 1.0 : -1.0) * sample.imag() -
-         (sample.imag()>0 ? 1.0 : -1.0) * sample.real());
+  float phase_error = ((sample.real()>0 ? 1.0 : -1.0) * sample.imag() -
+                      (sample.imag()>0 ? 1.0 : -1.0) * sample.real());
+  return phase_error;
 }
 
 float
@@ -130,77 +139,165 @@
 
 float gr_mpsk_receiver_cc::phase_error_detector_generic(gr_complex sample) 
const
 {
-  unsigned int min_m = 0;
-  gr_complex constpoint;
-  float min_s = 65535;
-  float error = 0;
-
-  // Develop all possible constellation points and find the one that minimizes
-  // the Euclidean distance (error) with the sample
-  for(unsigned int m=0; m < d_M; m++) {
-    constpoint = gr_expj((M_TWOPI/d_M)*m + d_theta);
-    gr_complex diff = (constpoint - sample)*(constpoint - sample);
-    
-    if(fabs(diff.real()) < min_s) {
-      min_s = fabs(diff.real());
-      min_m = m;
-    }
-  }
-  // Use this constellation point to derotate the sample
-  constpoint = gr_expj(-((M_TWOPI/d_M)*min_m + d_theta));
-  // get the error in the sample as the phase difference from the real axis
-  error = gr_fast_atan2f(sample*constpoint);
-
-  return error;
+  //return gr_fast_atan2f(sample*conj(d_constellation[d_current_const_point]));
+  return -arg(sample*conj(d_constellation[d_current_const_point]));
 }
 
-gr_complex
+unsigned int
 gr_mpsk_receiver_cc::decision_bpsk(gr_complex sample) const
 {
-  float real=-0.5, imag=0.0;
+  unsigned int index = 0;
 
   // Implements a 1-demensional slicer
   if(sample.real() > 0)
-    real = 0.5;
-  return gr_complex(real,imag);
+    index = 1;
+  return index;
 }
 
-gr_complex
+unsigned int
 gr_mpsk_receiver_cc::decision_qpsk(gr_complex sample) const
 {
-  float real=-0.5, imag=-0.5;
+  unsigned int index = 0;
 
   // Implements a simple slicer function
-  if(sample.real() > 0)
-    real = 0.5;
-  if(sample.imag() > 0)
-    imag = 0.5;
-  return gr_complex(real,imag);
+  if((sample.real() < 0) && (sample.imag() > 0))
+    index = 1;
+  else if((sample.real() < 0) && (sample.imag() < 0))
+    index = 2;
+  else
+    index = 3;
+  return index;
 }
 
-gr_complex
+unsigned int
 gr_mpsk_receiver_cc::decision_generic(gr_complex sample) const
 {
   unsigned int min_m = 0;
-  gr_complex constpoint;
   float min_s = 65535;
   
   // Develop all possible constellation points and find the one that minimizes
   // the Euclidean distance (error) with the sample
   for(unsigned int m=0; m < d_M; m++) {
-    constpoint = gr_expj((M_TWOPI/d_M)*m + d_theta);
-    gr_complex diff = (constpoint - sample)*(constpoint - sample);
+    gr_complex diff = ( d_constellation[m] - sample)*( d_constellation[m] - 
sample);
     
     if(fabs(diff.real()) < min_s) {
       min_s = fabs(diff.real());
       min_m = m;
     }
   }
-  // Return the constellation point that minimizes the error
-  constpoint = gr_expj(((M_TWOPI/d_M)*min_m + d_theta));
-  return constpoint;
+  // Return the index of the constellation point that minimizes the error
+  return min_m;
 }
 
+
+void
+gr_mpsk_receiver_cc::make_constellation()
+{
+  for(unsigned int m=0; m < d_M; m++) {
+    d_constellation.push_back(gr_expj((M_TWOPI/d_M)*m));
+  }
+}
+
+void
+gr_mpsk_receiver_cc::mm_sampler(const gr_complex symbol)
+{
+  gr_complex sample, nco;
+
+  d_mu--;             // skip a number of symbols between sampling
+  d_phase += d_freq;  // increment the phase based on the frequency of the 
rotation
+
+  // Keep phase clamped and not walk to infinity
+  while(d_phase>M_TWOPI)
+    d_phase -= M_TWOPI;
+  while(d_phase<-M_TWOPI)
+    d_phase += M_TWOPI;
+
+  nco = gr_expj(d_phase+d_theta);   // get the NCO value for derotating the 
current sample
+  sample = nco*symbol;      // get the downconverted symbol
+  
+  // Fill up the delay line for the interpolator
+  d_dl[d_dl_idx] = sample;
+  d_dl[(d_dl_idx + DLLEN)] = sample;  // put this in the second half of the 
buffer for overflows
+  d_dl_idx = (d_dl_idx+1) % DLLEN;    // Keep the delay line index in bounds
+}
+
+void
+gr_mpsk_receiver_cc::mm_error_tracking(gr_complex sample)
+{
+  gr_complex u, x, y;
+  float mm_error = 0;
+
+  // Make sample timing corrections
+
+  // set the delayed samples
+  d_p_2T = d_p_1T;
+  d_p_1T = d_p_0T;
+  d_p_0T = sample;
+  d_c_2T = d_c_1T;
+  d_c_1T = d_c_0T;
+  
+  d_current_const_point = (*this.*d_decision)(d_p_0T);  // make a decision on 
the sample value
+  d_c_0T = d_constellation[d_current_const_point];
+  
+  x = (d_c_0T - d_c_2T) * conj(d_p_1T);
+  y = (d_p_0T - d_p_2T) * conj(d_c_1T);
+  u = y - x;
+  mm_error = u.real();   // the error signal is in the real part
+  
+  // limit mm_val
+  if (mm_error > 1.0)
+    mm_error = 1.0;
+  else if (mm_error < -1.0)
+    mm_error = -1.0;
+  
+  d_omega = d_omega + d_gain_omega * mm_error;  // update omega based on loop 
error
+
+  // make sure we don't walk away
+  if (d_omega > d_max_omega)
+    d_omega = d_max_omega;
+  else if (d_omega < d_min_omega)
+    d_omega = d_min_omega;
+  
+  d_mu += d_omega + d_gain_mu * mm_error;   // update mu based on loop error
+  
+#if VERBOSE_MM
+  printf("mm: mu: %f   omega: %f  mm_error: %f  sample: %f+j%f  constellation: 
%f+j%f\n", 
+        d_mu, d_omega, mm_error, sample.real(), sample.imag(), 
+        d_constellation[d_current_const_point].real(), 
d_constellation[d_current_const_point].imag());
+#endif
+}
+
+
+void
+gr_mpsk_receiver_cc::phase_error_tracking(gr_complex sample)
+{
+  float phase_error = 0;
+
+  // Make phase and frequency corrections based on sampled value
+  phase_error = (*this.*d_phase_error_detector)(sample);
+  
+  d_freq += d_beta*phase_error;             // adjust frequency based on error
+  d_phase += d_freq + d_alpha*phase_error;  // adjust phase based on error
+  
+  // Make sure we stay within +-2pi
+  while(d_phase>M_TWOPI)
+    d_phase -= M_TWOPI;
+  while(d_phase<-M_TWOPI)
+    d_phase += M_TWOPI;
+  
+  // Limit the frequency range
+  if (d_freq > d_max_freq)
+    d_freq = d_max_freq;
+  else if (d_freq < d_min_freq)
+    d_freq = d_min_freq;
+
+#if VERBOSE_COSTAS
+  printf("cl: phase_error: %f  phase: %f  freq: %f  sample: %f+j%f  
constellation: %f+j%f\n",
+        phase_error, d_phase, d_freq, sample.real(), sample.imag(), 
+        d_constellation[d_current_const_point].real(), 
d_constellation[d_current_const_point].imag());
+#endif
+}
+
 int
 gr_mpsk_receiver_cc::general_work (int noutput_items,
                                   gr_vector_int &ninput_items,
@@ -210,88 +307,30 @@
   const gr_complex *in = (const gr_complex *) input_items[0];
   gr_complex *out = (gr_complex *) output_items[0];
 
-  float phase_error = 0, mm_error = 0;
-  gr_complex sample, nco, u, x, y;
-  
   int i=0, o=0;
 
-  while(i < ninput_items[0]) {    
+  //while(i < ninput_items[0]) {    
+  while(o < noutput_items) {
     while((d_mu > 1) && (i < ninput_items[0]))  {
-      d_mu--;             // skip a number of symbols between sampling
-      d_phase += d_freq;  // increment the phase based on the frequency of the 
rotation
-
-      // Keep phase clamped and not walk to infinity
-      while(d_phase>M_TWOPI)
-       d_phase -= M_TWOPI;
-      while(d_phase<-M_TWOPI)
-       d_phase += M_TWOPI;
-
-      nco = gr_expj(-d_phase); // get the NCO value for derotating the current 
sample
-      sample = nco*in[i];      // get the downconverted symbol
-
-      // Fill up the delay line for the interpolator
-      d_dl[d_dl_idx] = sample;
-      d_dl[(d_dl_idx + DLLEN)] = sample;  // put this in the second half of 
the buffer for overflows
-      d_dl_idx = (d_dl_idx+1) % DLLEN;    // Keep the delay line index in 
bounds
-
+      mm_sampler(in[i]);   // puts symbols into a buffer and adjusts d_mu
       i++;
     }
     
     if(i < ninput_items[0]) {
       gr_complex interp_sample = d_interp->interpolate(&d_dl[d_dl_idx], d_mu);
        
-      // Make sample timing corrections
-      d_p_2T = d_p_1T;
-      d_p_1T = d_p_0T;
-      d_p_0T = interp_sample;
+      mm_error_tracking(interp_sample);     // corrects M&M sample time
+      phase_error_tracking(interp_sample);  // corrects phase and frequency 
offsets
 
-      d_c_2T = d_c_1T;
-      d_c_1T = d_c_0T;
-      d_c_0T = (*this.*d_decision)(d_p_0T);
-      
-      x = (d_c_0T - d_c_2T) * conj(d_p_1T);
-      y = (d_p_0T - d_p_2T) * conj(d_c_1T);
-      u = y - x;
-      mm_error = u.real();
-
-      // limit mm_val
-      if (mm_error > 1.0)
-       mm_error = 1.0;
-      else if (mm_error < -1.0)
-       mm_error = -1.0;
-      
-      d_omega = d_omega + d_gain_omega * mm_error;
-      if (d_omega > d_max_omega)
-       d_omega = d_max_omega;
-      else if (d_omega < d_min_omega)
-       d_omega = d_min_omega;
-      
-      d_mu += d_omega + d_gain_mu * mm_error;
-      
-      #if 0
-      printf("samples: %f+j%f   mu: %f   omega: %f  mm_error: %f\n", 
-            interp_sample.real(), interp_sample.imag(), d_mu, d_omega, 
mm_error);
-      #endif
-
-      // Make phase and frequency corrections based on sampled value
-      phase_error = (*this.*d_phase_error_detector)(interp_sample);
-      d_freq += d_beta*phase_error;
-      d_phase += d_alpha*phase_error;
-
-      while(d_phase>M_TWOPI)
-       d_phase -= M_TWOPI;
-      while(d_phase<-M_TWOPI)
-       d_phase += M_TWOPI;
-
-      if (d_freq > d_max_freq)
-       d_freq = d_max_freq;
-      else if (d_freq < d_min_freq)
-       d_freq = d_min_freq;
-
-      out[o++] = d_p_0T;
+      out[o++] = interp_sample;
     }
   }
 
+  #if 0
+  printf("ninput_items: %d   noutput_items: %d   consuming: %d   returning: 
%d\n",
+        ninput_items[0], noutput_items, i, o);
+  #endif
+
   consume_each(i);
   return o;
 }

Modified: 
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.h
===================================================================
--- 
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.h
      2007-02-24 03:15:09 UTC (rev 4617)
+++ 
gnuradio/branches/developers/trondeau/digital-wip2/gnuradio-core/src/lib/general/gr_mpsk_receiver_cc.h
      2007-02-24 21:06:57 UTC (rev 4618)
@@ -157,11 +157,17 @@
    * value of M.
    */
   gr_mpsk_receiver_cc (unsigned int M, float theta, 
-                               float alpha, float beta,
-                               float fmin, float fmax,
-                               float mu, float gain_mu, 
-                               float omega, float gain_omega, float omega_rel);
+                      float alpha, float beta,
+                      float fmin, float fmax,
+                      float mu, float gain_mu, 
+                      float omega, float gain_omega, float omega_rel);
 
+  void make_constellation();
+  void mm_sampler(const gr_complex symbol);
+  void mm_error_tracking(gr_complex sample);
+  void phase_error_tracking(gr_complex sample);
+
+
 /*!
    * \brief Phase error detector for MPSK modulations.
    *
@@ -200,21 +206,20 @@
    * \returns the approximated phase error.
  */
   float phase_error_detector_qpsk(gr_complex sample) const;
- 
 
 
+
  /*!
    * \brief Decision maker for a generic MPSK constellation.
    *
    * \param sample   the baseband I&Q sample from which to make the decision
    *
-   * This decision maker is a generic implementation the creates a set of 
constellation point for MPSK
-   * and does a brute-force search for the constellation point that minimizes 
the error from the incoming
-   * signal.
+   * This decision maker is a generic implementation that does a brute-force 
search 
+   * for the constellation point that minimizes the error between it and the 
incoming signal.
    *
-   * \returns the constellation point that minimizes error.
+   * \returns the index to d_constellation that minimizes the error/
  */
-  gr_complex decision_generic(gr_complex sample) const;
+  unsigned int decision_generic(gr_complex sample) const;
 
 
  /*!
@@ -225,9 +230,9 @@
    * This decision maker is a simple slicer function that makes a decision on 
the symbol based on its
    * placement on the real axis of greater than 0 or less than 0; the 
quadrature component is always 0.
    *
-   * \returns the constellation point that minimizes error.
+   * \returns the index to d_constellation that minimizes the error/
  */
-  gr_complex decision_bpsk(gr_complex sample) const;
+  unsigned int decision_bpsk(gr_complex sample) const;
   
 
  /*!
@@ -238,9 +243,9 @@
    * This decision maker is a simple slicer function that makes a decision on 
the symbol based on its
    * placement versus both axes and returns which quadrant the symbol is in.
    *
-   * \returns the constellation point that minimizes error.
+   * \returns the index to d_constellation that minimizes the error/
  */
-  gr_complex decision_qpsk(gr_complex sample) const;
+  unsigned int decision_qpsk(gr_complex sample) const;
 
   private:
   unsigned int d_M;
@@ -259,10 +264,14 @@
    *
    * This is a function pointer that is set in the constructor to point to the 
proper decision function
    * for the specified constellation order.
+   *
+   * \return index into d_constellation point that is the closest to the 
recieved sample
  */
-  gr_complex (gr_mpsk_receiver_cc::*d_decision)(gr_complex sample) const; // 
pointer to decision function
+  unsigned int (gr_mpsk_receiver_cc::*d_decision)(gr_complex sample) const; // 
pointer to decision function
 
 
+  std::vector<gr_complex> d_constellation;
+  unsigned int d_current_const_point;
 
   // Members related to symbol timing
   float d_mu, d_gain_mu;