[Patch-gnuradio] gr_clock_recovery_mm_cc.cc/h

[Mattias Kjellsson]

Hi,

I rewrote gr_clock_recovery_cc to use either of the two slicers
depending on the new parameter "order" (2/4, 2 per default) as is done
in the costas- loop. Although I haven't incorporated code for rotating
the constellation, if it's supposed to look like a "+", instead of the
slicers assumed "x".

I read the article referenced in the source, and I don't think there
should be any problems, but it might be something to review by people
with a better understanding than me. I have tested the block in a radio-
receiver, and by looking at the constellation, it "looks like it might
do what it's supposed to" (making the received "o"- constellation shape,
into a "x"- shape).

Anyway, I thought I might post the code (although there might be room
for improvement, as mentioned above), and wait for reactions.

Best regards,
//Mattias Kjellsson

/* -*- c++ -*- */
/*
 * Copyright 2005,2006 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 <gr_io_signature.h>
#include <gr_prefs.h>
#include <gr_clock_recovery_mm_cc.h>
#include <gri_mmse_fir_interpolator_cc.h>
#include <stdexcept>
#include <cstdio>


// Public constructor


gr_clock_recovery_mm_cc_sptr 
gr_make_clock_recovery_mm_cc(float omega, float gain_omega, float mu, float 
gain_mu, float omega_relative_limit, int order){
  return gr_clock_recovery_mm_cc_sptr (new gr_clock_recovery_mm_cc (omega, 
                                                                    gain_omega, 
                                                                    mu,
                                                                    gain_mu,
                                                                    
omega_relative_limit,
                                                                                
order));
}

gr_clock_recovery_mm_cc::gr_clock_recovery_mm_cc (float omega, float 
gain_omega, float mu, float gain_mu, float omega_relative_limit, int order)
  : gr_block ("clock_recovery_mm_cc",
              gr_make_io_signature (1, 1, sizeof (gr_complex)),
              gr_make_io_signature (1, 2, sizeof (gr_complex))),
    d_mu (mu), d_omega(omega), d_gain_omega(gain_omega), 
    d_omega_relative_limit(omega_relative_limit), 
    d_gain_mu(gain_mu), d_last_sample(0), d_interp(new 
gri_mmse_fir_interpolator_cc()),
    d_verbose(gr_prefs::singleton()->get_bool("clock_recovery_mm_cc", 
"verbose", false)),
    d_p_2T(0), d_p_1T(0), d_p_0T(0), d_c_2T(0), d_c_1T(0), d_c_0T(0)
{
  if (omega <= 0.0)
    throw std::out_of_range ("clock rate must be > 0");
  if (gain_mu <  0  || gain_omega < 0)
    throw std::out_of_range ("Gains must be non-negative");

        switch(order) {
  case 2:
                d_order = 2;
    d_slicer = &gr_clock_recovery_mm_cc::slicer_0deg;
    break;

  case 4:
                d_order = 4;
    d_slicer = &gr_clock_recovery_mm_cc::slicer_45deg;
    break;

  default: 
    throw std::invalid_argument("order must be 2 or 4");
    break;
  }

  set_omega(omega);                     // also sets min and max omega
  set_relative_rate (1.0 / omega);
  set_history(3);                       // ensure 2 extra input sample is 
available
}

gr_clock_recovery_mm_cc::~gr_clock_recovery_mm_cc ()
{
  delete d_interp;
}

void
gr_clock_recovery_mm_cc::forecast(int noutput_items, gr_vector_int 
&ninput_items_required)
{
  unsigned ninputs = ninput_items_required.size();
  for (unsigned i=0; i < ninputs; i++)
    ninput_items_required[i] =
      (int) ceil((noutput_items * d_omega) + d_interp->ntaps());
}

gr_complex
gr_clock_recovery_mm_cc::slicer_0deg (gr_complex sample)
{
  float real=0, imag=0;

  if(sample.real() > 0)
    real = 1;
  if(sample.imag() > 0)
    imag = 1;
  return gr_complex(real,imag);
}

gr_complex
gr_clock_recovery_mm_cc::slicer_45deg (gr_complex sample)
{
  float real= -1, imag = -1;
  if(sample.real() > 0)
    real=1;
  if(sample.imag() > 0)
    imag = 1;
  return gr_complex(real,imag);
}

/*
  Modified Mueller and Muller clock recovery circuit
  Based:
     G. R. Danesfahani, T.G. Jeans, "Optimisation of modified Mueller and 
Muller 
     algorithm,"  Electronics Letters, Vol. 31, no. 13,  22 June 1995, pp. 1032 
- 1033.
*/

static const int FUDGE = 16;

int
gr_clock_recovery_mm_cc::general_work (int noutput_items,
                                       gr_vector_int &ninput_items,
                                       gr_vector_const_void_star &input_items,
                                       gr_vector_void_star &output_items)
{
  const gr_complex *in = (const gr_complex *) input_items[0];
  gr_complex *out = (gr_complex *) output_items[0];
  gr_complex *foptr = (gr_complex *) output_items[1];

  bool write_foptr = output_items.size() >= 2;
  
  int  ii = 0;                          // input index
  int  oo = 0;                          // output index
  int  ni = ninput_items[0] - d_interp->ntaps() - FUDGE;  // don't use more 
input than this

  assert(d_mu >= 0.0);
  assert(d_mu <= 1.0);

  float mm_val=0;
  gr_complex u, x, y;

  // This loop writes the error to the second output, if it exists
  if (write_foptr) {
    while(oo < noutput_items && ii < ni) {
      d_p_2T = d_p_1T;
      d_p_1T = d_p_0T;
      d_p_0T = d_interp->interpolate (&in[ii], d_mu);

      d_c_2T = d_c_1T;
      d_c_1T = d_c_0T;
      //d_c_0T = slicer_0deg(d_p_0T);
                        d_c_0T = (*this.*d_slicer)(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_val = u.real();
                        /*
                        if(d_order == 4){
                                //mm_val = 
sqrt(u.real()*u.real()+u.imag()*u.imag());
                                mm_val = abs(u);
                        }else{
                                if(d_order == 2){
                                        mm_val = u.real();      
                                }
                        }*/
      out[oo++] = d_p_0T;
      
      // limit mm_val
      mm_val = gr_branchless_clip(mm_val,1.0);
      d_omega = d_omega + d_gain_omega * mm_val;
      d_omega = d_omega_mid + gr_branchless_clip(d_omega-d_omega_mid, 
d_omega_relative_limit);   // make sure we don't walk away

      d_mu = d_mu + d_omega + d_gain_mu * mm_val;
      ii += (int)floor(d_mu);
      d_mu -= floor(d_mu);
            
      // write the error signal to the second output
      foptr[oo-1] = gr_complex(d_mu,0);
      
      if (ii < 0)       // clamp it.  This should only happen with bogus input
        ii = 0;
    }
  }
  // This loop does not write to the second output (ugly, but faster)
  else {
    while(oo < noutput_items && ii < ni) {
      d_p_2T = d_p_1T;
      d_p_1T = d_p_0T;
      d_p_0T = d_interp->interpolate (&in[ii], d_mu);

      d_c_2T = d_c_1T;
      d_c_1T = d_c_0T;
      //d_c_0T = slicer_0deg(d_p_0T);
                        d_c_0T = (*this.*d_slicer)(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_val = u.real();
                        /*
                        if(d_order == 4){
                                mm_val = abs(u);
                        }else{
                                if(d_order == 2){
                                        mm_val = u.real();      
                                }
                        }
                        */
      out[oo++] = d_p_0T;
      
      // limit mm_val
      mm_val = gr_branchless_clip(mm_val,1.0);
      
      d_omega = d_omega + d_gain_omega * mm_val;
      d_omega = d_omega_mid + gr_branchless_clip(d_omega-d_omega_mid, 
d_omega_relative_limit);   // make sure we don't walk away
      
      d_mu = d_mu + d_omega + d_gain_mu * mm_val;
      ii += (int)floor(d_mu);
      d_mu -= floor(d_mu);
      
      if(d_verbose) {
        printf("%f\t%f\n", d_omega, d_mu);
      }
            
      if (ii < 0)       // clamp it.  This should only happen with bogus input
        ii = 0;
    }
  }

  if (ii > 0){
    if (ii > ninput_items[0]){
      fprintf(stderr, "gr_clock_recovery_mm_cc: ii > ninput_items[0] (%d > 
%d)\n",
              ii, ninput_items[0]);
      assert(0);
    }
    consume_each (ii);
  }

  return oo;
}

/* -*- c++ -*- */
/*
 * Copyright 2004 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_GR_CLOCK_RECOVERY_MM_CC_H
#define INCLUDED_GR_CLOCK_RECOVERY_MM_CC_H

#include <gr_block.h>
#include <gr_complex.h>
#include <gr_math.h>

class gri_mmse_fir_interpolator_cc;

class gr_clock_recovery_mm_cc;
typedef boost::shared_ptr<gr_clock_recovery_mm_cc> gr_clock_recovery_mm_cc_sptr;

// public constructor
gr_clock_recovery_mm_cc_sptr 
gr_make_clock_recovery_mm_cc (float omega, float gain_omega, float mu, float 
gain_mu, float omega_relative_limit=0.001, int order = 2);

/*!
 * \brief Mueller and Müller (M&M) based clock recovery block with complex 
input, complex output.
 * \ingroup sync_blk
 *
 * This implements the Mueller and Müller (M&M) discrete-time error-tracking 
synchronizer.
 * The complex version here is based on:
 * Modified Mueller and Muller clock recovery circuit
 * Based:
 *    G. R. Danesfahani, T.G. Jeans, "Optimisation of modified Mueller and 
Muller 
 *    algorithm,"  Electronics Letters, Vol. 31, no. 13,  22 June 1995, pp. 
1032 - 1033.
 */
class gr_clock_recovery_mm_cc : public gr_block
{
 public:
  ~gr_clock_recovery_mm_cc ();
  void forecast(int noutput_items, gr_vector_int &ninput_items_required);
  int general_work (int noutput_items,
                    gr_vector_int &ninput_items,
                    gr_vector_const_void_star &input_items,
                    gr_vector_void_star &output_items);
  float mu() const { return d_mu;}
  float omega() const { return d_omega;}
  float gain_mu() const { return d_gain_mu;}
  float gain_omega() const { return d_gain_omega;}
  void set_verbose (bool verbose) { d_verbose = verbose; }

  void set_gain_mu (float gain_mu) { d_gain_mu = gain_mu; }
  void set_gain_omega (float gain_omega) { d_gain_omega = gain_omega; }
  void set_mu (float mu) { d_mu = mu; }
  void set_omega (float omega) { 
    d_omega = omega;
    d_min_omega = omega*(1.0 - d_omega_relative_limit);
    d_max_omega = omega*(1.0 + d_omega_relative_limit);
    d_omega_mid = 0.5*(d_min_omega+d_max_omega);
  }

protected:
  gr_clock_recovery_mm_cc (float omega, float gain_omega, float mu, float 
gain_mu, float omega_relative_limit, int order);

 private:
  float                         d_mu;
  float                         d_omega;
  float                         d_gain_omega;
  float                         d_min_omega;            // minimum allowed omega
  float                         d_max_omega;            // maximum allowed omeg
  float                         d_omega_relative_limit; // used to compute min 
and max omega
  float                         d_omega_mid;
  float                         d_gain_mu;
  gr_complex    d_last_sample;
  gri_mmse_fir_interpolator_cc  *d_interp;
  bool                          d_verbose;
        int                                     d_order;

  gr_complex                    d_p_2T;
  gr_complex                    d_p_1T;
  gr_complex                    d_p_0T;

  gr_complex                    d_c_2T;
  gr_complex                    d_c_1T;
  gr_complex                    d_c_0T;

  gr_complex slicer_0deg (gr_complex sample);
  gr_complex slicer_45deg (gr_complex sample);

        gr_complex (gr_clock_recovery_mm_cc::*d_slicer)(gr_complex sample);
        //gr_complex (gr_clock_recovery_mm_cc::*d_slicer)(gr_complex sample) 
const;
  friend gr_clock_recovery_mm_cc_sptr
  gr_make_clock_recovery_mm_cc (float omega, float gain_omega, float mu, float 
gain_mu, float omega_relative_limit, int order);
};

#endif