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