3
* Copyright 2009-2011 Free Software Foundation, Inc.
5
* This file is part of GNU Radio
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)
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.
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.
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>
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,
40
float max_rate_deviation,
43
return gnuradio::get_initial_sptr(new gr_pfb_clock_sync_ccf (sps, loop_bw, taps,
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,
56
float max_rate_deviation,
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)
65
d_nfilters = filter_size;
68
// Set the damping factor for a critically damped system
69
d_damping = sqrtf(2.0f)/2.0f;
71
// Set the bandwidth, which will then call update_gains()
72
set_loop_bandwidth(loop_bw);
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
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);
83
d_filters = std::vector<gr_fir_ccf*>(d_nfilters);
84
d_diff_filters = std::vector<gr_fir_ccf*>(d_nfilters);
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);
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);
100
gr_pfb_clock_sync_ccf::~gr_pfb_clock_sync_ccf ()
102
for(int i = 0; i < d_nfilters; i++) {
104
delete d_diff_filters[i];
109
gr_pfb_clock_sync_ccf::check_topology(int ninputs, int noutputs)
111
return noutputs == 1 || noutputs == 4;
116
/*******************************************************************
118
*******************************************************************/
122
gr_pfb_clock_sync_ccf::set_loop_bandwidth(float bw)
125
throw std::out_of_range ("gr_pfb_clock_sync_cc: invalid bandwidth. Must be >= 0.");
133
gr_pfb_clock_sync_ccf::set_damping_factor(float df)
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].");
144
gr_pfb_clock_sync_ccf::set_alpha(float alpha)
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].");
153
gr_pfb_clock_sync_ccf::set_beta(float beta)
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].");
161
/*******************************************************************
163
*******************************************************************/
167
gr_pfb_clock_sync_ccf::get_loop_bandwidth() const
173
gr_pfb_clock_sync_ccf::get_damping_factor() const
179
gr_pfb_clock_sync_ccf::get_alpha() const
185
gr_pfb_clock_sync_ccf::get_beta() const
191
gr_pfb_clock_sync_ccf::get_clock_rate() const
196
/*******************************************************************
197
*******************************************************************/
200
gr_pfb_clock_sync_ccf::update_gains()
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;
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)
215
unsigned int ntaps = newtaps.size();
216
d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_nfilters);
218
// Create d_numchan vectors to store each channel's taps
219
ourtaps.resize(d_nfilters);
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;
225
while((float)(tmp_taps.size()) < d_nfilters*d_taps_per_filter) {
226
tmp_taps.push_back(0.0);
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];
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]);
244
// Set the history to ensure enough input items for each filter
245
set_history (d_taps_per_filter + d_sps);
247
// Make sure there is enough output space for d_osps outputs/input.
248
set_output_multiple(d_osps);
254
gr_pfb_clock_sync_ccf::create_diff_taps(const std::vector<float> &newtaps,
255
std::vector<float> &difftaps)
257
std::vector<float> diff_filter(3);
263
difftaps.push_back(0);
264
for(unsigned int i = 0; i < newtaps.size()-2; i++) {
266
for(int j = 0; j < 3; j++) {
267
tap += diff_filter[j]*newtaps[i+j];
270
difftaps.push_back(tap);
272
difftaps.push_back(0);
274
for(unsigned int i = 0; i < difftaps.size(); i++) {
280
gr_pfb_clock_sync_ccf::get_taps_as_string()
283
std::stringstream str;
285
str.setf(std::ios::scientific);
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] << ", ";
293
str << d_taps[i][j] << "],";
295
str << " ]" << std::endl;
301
gr_pfb_clock_sync_ccf::get_diff_taps_as_string()
304
std::stringstream str;
306
str.setf(std::ios::scientific);
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] << ", ";
314
str << d_dtaps[i][j] << "],";
316
str << " ]" << std::endl;
321
std::vector< std::vector<float> >
322
gr_pfb_clock_sync_ccf::get_taps()
327
std::vector< std::vector<float> >
328
gr_pfb_clock_sync_ccf::get_diff_taps()
334
gr_pfb_clock_sync_ccf::get_channel_taps(int channel)
336
std::vector<float> taps;
337
for(int i = 0; i < d_taps_per_filter; i++) {
338
taps.push_back(d_taps[channel][i]);
344
gr_pfb_clock_sync_ccf::get_diff_channel_taps(int channel)
346
std::vector<float> taps;
347
for(int i = 0; i < d_taps_per_filter; i++) {
348
taps.push_back(d_dtaps[channel][i]);
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)
360
gr_complex *in = (gr_complex *) input_items[0];
361
gr_complex *out = (gr_complex *) output_items[0];
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];
372
return 0; // history requirements may have changed.
375
// We need this many to process one output
376
int nrequired = ninput_items[0] - d_taps_per_filter - d_osps;
378
int i = 0, count = 0;
379
float error_r, error_i;
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);
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) {
391
d_filtnum -= d_nfilters;
394
while(d_filtnum < 0) {
396
d_filtnum += d_nfilters;
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
404
if(output_items.size() == 4) {
406
outrate[i] = d_rate_f;
410
// We've run out of output items we can create; return now.
411
if(i+d_out_idx >= noutput_items) {
417
// reset here; if we didn't complete a full osps samples last time,
418
// the early return would take care of it.
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
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;
432
// Keep our rate within a good range
433
d_rate_f = gr_branchless_clip(d_rate_f, d_max_dev);
436
count += (int)floor(d_sps);