Adding Doxygen comments to new FFT filter blocks.
[debian/gnuradio] / gnuradio-core / src / lib / filter / gri_fft_filter_ccc.cc
1 /* -*- c++ -*- */
2 /*
3  * Copyright 2010 Free Software Foundation, Inc.
4  * 
5  * This file is part of GNU Radio
6  * 
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)
10  * any later version.
11  * 
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.
16  * 
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.
21  */
22
23 #ifdef HAVE_CONFIG_H
24 #include "config.h"
25 #endif
26
27 #include <gri_fft_filter_ccc.h>
28 #include <gri_fft.h>
29 #include <assert.h>
30 #include <stdexcept>
31 #include <cstdio>
32
33 gri_fft_filter_ccc::gri_fft_filter_ccc (int decimation, const std::vector<gr_complex> &taps)
34   : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), d_invfft(0)
35 {
36   set_taps(taps);
37 }
38
39 gri_fft_filter_ccc::~gri_fft_filter_ccc ()
40 {
41   delete d_fwdfft;
42   delete d_invfft;
43 }
44
45 #if 0
46 static void 
47 print_vector_complex(const std::string label, const std::vector<gr_complex> &x)
48 {
49   std::cout << label;
50   for (unsigned i = 0; i < x.size(); i++)
51     std::cout << x[i] << " ";
52   std::cout << "\n";
53 }
54 #endif
55
56
57 /*
58  * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps
59  */
60 int
61 gri_fft_filter_ccc::set_taps (const std::vector<gr_complex> &taps)
62 {
63   int i = 0;
64   compute_sizes(taps.size());
65
66   d_tail.resize(tailsize());
67   for (i = 0; i < tailsize(); i++)
68     d_tail[i] = 0;
69
70   gr_complex *in = d_fwdfft->get_inbuf();
71   gr_complex *out = d_fwdfft->get_outbuf();
72
73   float scale = 1.0 / d_fftsize;
74   
75   // Compute forward xform of taps.
76   // Copy taps into first ntaps slots, then pad with zeros
77   for (i = 0; i < d_ntaps; i++)
78     in[i] = taps[i] * scale;
79
80   for (; i < d_fftsize; i++)
81     in[i] = 0;
82
83   d_fwdfft->execute();          // do the xform
84
85   // now copy output to d_xformed_taps
86   for (i = 0; i < d_fftsize; i++)
87     d_xformed_taps[i] = out[i];
88   
89   return d_nsamples;
90 }
91
92 // determine and set d_ntaps, d_nsamples, d_fftsize
93
94 void
95 gri_fft_filter_ccc::compute_sizes(int ntaps)
96 {
97   int old_fftsize = d_fftsize;
98   d_ntaps = ntaps;
99   d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2))));
100   d_nsamples = d_fftsize - d_ntaps + 1;
101
102   if (0)
103     fprintf(stderr, "gri_fft_filter_ccc: ntaps = %d, fftsize = %d, nsamples = %d\n",
104             d_ntaps, d_fftsize, d_nsamples);
105
106   assert(d_fftsize == d_ntaps + d_nsamples -1 );
107
108   if (d_fftsize != old_fftsize){        // compute new plans
109     delete d_fwdfft;
110     delete d_invfft;
111     d_fwdfft = new gri_fft_complex(d_fftsize, true);
112     d_invfft = new gri_fft_complex(d_fftsize, false);
113     d_xformed_taps.resize(d_fftsize);
114   }
115 }
116
117 int
118 gri_fft_filter_ccc::filter (int nitems, const gr_complex *input, gr_complex *output)
119 {
120   int dec_ctr = 0;
121   int j = 0;
122   int ninput_items = nitems * d_decimation;
123
124   for (int i = 0; i < ninput_items; i += d_nsamples){
125     
126     memcpy(d_fwdfft->get_inbuf(), &input[i], d_nsamples * sizeof(gr_complex));
127
128     for (j = d_nsamples; j < d_fftsize; j++)
129       d_fwdfft->get_inbuf()[j] = 0;
130
131     d_fwdfft->execute();        // compute fwd xform
132
133     gr_complex *a = d_fwdfft->get_outbuf();
134     gr_complex *b = &d_xformed_taps[0];
135     gr_complex *c = d_invfft->get_inbuf();
136
137     for (j = 0; j < d_fftsize; j++)     // filter in the freq domain
138       c[j] = a[j] * b[j];
139     
140     d_invfft->execute();        // compute inv xform
141
142     // add in the overlapping tail
143
144     for (j = 0; j < tailsize(); j++)
145       d_invfft->get_outbuf()[j] += d_tail[j];
146
147     // copy nsamples to output
148     j = dec_ctr;
149     while (j < d_nsamples) {
150       *output++ = d_invfft->get_outbuf()[j];
151       j += d_decimation;
152     }
153     dec_ctr = (j - d_nsamples);
154
155     // stash the tail
156     memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples,
157            tailsize() * sizeof(gr_complex));
158   }
159
160   assert(dec_ctr == 0);
161
162   return nitems;
163 }