3 * Copyright 2009 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 gr_pfb_arb_resampler_ccf_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 /* The number of filters is specified by the user as the filter size;
51 this is also the interpolation rate of the filter. We use it and the
52 rate provided to determine the decimation rate. This acts as a
53 rational resampler. The flt_rate is calculated as the residual
54 between the integer decimation rate and the real decimation rate and
55 will be used to determine to interpolation point of the resampling
58 d_int_rate = filter_size;
59 d_dec_rate = (unsigned int)floor(d_int_rate/rate);
60 d_flt_rate = (d_int_rate/rate) - d_dec_rate;
62 // The accumulator keeps track of overflow to increment the stride correctly.
65 // Store the last filter between calls to work
68 d_filters = std::vector<gr_fir_ccf*>(d_int_rate);
70 // Create an FIR filter for each channel and zero out the taps
71 std::vector<float> vtaps(0, d_int_rate);
72 for(unsigned int i = 0; i < d_int_rate; i++) {
73 d_filters[i] = gr_fir_util::create_gr_fir_ccf(vtaps);
76 // Now, actually set the filters' taps
80 gr_pfb_arb_resampler_ccf::~gr_pfb_arb_resampler_ccf ()
82 for(unsigned int i = 0; i < d_int_rate; i++) {
88 gr_pfb_arb_resampler_ccf::set_taps (const std::vector<float> &taps)
92 unsigned int ntaps = taps.size();
93 d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate);
95 // Create d_numchan vectors to store each channel's taps
96 d_taps.resize(d_int_rate);
98 // Make a vector of the taps plus fill it out with 0's to fill
99 // each polyphase filter with exactly d_taps_per_filter
100 std::vector<float> tmp_taps;
102 while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
103 tmp_taps.push_back(0.0);
106 // Partition the filter
107 for(i = 0; i < d_int_rate; i++) {
108 // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
109 d_taps[i] = std::vector<float>(d_taps_per_filter, 0);
110 for(j = 0; j < d_taps_per_filter; j++) {
111 d_taps[i][j] = tmp_taps[i + j*d_int_rate]; // add taps to channels in reverse order
114 // Build a filter for each channel and add it's taps to it
115 d_filters[i]->set_taps(d_taps[i]);
118 // Set the history to ensure enough input items for each filter
119 set_history (d_taps_per_filter);
125 gr_pfb_arb_resampler_ccf::print_taps()
128 for(i = 0; i < d_int_rate; i++) {
129 printf("filter[%d]: [", i);
130 for(j = 0; j < d_taps_per_filter; j++) {
131 printf(" %.4e", d_taps[i][j]);
138 gr_pfb_arb_resampler_ccf::general_work (int noutput_items,
139 gr_vector_int &ninput_items,
140 gr_vector_const_void_star &input_items,
141 gr_vector_void_star &output_items)
143 gr_complex *in = (gr_complex *) input_items[0];
144 gr_complex *out = (gr_complex *) output_items[0];
148 return 0; // history requirements may have changed.
151 int i = 0, j, count = 0;
154 // Restore the last filter position
157 // produce output as long as we can and there are enough input samples
158 while((i < noutput_items) && (count < ninput_items[0]-1)) {
160 // start j by wrapping around mod the number of channels
161 while((j < d_int_rate) && (i < noutput_items)) {
162 // Take the current filter output
163 o0 = d_filters[j]->filter(&in[count]);
165 // Take the next filter output; wrap around to 0 if necessary
166 if(j+1 == d_int_rate)
167 // Use the sample of the next input item through the first filter
168 o1 = d_filters[0]->filter(&in[count+1]);
170 // Use the sample from the current input item through the nex filter
171 o1 = d_filters[j+1]->filter(&in[count]);
174 //out[i] = o0; // nearest-neighbor approach
175 out[i] = o0 + (o1 - o0)*d_acc; // linearly interpolate between samples
178 // Accumulate the position in the stream for the interpolated point.
179 // If it goes above 1, roll around to zero and increment the stride
180 // length this time by the decimation rate plus 1 for the increment
181 // due to the acculated position.
183 j += d_dec_rate + (int)floor(d_acc);
184 d_acc = fmodf(d_acc, 1.0);
186 if(i < noutput_items) { // keep state for next entry
187 float ss = (int)(j / d_int_rate); // number of items to skip ahead by
188 count += ss; // we have fully consumed another input
189 j = j % d_int_rate; // roll filter around
193 // Store the current filter position