iir_filter.h

Go to the documentation of this file.

/* -*- c++ -*- */
/*
 * Copyright 2002,2012 Free Software Foundation, Inc.
 *
 * This file is part of GNU Radio
 *
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3, or (at your option)
 * any later version.
 *
 * GNU Radio is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifndef INCLUDED_IIR_FILTER_H
#define INCLUDED_IIR_FILTER_H

#include <gnuradio/filter/api.h>
#include <gnuradio/gr_complex.h>
#include <vector>
#include <stdexcept>

namespace gr {
  namespace filter {
    namespace kernel {

      /*!
       * \brief Base class template for Infinite Impulse Response filter (IIR)
       *
       * \details
       *
       * This class provides a templated kernel for IIR filters. These
       * iir_filters can be instantiated with a set of feed-forward
       * and feed-back taps in the constructor. We then call the
       * iir_filter::filter function to add a new sample to the
       * filter, or iir_filter::filter_n to add a vector of samples to
       * be filtered.
       *
       * Instantiating a filter means defining the templates for the
       * data types being processed by the filter. There are four templates:
       *
       * \li i_type the data type of the input data (i.e., float).
       * \li o_type the data type of the output data (i.e., float).
       * \li tap_type the data type of the filter taps (i.e., double).
       * \li acc_type the data type of the internal accumulator (i.e., double).
       *
       * The acc_type is specified to control how data is handled
       * internally in the filter. This should always be the highest
       * precision data type of any of the first three. Often, IIR
       * filters require double-precision values in the taps for
       * stability, and so the internal accumulator should also be
       * double precision.
       *
       * Example:
       *
       * \code
       * gr::filter::kernel::iir_filter<float,float,double,double> iir_filt(fftaps, fbtaps);
       * ...
       * float y = iir_filt.filter(x);
       *
       * <or>
       *
       * iir_filt.filter(y, x, N); // y and x are float arrays
       * \endcode
       *
       * Another example for handling complex samples with
       * double-precision taps (see filter::iir_filter_ccz):
       *
       * \code
       * gr:;filter::kernel::iir_filter<gr_complex, gr_complex, gr_complexd, gr_complexd> iir_filt(fftaps, fbtaps);
       * \endcode
       */
      template<class i_type, class o_type, class tap_type, class acc_type>
      class FILTER_API iir_filter
      {
      public:
        /*!
         * \brief Construct an IIR with the given taps.
         *
         * This filter uses the Direct Form I implementation, where
         * \p fftaps contains the feed-forward taps, and \p fbtaps the feedback ones.
         *
         * \p oldstyle: The old style of the IIR filter uses feedback
         * taps that are negative of what most definitions use (scipy
         * and Matlab among them). This parameter keeps using the old
         * GNU Radio style and is set to TRUE by default. When taps
         * generated from scipy, Matlab, or gr_filter_design, use the
         * new style by setting this to FALSE.
         *
         * When \p oldstyle is set FALSE, the input and output satisfy a
         * difference equation of the form

         \f[
         y[n] + \sum_{k=1}^{M} a_k y[n-k] = \sum_{k=0}^{N} b_k x[n-k]
         \f]

         * with the corresponding rational system function

         \f[
         H(z) = \frac{\sum_{k=0}^{N} b_k z^{-k}}{1 + \sum_{k=1}^{M} a_k z^{-k}}
         \f]
         * where:
         * \f$x\f$ - input signal,
         * \f$y\f$ - output signal,
         * \f$a_k\f$ - k-th feedback tap,
         * \f$b_k\f$ - k-th feed-forward tap,
         * \f$M\f$ - \p len(fbtaps)-1,
         * \f$N\f$ - \p len(fftaps)-1.

         *  \f$a_0\f$, that is \p fbtaps[0], is ignored.
         */
        iir_filter(const std::vector<tap_type>& fftaps,
                   const std::vector<tap_type>& fbtaps,
                   bool oldstyle=true)
          throw (std::invalid_argument)
        {
          d_oldstyle = oldstyle;
          set_taps(fftaps, fbtaps);
        }

        iir_filter() : d_latest_n(0),d_latest_m(0) { }

        ~iir_filter() {}

        /*!
         * \brief compute a single output value.
         * \returns the filtered input value.
         */
        o_type filter(const i_type input);

        /*!
         * \brief compute an array of N output values.
         * \p input must have N valid entries.
         */
        void filter_n(o_type output[], const i_type input[], long n);

        /*!
         * \return number of taps in filter.
         */
        unsigned ntaps_ff() const { return d_fftaps.size(); }
        unsigned ntaps_fb() const { return d_fbtaps.size(); }

        /*!
         * \brief install new taps.
         */
        void set_taps(const std::vector<tap_type> &fftaps,
                      const std::vector<tap_type> &fbtaps)
          throw (std::invalid_argument)
        {
          d_latest_n = 0;
          d_latest_m = 0;
          d_fftaps = fftaps;

          if(d_oldstyle) {
            d_fbtaps = fbtaps;
          }
          else {
            // New style negates taps a[1:N-1] to fit with most IIR
            // tap generating programs.
            d_fbtaps.resize(fbtaps.size());
            d_fbtaps[0] = fbtaps[0];
            for(size_t i = 1; i < fbtaps.size(); i++) {
              d_fbtaps[i] = -fbtaps[i];
            }
          }

          int n = fftaps.size();
          int m = fbtaps.size();
          d_prev_input.clear();
          d_prev_output.clear();
          d_prev_input.resize(2 * n, 0);
          d_prev_output.resize(2 * m, 0);
        }

      protected:
        bool                    d_oldstyle;
        std::vector<tap_type>   d_fftaps;
        std::vector<tap_type>   d_fbtaps;
        int                     d_latest_n;
        int                     d_latest_m;
        std::vector<acc_type>   d_prev_output;
        std::vector<i_type>     d_prev_input;
      };

      //
      // general case.  We may want to specialize this
      //
      template<class i_type, class o_type, class tap_type, class acc_type>
      o_type
      iir_filter<i_type, o_type, tap_type, acc_type>::filter(const i_type input)
      {
        acc_type acc;
        unsigned i = 0;
        unsigned n = ntaps_ff();
        unsigned m = ntaps_fb();

        if(n == 0)
          return (o_type)0;

        int latest_n = d_latest_n;
        int latest_m = d_latest_m;

        acc = d_fftaps[0] * input;
        for(i = 1; i < n; i ++)
          acc += (d_fftaps[i] * d_prev_input[latest_n + i]);
        for(i = 1; i < m; i ++)
          acc += (d_fbtaps[i] * d_prev_output[latest_m + i]);

        // store the values twice to avoid having to handle wrap-around in the loop
        d_prev_output[latest_m] = acc;
        d_prev_output[latest_m+m] = acc;
        d_prev_input[latest_n] = input;
        d_prev_input[latest_n+n] = input;

        latest_n--;
        latest_m--;
        if(latest_n < 0)
          latest_n += n;
        if(latest_m < 0)
          latest_m += m;

        d_latest_m = latest_m;
        d_latest_n = latest_n;
        return (o_type)acc;
      }

      template<class i_type, class o_type, class tap_type, class acc_type>
      void
      iir_filter<i_type, o_type, tap_type, acc_type>::filter_n(o_type output[],
                                                     const i_type input[],
                                                     long n)
      {
        for(int i = 0; i < n; i++)
          output[i] = filter(input[i]);
      }

      template<>
      gr_complex
      iir_filter<gr_complex, gr_complex, float, gr_complex>::filter(const gr_complex input);

      template<>
      gr_complex
      iir_filter<gr_complex, gr_complex, double, gr_complexd>::filter(const gr_complex input);

      template<>
      gr_complex
      iir_filter<gr_complex, gr_complex, gr_complexd, gr_complexd>::filter(const gr_complex input);

    } /* namespace kernel */
  } /* namespace filter */
} /* namespace gr */

#endif /* INCLUDED_IIR_FILTER_H */