GNU Radio 3.5.1 C++ API
|
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