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 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 set_taps(taps, d_taps, d_filters);
83 set_taps(dtaps, d_dtaps, d_diff_filters);
86 gr_pfb_arb_resampler_ccf::~gr_pfb_arb_resampler_ccf ()
88 for(unsigned int i = 0; i < d_int_rate; i++) {
94 gr_pfb_arb_resampler_ccf::set_taps (const std::vector<float> &newtaps,
95 std::vector< std::vector<float> > &ourtaps,
96 std::vector<gr_fir_ccf*> &ourfilter)
100 unsigned int ntaps = newtaps.size();
101 d_taps_per_filter = (unsigned int)ceil((double)ntaps/(double)d_int_rate);
103 // Create d_numchan vectors to store each channel's taps
104 ourtaps.resize(d_int_rate);
106 // Make a vector of the taps plus fill it out with 0's to fill
107 // each polyphase filter with exactly d_taps_per_filter
108 std::vector<float> tmp_taps;
110 while((float)(tmp_taps.size()) < d_int_rate*d_taps_per_filter) {
111 tmp_taps.push_back(0.0);
114 // Partition the filter
115 for(i = 0; i < d_int_rate; i++) {
116 // Each channel uses all d_taps_per_filter with 0's if not enough taps to fill out
117 ourtaps[d_int_rate-1-i] = std::vector<float>(d_taps_per_filter, 0);
118 for(j = 0; j < d_taps_per_filter; j++) {
119 ourtaps[d_int_rate - 1 - i][j] = tmp_taps[i + j*d_int_rate];
122 // Build a filter for each channel and add it's taps to it
123 ourfilter[i]->set_taps(ourtaps[d_int_rate-1-i]);
126 // Set the history to ensure enough input items for each filter
127 set_history (d_taps_per_filter + 1);
133 gr_pfb_arb_resampler_ccf::create_diff_taps(const std::vector<float> &newtaps,
134 std::vector<float> &difftaps)
136 // Calculate the differential taps (derivative filter) by taking the difference
137 // between two taps. Duplicate the last one to make both filters the same length.
140 for(unsigned int i = 0; i < newtaps.size()-1; i++) {
141 tap = newtaps[i+1] - newtaps[i];
142 difftaps.push_back(tap);
144 difftaps.push_back(tap);
148 gr_pfb_arb_resampler_ccf::print_taps()
151 for(i = 0; i < d_int_rate; i++) {
152 printf("filter[%d]: [", i);
153 for(j = 0; j < d_taps_per_filter; j++) {
154 printf(" %.4e", d_taps[i][j]);
161 gr_pfb_arb_resampler_ccf::general_work (int noutput_items,
162 gr_vector_int &ninput_items,
163 gr_vector_const_void_star &input_items,
164 gr_vector_void_star &output_items)
166 gr_complex *in = (gr_complex *) input_items[0];
167 gr_complex *out = (gr_complex *) output_items[0];
171 return 0; // history requirements may have changed.
174 int i = 0, j, count = d_start_index;
177 // Restore the last filter position
180 // produce output as long as we can and there are enough input samples
181 while((i < noutput_items) && (count < ninput_items[0]-1)) {
183 // start j by wrapping around mod the number of channels
184 while((j < d_int_rate) && (i < noutput_items)) {
185 // Take the current filter and derivative filter output
186 o0 = d_filters[j]->filter(&in[count]);
187 o1 = d_diff_filters[j]->filter(&in[count]);
189 out[i] = o0 + o1*d_acc; // linearly interpolate between samples
192 // Adjust accumulator and index into filterbank
194 j += d_dec_rate + (int)floor(d_acc);
195 d_acc = fmodf(d_acc, 1.0);
197 if(i < noutput_items) { // keep state for next entry
198 float ss = (int)(j / d_int_rate); // number of items to skip ahead by
199 count += ss; // we have fully consumed another input
200 j = j % d_int_rate; // roll filter around
204 // Store the current filter position and start of next sample
206 d_start_index = std::max(0, count - ninput_items[0]);
208 // consume all we've processed but no more than we can
209 consume_each(std::min(count, ninput_items[0]));