3 * Copyright 2009,2010 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.
27 #include <gr_pfb_arb_resampler_ccf.h>
28 #include <gr_fir_ccf.h>
29 #include <gr_fir_util.h>
30 #include <gr_io_signature.h>
33 gr_pfb_arb_resampler_ccf_sptr gr_make_pfb_arb_resampler_ccf (float rate,
34 const std::vector<float> &taps,
35 unsigned int filter_size)
37 return gnuradio::get_initial_sptr(new gr_pfb_arb_resampler_ccf (rate, taps,
42 gr_pfb_arb_resampler_ccf::gr_pfb_arb_resampler_ccf (float rate,
43 const std::vector<float> &taps,
44 unsigned int filter_size)
45 : gr_block ("pfb_arb_resampler_ccf",
46 gr_make_io_signature (1, 1, sizeof(gr_complex)),
47 gr_make_io_signature (1, 1, sizeof(gr_complex))),
50 d_acc = 0; // start accumulator at 0
52 /* The number of filters is specified by the user as the filter size;
53 this is also the interpolation rate of the filter. We use it and the
54 rate provided to determine the decimation rate. This acts as a
55 rational resampler. The flt_rate is calculated as the residual
56 between the integer decimation rate and the real decimation rate and
57 will be used to determine to interpolation point of the resampling
60 d_int_rate = filter_size;
61 d_dec_rate = (unsigned int)floor(d_int_rate/rate);
62 d_flt_rate = (d_int_rate/rate) - d_dec_rate;
64 // Store the last filter between calls to work
69 d_filters = std::vector<gr_fir_ccf*>(d_int_rate);
70 d_diff_filters = std::vector<gr_fir_ccf*>(d_int_rate);
72 // Create an FIR filter for each channel and zero out the taps
73 std::vector<float> vtaps(0, d_int_rate);
74 for(int i = 0; i < d_int_rate; i++) {
75 d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps);
76 d_diff_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps);
79 // Now, actually set the filters' taps
80 std::vector<float> dtaps;
81 create_diff_taps(taps, dtaps);
82 create_taps(taps, d_taps, d_filters);
83 create_taps(dtaps, d_dtaps, d_diff_filters);
85 set_relative_rate(rate);
88 gr_pfb_arb_resampler_ccf::~gr_pfb_arb_resampler_ccf ()
90 for(unsigned int i = 0; i < d_int_rate; i++) {
96 gr_pfb_arb_resampler_ccf::create_taps (const std::vector<float> &newtaps,
97 std::vector< std::vector<float> > &ourtaps,
98 std::vector<gr_fir_ccf*> &ourfilter)
102 unsigned int ntaps = newtaps.size();
103 d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate);
105 // Create d_numchan vectors to store each channel's taps
106 ourtaps.resize(d_int_rate);
108 // Make a vector of the taps plus fill it out with 0's to fill
109 // each polyphase filter with exactly d_taps_per_filter
110 std::vector<float> tmp_taps;
112 while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
113 tmp_taps.push_back(0.0);
116 // Partition the filter
117 for(i = 0; i < d_int_rate; i++) {
118 // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
119 ourtaps[d_int_rate-1-i] = std::vector<float>(d_taps_per_filter, 0);
120 for(j = 0; j < d_taps_per_filter; j++) {
121 ourtaps[d_int_rate - 1 - i][j] = tmp_taps[i + j*d_int_rate];
124 // Build a filter for each channel and add it's taps to it
125 ourfilter[i]->set_taps(ourtaps[d_int_rate-1-i]);
128 // Set the history to ensure enough input items for each filter
129 set_history (d_taps_per_filter + 1);
135 gr_pfb_arb_resampler_ccf::create_diff_taps(const std::vector<float> &newtaps,
136 std::vector<float> &difftaps)
138 // Calculate the differential taps (derivative filter) by taking the difference
139 // between two taps. Duplicate the last one to make both filters the same length.
142 for(unsigned int i = 0; i < newtaps.size()-1; i++) {
143 tap = newtaps[i+1] - newtaps[i];
144 difftaps.push_back(tap);
146 difftaps.push_back(tap);
150 gr_pfb_arb_resampler_ccf::print_taps()
153 for(i = 0; i < d_int_rate; i++) {
154 printf("filter[%d]: [", i);
155 for(j = 0; j < d_taps_per_filter; j++) {
156 printf(" %.4e", d_taps[i][j]);
163 gr_pfb_arb_resampler_ccf::general_work (int noutput_items,
164 gr_vector_int &ninput_items,
165 gr_vector_const_void_star &input_items,
166 gr_vector_void_star &output_items)
168 gr_complex *in = (gr_complex *) input_items[0];
169 gr_complex *out = (gr_complex *) output_items[0];
173 return 0; // history requirements may have changed.
176 int i = 0, j, count = d_start_index;
179 // Restore the last filter position
182 // produce output as long as we can and there are enough input samples
183 while((i < noutput_items) && (count < ninput_items[0]-1)) {
185 // start j by wrapping around mod the number of channels
186 while((j < d_int_rate) && (i < noutput_items)) {
187 // Take the current filter and derivative filter output
188 o0 = d_filters[j]->filter(&in[count]);
189 o1 = d_diff_filters[j]->filter(&in[count]);
191 out[i] = o0 + o1*d_acc; // linearly interpolate between samples
194 // Adjust accumulator and index into filterbank
196 j += d_dec_rate + (int)floor(d_acc);
197 d_acc = fmodf(d_acc, 1.0);
199 if(i < noutput_items) { // keep state for next entry
200 float ss = (int)(j / d_int_rate); // number of items to skip ahead by
201 count += ss; // we have fully consumed another input
202 j = j % d_int_rate; // roll filter around
206 // Store the current filter position and start of next sample
208 d_start_index = std::max(0, count - ninput_items[0]);
210 // consume all we've processed but no more than we can
211 consume_each(std::min(count, ninput_items[0]));