3 * Copyright 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 <gri_fft_filter_fff_generic.h>
33 gri_fft_filter_fff_generic::gri_fft_filter_fff_generic (int decimation,
34 const std::vector<float> &taps)
35 : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), d_invfft(0)
40 gri_fft_filter_fff_generic::~gri_fft_filter_fff_generic ()
47 * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps
50 gri_fft_filter_fff_generic::set_taps (const std::vector<float> &taps)
53 compute_sizes(taps.size());
55 d_tail.resize(tailsize());
56 for (i = 0; i < tailsize(); i++)
59 float *in = d_fwdfft->get_inbuf();
60 gr_complex *out = d_fwdfft->get_outbuf();
62 float scale = 1.0 / d_fftsize;
64 // Compute forward xform of taps.
65 // Copy taps into first ntaps slots, then pad with zeros
66 for (i = 0; i < d_ntaps; i++)
67 in[i] = taps[i] * scale;
69 for (; i < d_fftsize; i++)
72 d_fwdfft->execute(); // do the xform
74 // now copy output to d_xformed_taps
75 for (i = 0; i < d_fftsize/2+1; i++)
76 d_xformed_taps[i] = out[i];
81 // determine and set d_ntaps, d_nsamples, d_fftsize
84 gri_fft_filter_fff_generic::compute_sizes(int ntaps)
86 int old_fftsize = d_fftsize;
88 d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2))));
89 d_nsamples = d_fftsize - d_ntaps + 1;
92 fprintf(stderr, "gri_fft_filter_fff_generic: ntaps = %d, fftsize = %d, nsamples = %d\n",
93 d_ntaps, d_fftsize, d_nsamples);
95 assert(d_fftsize == d_ntaps + d_nsamples -1 );
97 if (d_fftsize != old_fftsize){ // compute new plans
100 d_fwdfft = new gri_fft_real_fwd(d_fftsize);
101 d_invfft = new gri_fft_real_rev(d_fftsize);
102 d_xformed_taps.resize(d_fftsize/2+1);
107 gri_fft_filter_fff_generic::filter (int nitems, const float *input, float *output)
111 int ninput_items = nitems * d_decimation;
113 for (int i = 0; i < ninput_items; i += d_nsamples){
115 memcpy(d_fwdfft->get_inbuf(), &input[i], d_nsamples * sizeof(float));
117 for (j = d_nsamples; j < d_fftsize; j++)
118 d_fwdfft->get_inbuf()[j] = 0;
120 d_fwdfft->execute(); // compute fwd xform
122 gr_complex *a = d_fwdfft->get_outbuf();
123 gr_complex *b = &d_xformed_taps[0];
124 gr_complex *c = d_invfft->get_inbuf();
126 for (j = 0; j < d_fftsize/2+1; j++) { // filter in the freq domain
130 d_invfft->execute(); // compute inv xform
132 // add in the overlapping tail
134 for (j = 0; j < tailsize(); j++)
135 d_invfft->get_outbuf()[j] += d_tail[j];
137 // copy nsamples to output
139 //memcpy(out, d_invfft->get_outbuf(), d_nsamples * sizeof(float));
143 while (j < d_nsamples) {
144 *output++ = d_invfft->get_outbuf()[j];
147 dec_ctr = (j - d_nsamples);
150 memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples,
151 tailsize() * sizeof(float));
154 assert(dec_ctr == 0);