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 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 #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}^{M} a_k y[n-k] = \sum_{k=0}^{N} 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}^{N} b_k z^{-k}}{1 - \sum_{k=1}^{M} 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_n(0),d_latest_m(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_ff () const { return d_fftaps.size (); }
  84   unsigned ntaps_fb () const { return d_fbtaps.size (); }
  85
  86   /*!
  87    * \brief install new taps.
  88    */
  89   void set_taps (const std::vector<tap_type> &fftaps,
  90                  const std::vector<tap_type> &fbtaps) throw (std::invalid_argument)
  91   {
  92
  93
  94     d_latest_n = 0;
  95     d_latest_m = 0;
  96     d_fftaps = fftaps;
  97     d_fbtaps = fbtaps;
  98
  99     int n = fftaps.size ();
 100     int m = fbtaps.size ();
 101     d_prev_input.resize (2 * n);
 102     d_prev_output.resize (2 * m);
 103
 104     for (int i = 0; i < 2 * n; i++){
 105       d_prev_input[i] = 0;
 106      }
 107     for (int i = 0; i < 2 * m; i++){
 108       d_prev_output[i] = 0;
 109     }
 110   }
 111
 112 protected:
 113   std::vector<tap_type> d_fftaps;
 114   std::vector<tap_type> d_fbtaps;
 115   int                   d_latest_n;
 116   int                   d_latest_m;
 117   std::vector<tap_type> d_prev_output;
 118   std::vector<i_type>   d_prev_input;
 119 };
 120
 121
 122 //
 123 // general case.  We may want to specialize this
 124 //
 125 template<class i_type, class o_type, class tap_type>
 126 o_type
 127 gri_iir<i_type, o_type, tap_type>::filter (const i_type input)
 128 {
 129   tap_type      acc;
 130   unsigned      i = 0;
 131   unsigned      n = ntaps_ff ();
 132   unsigned      m = ntaps_fb ();
 133
 134   if (n == 0)
 135     return (o_type) 0;
 136
 137   int latest_n = d_latest_n;
 138   int latest_m = d_latest_m;
 139
 140   acc = d_fftaps[0] * input;
 141   for (i = 1; i < n; i ++)
 142     acc += (d_fftaps[i] * d_prev_input[latest_n + i]);
 143   for (i = 1; i < m; i ++)
 144     acc += (d_fbtaps[i] * d_prev_output[latest_m + i]);
 145
 146   // store the values twice to avoid having to handle wrap-around in the loop
 147   d_prev_output[latest_m] = acc;
 148   d_prev_output[latest_m+m] = acc;
 149   d_prev_input[latest_n] = input;
 150   d_prev_input[latest_n+n] = input;
 151
 152   latest_n--;
 153   latest_m--;
 154   if (latest_n < 0)
 155     latest_n += n;
 156   if (latest_m < 0)
 157     latest_m += m;
 158
 159   d_latest_m = latest_m;
 160   d_latest_n = latest_n;
 161   return (o_type) acc;
 162 }
 163
 164
 165 template<class i_type, class o_type, class tap_type>
 166 void
 167 gri_iir<i_type, o_type, tap_type>::filter_n (o_type output[],
 168                                              const i_type input[],
 169                                              long n)
 170 {
 171   for (int i = 0; i < n; i++)
 172     output[i] = filter (input[i]);
 173 }
 174
 175 #endif /* INCLUDED_GRI_IIR_H */
 176