Updated license from GPL version 2 or later to GPL version 3 or later.
[debian/gnuradio] / gnuradio-core / src / lib / general / gr_clock_recovery_mm_cc.cc
1 /* -*- c++ -*- */
2 /*
3  * Copyright 2005,2006 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_clock_recovery_mm_cc.h>
30 #include <gri_mmse_fir_interpolator_cc.h>
31 #include <stdexcept>
32
33 // Public constructor
34
35
36 gr_clock_recovery_mm_cc_sptr 
37 gr_make_clock_recovery_mm_cc(float omega, float gain_omega, float mu, float gain_mu,
38                              float omega_relative_limit)
39 {
40   return gr_clock_recovery_mm_cc_sptr (new gr_clock_recovery_mm_cc (omega, 
41                                                                     gain_omega, 
42                                                                     mu,
43                                                                     gain_mu,
44                                                                     omega_relative_limit));
45 }
46
47 gr_clock_recovery_mm_cc::gr_clock_recovery_mm_cc (float omega, float gain_omega, float mu, 
48                                                   float gain_mu, float omega_relative_limit)
49   : gr_block ("clock_recovery_mm_cc",
50               gr_make_io_signature (1, 1, sizeof (gr_complex)),
51               gr_make_io_signature (1, 2, sizeof (gr_complex))),
52     d_mu (mu), d_omega(omega), d_gain_omega(gain_omega), 
53     d_omega_relative_limit(omega_relative_limit), 
54     d_gain_mu(gain_mu), d_last_sample(0), d_interp(new gri_mmse_fir_interpolator_cc()),
55     d_verbose(gr_prefs::singleton()->get_bool("clock_recovery_mm_cc", "verbose", false)),
56     d_p_2T(0), d_p_1T(0), d_p_0T(0), d_c_2T(0), d_c_1T(0), d_c_0T(0)
57 {
58   if (omega <= 0.0)
59     throw std::out_of_range ("clock rate must be > 0");
60   if (gain_mu <  0  || gain_omega < 0)
61     throw std::out_of_range ("Gains must be non-negative");
62
63   set_omega(omega);                     // also sets min and max omega
64   set_relative_rate (1.0 / omega);
65   set_history(3);                       // ensure 2 extra input sample is available
66 }
67
68 gr_clock_recovery_mm_cc::~gr_clock_recovery_mm_cc ()
69 {
70   delete d_interp;
71 }
72
73 void
74 gr_clock_recovery_mm_cc::forecast(int noutput_items, gr_vector_int &ninput_items_required)
75 {
76   unsigned ninputs = ninput_items_required.size();
77   for (unsigned i=0; i < ninputs; i++)
78     ninput_items_required[i] =
79       (int) ceil((noutput_items * d_omega) + d_interp->ntaps());
80 }
81
82 gr_complex
83 gr_clock_recovery_mm_cc::slicer_0deg (gr_complex sample)
84 {
85   float real=0, imag=0;
86
87   if(sample.real() > 0)
88     real = 1;
89   if(sample.imag() > 0)
90     imag = 1;
91   return gr_complex(real,imag);
92 }
93
94 gr_complex
95 gr_clock_recovery_mm_cc::slicer_45deg (gr_complex sample)
96 {
97   float real= -1, imag = -1;
98   if(sample.real() > 0)
99     real=1;
100   if(sample.imag() > 0)
101     imag = 1;
102   return gr_complex(real,imag);
103 }
104
105 /*
106   Modified Mueller and Muller clock recovery circuit
107   Based:
108      G. R. Danesfahani, T.G. Jeans, "Optimisation of modified Mueller and Muller 
109      algorithm,"  Electronics Letters, Vol. 31, no. 13,  22 June 1995, pp. 1032 - 1033.
110 */
111
112 static const int FUDGE = 16;
113
114 int
115 gr_clock_recovery_mm_cc::general_work (int noutput_items,
116                                        gr_vector_int &ninput_items,
117                                        gr_vector_const_void_star &input_items,
118                                        gr_vector_void_star &output_items)
119 {
120   const gr_complex *in = (const gr_complex *) input_items[0];
121   gr_complex *out = (gr_complex *) output_items[0];
122   gr_complex *foptr = (gr_complex *) output_items[1];
123
124   bool write_foptr = output_items.size() >= 2;
125   
126   int  ii = 0;                          // input index
127   int  oo = 0;                          // output index
128   int  ni = ninput_items[0] - d_interp->ntaps() - FUDGE;  // don't use more input than this
129
130   assert(d_mu >= 0.0);
131   assert(d_mu <= 1.0);
132
133   float mm_val=0;
134   gr_complex u, x, y;
135
136   // This loop writes the error to the second output, if it exists
137   if (write_foptr) {
138     while(oo < noutput_items && ii < ni) {
139       d_p_2T = d_p_1T;
140       d_p_1T = d_p_0T;
141       d_p_0T = d_interp->interpolate (&in[ii], d_mu);
142
143       d_c_2T = d_c_1T;
144       d_c_1T = d_c_0T;
145       d_c_0T = slicer_0deg(d_p_0T);
146       
147       x = (d_c_0T - d_c_2T) * conj(d_p_1T);
148       y = (d_p_0T - d_p_2T) * conj(d_c_1T);
149       u = y - x;
150       mm_val = u.real();
151       out[oo++] = d_p_0T;
152       
153       // limit mm_val
154       if (mm_val > 1.0)
155         mm_val = 1.0;
156       else if (mm_val < -1.0)
157         mm_val = -1.0;
158       
159       d_omega = d_omega + d_gain_omega * mm_val;
160       if (d_omega > d_max_omega)
161         d_omega = d_max_omega;
162       else if (d_omega < d_min_omega)
163         d_omega = d_min_omega;
164       
165       d_mu = d_mu + d_omega + d_gain_mu * mm_val;
166       ii += (int)floor(d_mu);
167       d_mu -= floor(d_mu);
168       
169       #if 0
170       printf("%f\t%f\n", d_omega, d_mu);
171       #endif
172       
173       // write the error signal to the second output
174       foptr[oo-1] = gr_complex(d_mu,0);
175       
176       if (ii < 0)       // clamp it.  This should only happen with bogus input
177         ii = 0;
178     }
179   }
180   // This loop does not write to the second output (ugly, but faster)
181   else {
182     while(oo < noutput_items && ii < ni) {
183       d_p_2T = d_p_1T;
184       d_p_1T = d_p_0T;
185       d_p_0T = d_interp->interpolate (&in[ii], d_mu);
186
187       d_c_2T = d_c_1T;
188       d_c_1T = d_c_0T;
189       d_c_0T = slicer_0deg(d_p_0T);
190       
191       x = (d_c_0T - d_c_2T) * conj(d_p_1T);
192       y = (d_p_0T - d_p_2T) * conj(d_c_1T);
193       u = y - x;
194       mm_val = u.real();
195       out[oo++] = d_p_0T;
196       
197       // limit mm_val
198       if (mm_val > 1.0)
199         mm_val = 1.0;
200       else if (mm_val < -1.0)
201         mm_val = -1.0;
202       
203       d_omega = d_omega + d_gain_omega * mm_val;
204       if (d_omega > d_max_omega)
205         d_omega = d_max_omega;
206       else if (d_omega < d_min_omega)
207         d_omega = d_min_omega;
208       
209       d_mu = d_mu + d_omega + d_gain_mu * mm_val;
210       ii += (int)floor(d_mu);
211       d_mu -= floor(d_mu);
212       
213       if(d_verbose) {
214         printf("%f\t%f\n", d_omega, d_mu);
215       }
216             
217       if (ii < 0)       // clamp it.  This should only happen with bogus input
218         ii = 0;
219     }
220   }
221
222   if (ii > 0){
223     if (ii > ninput_items[0]){
224       fprintf(stderr, "gr_clock_recovery_mm_cc: ii > ninput_items[0] (%d > %d)\n",
225               ii, ninput_items[0]);
226       assert(0);
227     }
228     consume_each (ii);
229   }
230
231   return oo;
232 }