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