GNU Radio 3.5.1 C++ API
gri_iir.h
Go to the documentation of this file.
00001 /* -*- c++ -*- */
00002 /*
00003  * Copyright 2002 Free Software Foundation, Inc.
00004  * 
00005  * This file is part of GNU Radio
00006  * 
00007  * GNU Radio is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 3, or (at your option)
00010  * any later version.
00011  * 
00012  * GNU Radio is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  * 
00017  * You should have received a copy of the GNU General Public License
00018  * along with GNU Radio; see the file COPYING.  If not, write to
00019  * the Free Software Foundation, Inc., 51 Franklin Street,
00020  * Boston, MA 02110-1301, USA.
00021  */
00022 
00023 #ifndef INCLUDED_GRI_IIR_H
00024 #define INCLUDED_GRI_IIR_H
00025 
00026 #include <gr_core_api.h>
00027 #include <vector>
00028 #include <stdexcept>
00029 
00030 /*!
00031  * \brief base class template for Infinite Impulse Response filter (IIR)
00032  */
00033 template<class i_type, class o_type, class tap_type> 
00034 class gri_iir {
00035 public:
00036   /*!
00037    * \brief Construct an IIR with the given taps.
00038    *
00039    * This filter uses the Direct Form I implementation, where
00040    * \p fftaps contains the feed-forward taps, and \p fbtaps the feedback ones.
00041    *
00042    * \p fftaps and \p fbtaps must have equal numbers of taps
00043    *
00044    * The input and output satisfy a difference equation of the form
00045 
00046    \f[
00047    y[n] - \sum_{k=1}^{M} a_k y[n-k] = \sum_{k=0}^{N} b_k x[n-k]
00048    \f]
00049 
00050    * with the corresponding rational system function
00051 
00052    \f[
00053    H(z) = \frac{\sum_{k=0}^{N} b_k z^{-k}}{1 - \sum_{k=1}^{M} a_k z^{-k}}
00054    \f]
00055 
00056    * Note that some texts define the system function with a + in the denominator.
00057    * If you're using that convention, you'll need to negate the feedback taps.
00058    */
00059   gri_iir (const std::vector<tap_type>& fftaps,
00060           const std::vector<tap_type>& fbtaps) throw (std::invalid_argument)
00061   {
00062     set_taps (fftaps, fbtaps);
00063   }
00064 
00065   gri_iir () : d_latest_n(0),d_latest_m(0) { }
00066 
00067   ~gri_iir () {}
00068 
00069   /*!
00070    * \brief compute a single output value.
00071    * \returns the filtered input value.
00072    */
00073   o_type filter (const i_type input);
00074 
00075   /*!
00076    * \brief compute an array of N output values.
00077    * \p input must have N valid entries.
00078    */
00079   void filter_n (o_type output[], const i_type input[], long n);
00080 
00081   /*!
00082    * \return number of taps in filter.
00083    */
00084   unsigned ntaps_ff () const { return d_fftaps.size (); }
00085   unsigned ntaps_fb () const { return d_fbtaps.size (); }
00086 
00087   /*!
00088    * \brief install new taps.
00089    */
00090   void set_taps (const std::vector<tap_type> &fftaps, 
00091                  const std::vector<tap_type> &fbtaps) throw (std::invalid_argument)
00092   { 
00093 
00094 
00095     d_latest_n = 0;
00096     d_latest_m = 0;
00097     d_fftaps = fftaps; 
00098     d_fbtaps = fbtaps; 
00099 
00100     int n = fftaps.size ();
00101     int m = fbtaps.size ();
00102     d_prev_input.resize (2 * n);
00103     d_prev_output.resize (2 * m);
00104 
00105     for (int i = 0; i < 2 * n; i++){
00106       d_prev_input[i] = 0;
00107      }
00108     for (int i = 0; i < 2 * m; i++){
00109       d_prev_output[i] = 0;
00110     }
00111   }
00112 
00113 protected:
00114   std::vector<tap_type> d_fftaps;
00115   std::vector<tap_type> d_fbtaps;
00116   int                   d_latest_n;
00117   int                   d_latest_m;
00118   std::vector<tap_type> d_prev_output;
00119   std::vector<i_type>   d_prev_input;
00120 };
00121 
00122 
00123 //
00124 // general case.  We may want to specialize this
00125 //
00126 template<class i_type, class o_type, class tap_type> 
00127 o_type
00128 gri_iir<i_type, o_type, tap_type>::filter (const i_type input)
00129 {
00130   tap_type      acc;
00131   unsigned      i = 0;
00132   unsigned      n = ntaps_ff ();
00133   unsigned      m = ntaps_fb ();
00134 
00135   if (n == 0)
00136     return (o_type) 0;
00137 
00138   int latest_n = d_latest_n;
00139   int latest_m = d_latest_m;
00140   
00141   acc = d_fftaps[0] * input;
00142   for (i = 1; i < n; i ++)
00143     acc += (d_fftaps[i] * d_prev_input[latest_n + i]);
00144   for (i = 1; i < m; i ++)
00145     acc += (d_fbtaps[i] * d_prev_output[latest_m + i]);
00146 
00147   // store the values twice to avoid having to handle wrap-around in the loop
00148   d_prev_output[latest_m] = acc;
00149   d_prev_output[latest_m+m] = acc;
00150   d_prev_input[latest_n] = input;
00151   d_prev_input[latest_n+n] = input;
00152 
00153   latest_n--;
00154   latest_m--;
00155   if (latest_n < 0)
00156     latest_n += n;
00157   if (latest_m < 0)
00158     latest_m += m;
00159 
00160   d_latest_m = latest_m;
00161   d_latest_n = latest_n;
00162   return (o_type) acc;
00163 }
00164 
00165 
00166 template<class i_type, class o_type, class tap_type> 
00167 void 
00168 gri_iir<i_type, o_type, tap_type>::filter_n (o_type output[],
00169                                              const i_type input[],
00170                                              long n)
00171 {
00172   for (int i = 0; i < n; i++)
00173     output[i] = filter (input[i]);
00174 }
00175 
00176 #endif /* INCLUDED_GRI_IIR_H */
00177