git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/filter/gri_iir.h

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002 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 2, 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 #ifndef INCLUDED_GRI_IIR_H
  24 #define INCLUDED_GRI_IIR_H
  25
  26 #include <vector>
  27 #include <stdexcept>
  28
  29 /*!
  30  * \brief base class template for Infinite Impulse Response filter (IIR)
  31  */
  32 template<class i_type, class o_type, class tap_type>
  33 class gri_iir {
  34 public:
  35   /*!
  36    * \brief Construct an IIR with the given taps.
  37    *
  38    * This filter uses the Direct Form I implementation, where
  39    * \p fftaps contains the feed-forward taps, and \p fbtaps the feedback ones.
  40    *
  41    * \p fftaps and \p fbtaps must have equal numbers of taps
  42    *
  43    * The input and output satisfy a difference equation of the form
  44
  45    \f[
  46    y[n] - \sum_{k=1}^{N} a_k y[n-k] = \sum_{k=0}^{M} b_k x[n-k]
  47    \f]
  48
  49    * with the corresponding rational system function
  50
  51    \f[
  52    H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 - \sum_{k=1}^{N} a_k z^{-k}}
  53    \f]
  54
  55    * Note that some texts define the system function with a + in the denominator.
  56    * If you're using that convention, you'll need to negate the feedback taps.
  57    */
  58   gri_iir (const std::vector<tap_type>& fftaps,
  59           const std::vector<tap_type>& fbtaps) throw (std::invalid_argument)
  60   {
  61     set_taps (fftaps, fbtaps);
  62   }
  63
  64   gri_iir () : d_latest(0) { }
  65
  66   ~gri_iir () {}
  67
  68   /*!
  69    * \brief compute a single output value.
  70    * \returns the filtered input value.
  71    */
  72   o_type filter (const i_type input);
  73
  74   /*!
  75    * \brief compute an array of N output values.
  76    * \p input must have N valid entries.
  77    */
  78   void filter_n (o_type output[], const i_type input[], long n);
  79
  80   /*!
  81    * \return number of taps in filter.
  82    */
  83   unsigned ntaps () const { return d_fftaps.size (); }
  84
  85   /*!
  86    * \brief install new taps.
  87    */
  88   void set_taps (const std::vector<tap_type> &fftaps,
  89                  const std::vector<tap_type> &fbtaps) throw (std::invalid_argument)
  90   {
  91     if (fftaps.size () != fbtaps.size ())
  92       throw std::invalid_argument ("gri_iir::set_taps");
  93
  94     d_latest = 0;
  95     d_fftaps = fftaps;
  96     d_fbtaps = fbtaps;
  97
  98     int n = fftaps.size ();
  99     d_prev_input.resize (2 * n);
 100     d_prev_output.resize (2 * n);
 101
 102     for (int i = 0; i < 2 * n; i++){
 103       d_prev_input[i] = 0;
 104       d_prev_output[i] = 0;
 105     }
 106   }
 107
 108 protected:
 109   std::vector<tap_type> d_fftaps;
 110   std::vector<tap_type> d_fbtaps;
 111   int                   d_latest;
 112   std::vector<tap_type> d_prev_output;
 113   std::vector<i_type>   d_prev_input;
 114 };
 115
 116
 117 //
 118 // general case.  We may want to specialize this
 119 //
 120 template<class i_type, class o_type, class tap_type>
 121 o_type
 122 gri_iir<i_type, o_type, tap_type>::filter (const i_type input)
 123 {
 124   tap_type      acc;
 125   unsigned      i = 0;
 126   unsigned      n = ntaps ();
 127
 128   if (n == 0)
 129     return (o_type) 0;
 130
 131   int latest = d_latest;
 132
 133   acc = d_fftaps[0] * input;
 134   for (i = 1; i < n; i ++)
 135     acc += (d_fftaps[i] * d_prev_input[latest + i]
 136             + d_fbtaps[i] * d_prev_output[latest + i]);
 137
 138   // store the values twice to avoid having to handle wrap-around in the loop
 139   d_prev_output[latest] = acc;
 140   d_prev_output[latest+n] = acc;
 141   d_prev_input[latest] = input;
 142   d_prev_input[latest+n] = input;
 143
 144   latest--;
 145   if (latest < 0)
 146     latest += n;
 147
 148   d_latest = latest;
 149   return (o_type) acc;
 150 }
 151
 152
 153 template<class i_type, class o_type, class tap_type>
 154 void
 155 gri_iir<i_type, o_type, tap_type>::filter_n (o_type output[],
 156                                              const i_type input[],
 157                                              long n)
 158 {
 159   for (int i = 0; i < n; i++)
 160     output[i] = filter (input[i]);
 161 }
 162
 163 #endif /* INCLUDED_GRI_IIR_H */
 164