~ubuntu-branches/ubuntu/trusty/gnuradio/trusty-updates

« back to all changes in this revision

Viewing changes to gnuradio-core/src/lib/filter/gr_pfb_clock_sync_ccf.cc

  • Committer: Package Import Robot
  • Author(s): A. Maitland Bottoms
  • Date: 2012-02-26 21:26:16 UTC
  • mfrom: (1.1.4)
  • Revision ID: package-import@ubuntu.com-20120226212616-vsfkbi1158xshdql
Tags: 3.5.1-1
* new upstream version, re-packaged from scratch with modern tools
    closes: #642716, #645332, #394849, #616832, #590048, #642580,
    #647018, #557050, #559640, #631863
* CMake build

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* -*- c++ -*- */
 
2
/*
 
3
 * Copyright 2009-2011 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 <cstdio>
 
28
#include <cmath>
 
29
 
 
30
#include <gr_pfb_clock_sync_ccf.h>
 
31
#include <gr_fir_ccf.h>
 
32
#include <gr_fir_util.h>
 
33
#include <gr_io_signature.h>
 
34
#include <gr_math.h>
 
35
 
 
36
gr_pfb_clock_sync_ccf_sptr gr_make_pfb_clock_sync_ccf (double sps, float loop_bw,
 
37
                                                       const std::vector<float> &taps,
 
38
                                                       unsigned int filter_size,
 
39
                                                       float init_phase,
 
40
                                                       float max_rate_deviation,
 
41
                                                       int osps)
 
42
{
 
43
  return gnuradio::get_initial_sptr(new gr_pfb_clock_sync_ccf (sps, loop_bw, taps,
 
44
                                                               filter_size,
 
45
                                                               init_phase,
 
46
                                                               max_rate_deviation,
 
47
                                                               osps));
 
48
}
 
49
 
 
50
static int ios[] = {sizeof(gr_complex), sizeof(float), sizeof(float), sizeof(float)};
 
51
static std::vector<int> iosig(ios, ios+sizeof(ios)/sizeof(int));
 
52
gr_pfb_clock_sync_ccf::gr_pfb_clock_sync_ccf (double sps, float loop_bw,
 
53
                                              const std::vector<float> &taps,
 
54
                                              unsigned int filter_size,
 
55
                                              float init_phase,
 
56
                                              float max_rate_deviation,
 
57
                                              int osps)
 
58
  : gr_block ("pfb_clock_sync_ccf",
 
59
              gr_make_io_signature (1, 1, sizeof(gr_complex)),
 
60
              gr_make_io_signaturev (1, 4, iosig)),
 
61
    d_updated (false), d_nfilters(filter_size),
 
62
    d_max_dev(max_rate_deviation),
 
63
    d_osps(osps), d_error(0), d_out_idx(0)
 
64
{
 
65
  d_nfilters = filter_size;
 
66
  d_sps = floor(sps);
 
67
 
 
68
  // Set the damping factor for a critically damped system
 
69
  d_damping = sqrtf(2.0f)/2.0f;
 
70
  
 
71
  // Set the bandwidth, which will then call update_gains()
 
72
  set_loop_bandwidth(loop_bw);
 
73
  
 
74
  // Store the last filter between calls to work
 
75
  // The accumulator keeps track of overflow to increment the stride correctly.
 
76
  // set it here to the fractional difference based on the initial phaes
 
77
  d_k = init_phase;
 
78
  d_rate = (sps-floor(sps))*(double)d_nfilters;
 
79
  d_rate_i = (int)floor(d_rate);
 
80
  d_rate_f = d_rate - (float)d_rate_i;
 
81
  d_filtnum = (int)floor(d_k);
 
82
 
 
83
  d_filters = std::vector<gr_fir_ccf*>(d_nfilters);
 
84
  d_diff_filters = std::vector<gr_fir_ccf*>(d_nfilters);
 
85
 
 
86
  // Create an FIR filter for each channel and zero out the taps
 
87
  std::vector<float> vtaps(0, d_nfilters);
 
88
  for(int i = 0; i < d_nfilters; i++) {
 
89
    d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps);
 
90
    d_diff_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps);
 
91
  }
 
92
 
 
93
  // Now, actually set the filters' taps
 
94
  std::vector<float> dtaps;
 
95
  create_diff_taps(taps, dtaps);
 
96
  set_taps(taps, d_taps, d_filters);
 
97
  set_taps(dtaps, d_dtaps, d_diff_filters);
 
98
}
 
99
 
 
100
gr_pfb_clock_sync_ccf::~gr_pfb_clock_sync_ccf ()
 
101
{
 
102
  for(int i = 0; i < d_nfilters; i++) {
 
103
    delete d_filters[i];
 
104
    delete d_diff_filters[i];
 
105
  }
 
106
}
 
107
 
 
108
bool
 
109
gr_pfb_clock_sync_ccf::check_topology(int ninputs, int noutputs)
 
110
{
 
111
  return noutputs == 1 || noutputs == 4;
 
112
}
 
113
 
 
114
 
 
115
 
 
116
/*******************************************************************
 
117
    SET FUNCTIONS
 
118
*******************************************************************/
 
119
 
 
120
 
 
121
void
 
122
gr_pfb_clock_sync_ccf::set_loop_bandwidth(float bw) 
 
123
{
 
124
  if(bw < 0) {
 
125
    throw std::out_of_range ("gr_pfb_clock_sync_cc: invalid bandwidth. Must be >= 0.");
 
126
  }
 
127
  
 
128
  d_loop_bw = bw;
 
129
  update_gains();
 
130
}
 
131
 
 
132
void
 
133
gr_pfb_clock_sync_ccf::set_damping_factor(float df) 
 
134
{
 
135
  if(df < 0 || df > 1.0) {
 
136
    throw std::out_of_range ("gr_pfb_clock_sync_cc: invalid damping factor. Must be in [0,1].");
 
137
  }
 
138
  
 
139
  d_damping = df;
 
140
  update_gains();
 
141
}
 
142
 
 
143
void
 
144
gr_pfb_clock_sync_ccf::set_alpha(float alpha)
 
145
{
 
146
  if(alpha < 0 || alpha > 1.0) {
 
147
    throw std::out_of_range ("gr_pfb_clock_sync_cc: invalid alpha. Must be in [0,1].");
 
148
  }
 
149
  d_alpha = alpha;
 
150
}
 
151
 
 
152
void
 
153
gr_pfb_clock_sync_ccf::set_beta(float beta)
 
154
{
 
155
  if(beta < 0 || beta > 1.0) {
 
156
    throw std::out_of_range ("gr_pfb_clock_sync_cc: invalid beta. Must be in [0,1].");
 
157
  }
 
158
  d_beta = beta;
 
159
}
 
160
 
 
161
/*******************************************************************
 
162
    GET FUNCTIONS
 
163
*******************************************************************/
 
164
 
 
165
 
 
166
float
 
167
gr_pfb_clock_sync_ccf::get_loop_bandwidth() const
 
168
{
 
169
  return d_loop_bw;
 
170
}
 
171
 
 
172
float
 
173
gr_pfb_clock_sync_ccf::get_damping_factor() const
 
174
{
 
175
  return d_damping;
 
176
}
 
177
 
 
178
float
 
179
gr_pfb_clock_sync_ccf::get_alpha() const
 
180
{
 
181
  return d_alpha;
 
182
}
 
183
 
 
184
float
 
185
gr_pfb_clock_sync_ccf::get_beta() const
 
186
{
 
187
  return d_beta;
 
188
}
 
189
 
 
190
float
 
191
gr_pfb_clock_sync_ccf::get_clock_rate() const
 
192
{
 
193
  return d_rate_f;
 
194
}
 
195
 
 
196
/*******************************************************************
 
197
*******************************************************************/
 
198
 
 
199
void
 
200
gr_pfb_clock_sync_ccf::update_gains()
 
201
{
 
202
  float denom = (1.0 + 2.0*d_damping*d_loop_bw + d_loop_bw*d_loop_bw);
 
203
  d_alpha = (4*d_damping*d_loop_bw) / denom;
 
204
  d_beta = (4*d_loop_bw*d_loop_bw) / denom;
 
205
}
 
206
 
 
207
 
 
208
void
 
209
gr_pfb_clock_sync_ccf::set_taps (const std::vector<float> &newtaps,
 
210
                                 std::vector< std::vector<float> > &ourtaps,
 
211
                                 std::vector<gr_fir_ccf*> &ourfilter)
 
212
{
 
213
  int i,j;
 
214
 
 
215
  unsigned int ntaps = newtaps.size();
 
216
  d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_nfilters);
 
217
 
 
218
  // Create d_numchan vectors to store each channel's taps
 
219
  ourtaps.resize(d_nfilters);
 
220
  
 
221
  // Make a vector of the taps plus fill it out with 0's to fill
 
222
  // each polyphase filter with exactly d_taps_per_filter
 
223
  std::vector<float> tmp_taps;
 
224
  tmp_taps = newtaps;
 
225
  while((float)(tmp_taps.size()) < d_nfilters*d_taps_per_filter) {
 
226
    tmp_taps.push_back(0.0);
 
227
  }
 
228
  
 
229
  // Partition the filter
 
230
  for(i = 0; i < d_nfilters; i++) {
 
231
    // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
 
232
    //ourtaps[d_nfilters-1-i] = std::vector<float>(d_taps_per_filter, 0);
 
233
    ourtaps[i] = std::vector<float>(d_taps_per_filter, 0);
 
234
    for(j = 0; j < d_taps_per_filter; j++) {
 
235
      //ourtaps[d_nfilters - 1 - i][j] = tmp_taps[i + j*d_nfilters];
 
236
      ourtaps[i][j] = tmp_taps[i + j*d_nfilters];
 
237
    }
 
238
    
 
239
    // Build a filter for each channel and add it's taps to it
 
240
    //ourfilter[i]->set_taps(ourtaps[d_nfilters-1-i]);
 
241
    ourfilter[i]->set_taps(ourtaps[i]);
 
242
  }
 
243
 
 
244
  // Set the history to ensure enough input items for each filter
 
245
  set_history (d_taps_per_filter + d_sps);
 
246
 
 
247
  // Make sure there is enough output space for d_osps outputs/input.
 
248
  set_output_multiple(d_osps);
 
249
 
 
250
  d_updated = true;
 
251
}
 
252
 
 
253
void
 
254
gr_pfb_clock_sync_ccf::create_diff_taps(const std::vector<float> &newtaps,
 
255
                                        std::vector<float> &difftaps)
 
256
{
 
257
  std::vector<float> diff_filter(3);
 
258
  diff_filter[0] = -1;
 
259
  diff_filter[1] = 0;
 
260
  diff_filter[2] = 1;
 
261
 
 
262
  float pwr = 0;
 
263
  difftaps.push_back(0);
 
264
  for(unsigned int i = 0; i < newtaps.size()-2; i++) {
 
265
    float tap = 0;
 
266
    for(int j = 0; j < 3; j++) {
 
267
      tap += diff_filter[j]*newtaps[i+j];
 
268
      pwr += fabsf(tap);
 
269
    }
 
270
    difftaps.push_back(tap);
 
271
  }
 
272
  difftaps.push_back(0);
 
273
 
 
274
  for(unsigned int i = 0; i < difftaps.size(); i++) {
 
275
    difftaps[i] *= pwr;
 
276
  }
 
277
}
 
278
 
 
279
std::string
 
280
gr_pfb_clock_sync_ccf::get_taps_as_string()
 
281
{
 
282
  int i, j;
 
283
  std::stringstream str;
 
284
  str.precision(4);
 
285
  str.setf(std::ios::scientific);
 
286
 
 
287
  str << "[ ";
 
288
  for(i = 0; i < d_nfilters; i++) {
 
289
    str << "[" << d_taps[i][0] << ", ";
 
290
    for(j = 1; j < d_taps_per_filter-1; j++) {
 
291
      str << d_taps[i][j] << ", ";
 
292
    }
 
293
    str << d_taps[i][j] << "],";
 
294
  }
 
295
  str << " ]" << std::endl;
 
296
  
 
297
  return str.str();
 
298
}
 
299
 
 
300
std::string
 
301
gr_pfb_clock_sync_ccf::get_diff_taps_as_string()
 
302
{
 
303
  int i, j;
 
304
  std::stringstream str;
 
305
  str.precision(4);
 
306
  str.setf(std::ios::scientific);
 
307
 
 
308
  str << "[ ";
 
309
  for(i = 0; i < d_nfilters; i++) {
 
310
    str << "[" << d_dtaps[i][0] << ", ";
 
311
    for(j = 1; j < d_taps_per_filter-1; j++) {
 
312
      str << d_dtaps[i][j] << ", ";
 
313
    }
 
314
    str << d_dtaps[i][j] << "],";
 
315
  }
 
316
  str << " ]" << std::endl;
 
317
 
 
318
  return str.str();
 
319
}
 
320
 
 
321
std::vector< std::vector<float> > 
 
322
gr_pfb_clock_sync_ccf::get_taps()
 
323
{
 
324
  return d_taps;
 
325
}
 
326
 
 
327
std::vector< std::vector<float> > 
 
328
gr_pfb_clock_sync_ccf::get_diff_taps()
 
329
{
 
330
  return d_dtaps;
 
331
}
 
332
 
 
333
std::vector<float>
 
334
gr_pfb_clock_sync_ccf::get_channel_taps(int channel)
 
335
{
 
336
  std::vector<float> taps;
 
337
  for(int i = 0; i < d_taps_per_filter; i++) {
 
338
    taps.push_back(d_taps[channel][i]);
 
339
  }
 
340
  return taps;
 
341
}
 
342
 
 
343
std::vector<float>
 
344
gr_pfb_clock_sync_ccf::get_diff_channel_taps(int channel)
 
345
{
 
346
  std::vector<float> taps;
 
347
  for(int i = 0; i < d_taps_per_filter; i++) {
 
348
    taps.push_back(d_dtaps[channel][i]);
 
349
  }
 
350
  return taps;
 
351
}
 
352
 
 
353
 
 
354
int
 
355
gr_pfb_clock_sync_ccf::general_work (int noutput_items,
 
356
                                     gr_vector_int &ninput_items,
 
357
                                     gr_vector_const_void_star &input_items,
 
358
                                     gr_vector_void_star &output_items)
 
359
{
 
360
  gr_complex *in = (gr_complex *) input_items[0];
 
361
  gr_complex *out = (gr_complex *) output_items[0];
 
362
 
 
363
  float *err = NULL, *outrate = NULL, *outk = NULL;
 
364
  if(output_items.size() == 4) {
 
365
    err = (float *) output_items[1];
 
366
    outrate = (float*)output_items[2];
 
367
    outk = (float*)output_items[3];
 
368
  }
 
369
  
 
370
  if (d_updated) {
 
371
    d_updated = false;
 
372
    return 0;                // history requirements may have changed.
 
373
  }
 
374
 
 
375
  // We need this many to process one output
 
376
  int nrequired = ninput_items[0] - d_taps_per_filter - d_osps;
 
377
 
 
378
  int i = 0, count = 0;
 
379
  float error_r, error_i;
 
380
 
 
381
  // produce output as long as we can and there are enough input samples
 
382
  while((i < noutput_items) && (count < nrequired)) {
 
383
    while(d_out_idx < d_osps) {
 
384
      d_filtnum = (int)floor(d_k);
 
385
      
 
386
      // Keep the current filter number in [0, d_nfilters]
 
387
      // If we've run beyond the last filter, wrap around and go to next sample
 
388
      // If we've go below 0, wrap around and go to previous sample
 
389
      while(d_filtnum >= d_nfilters) {
 
390
        d_k -= d_nfilters;
 
391
        d_filtnum -= d_nfilters;
 
392
        count += 1;
 
393
      }
 
394
      while(d_filtnum < 0) {
 
395
        d_k += d_nfilters;
 
396
        d_filtnum += d_nfilters;
 
397
        count -= 1;
 
398
      }
 
399
      
 
400
      out[i+d_out_idx] = d_filters[d_filtnum]->filter(&in[count+d_out_idx]);
 
401
      d_k = d_k + d_rate_i + d_rate_f; // update phase
 
402
      d_out_idx++;
 
403
      
 
404
      if(output_items.size() == 4) {
 
405
        err[i] = d_error;
 
406
        outrate[i] = d_rate_f;
 
407
        outk[i] = d_k;
 
408
      }
 
409
 
 
410
      // We've run out of output items we can create; return now.
 
411
      if(i+d_out_idx >= noutput_items) {
 
412
        consume_each(count);
 
413
        return i;
 
414
      }
 
415
    }
 
416
 
 
417
    // reset here; if we didn't complete a full osps samples last time,
 
418
    // the early return would take care of it.
 
419
    d_out_idx = 0;
 
420
 
 
421
    // Update the phase and rate estimates for this symbol
 
422
    gr_complex diff = d_diff_filters[d_filtnum]->filter(&in[count]);
 
423
    error_r = out[i].real() * diff.real();
 
424
    error_i = out[i].imag() * diff.imag();
 
425
    d_error = (error_i + error_r) / 2.0;       // average error from I&Q channel
 
426
    
 
427
    // Run the control loop to update the current phase (k) and
 
428
    // tracking rate estimates based on the error value
 
429
    d_rate_f = d_rate_f + d_beta*d_error;
 
430
    d_k = d_k + d_alpha*d_error; 
 
431
    
 
432
    // Keep our rate within a good range
 
433
    d_rate_f = gr_branchless_clip(d_rate_f, d_max_dev);
 
434
 
 
435
    i+=d_osps;
 
436
    count += (int)floor(d_sps);
 
437
  }
 
438
 
 
439
  consume_each(count);
 
440
  return i;
 
441
}