[Commit-gnuradio] r11527 - gnuradio/branches/developers/n4hy/iir_pfb/gnuradio-core/src/lib/filter

[n4hy]

Author: n4hy
Date: 2009-08-02 23:19:46 -0600 (Sun, 02 Aug 2009)
New Revision: 11527

Modified:
   
gnuradio/branches/developers/n4hy/iir_pfb/gnuradio-core/src/lib/filter/gri_sos_iir.h
Log:
Implementing second order sections as Transposed Direct Form II.  This has the 
advantage of applying the zeros first and poles later increasing stability.

Modified: 
gnuradio/branches/developers/n4hy/iir_pfb/gnuradio-core/src/lib/filter/gri_sos_iir.h
===================================================================
--- 
gnuradio/branches/developers/n4hy/iir_pfb/gnuradio-core/src/lib/filter/gri_sos_iir.h
        2009-08-02 15:45:12 UTC (rev 11526)
+++ 
gnuradio/branches/developers/n4hy/iir_pfb/gnuradio-core/src/lib/filter/gri_sos_iir.h
        2009-08-03 05:19:46 UTC (rev 11527)
@@ -30,13 +30,14 @@
 public:
   /*!
    * \brief Construct an SOS_IIR with the given taps.
-   *
-   * This filter uses the Direct Form I implementation, where
+   * 
+   * This filter uses the Transposed Direct Form II implementation, where
    * \p fftaps contains the feed-forward taps, and \p fbtaps the feedback ones.
    *
    * \p fftaps and \p fbtaps must have equal numbers of taps
    *
    * The input and output satisfy a difference equation of the form
+   *  As gri_iir but M == 2.
 
    \f[
    y[n] - \sum_{k=1}^{M} a_k y[n-k] = \sum_{k=0}^{N} b_k x[n-k]
@@ -50,6 +51,8 @@
 
    * Note that some texts define the system function with a + in the 
denominator.
    * If you're using that convention, you'll need to negate the feedback taps.
+
+   * Reference: 
http://ccrma-www.stanford.edu/~jos/fp/Transposed_Direct_Forms.html
    */
   gri_sos_iir (const std::vector<float>& fftaps,
          const std::vector<float>& fbtaps) throw (std::invalid_argument)
@@ -57,7 +60,7 @@
     set_taps (fftaps, fbtaps);
   }
 
-  gri_sos_iir () : d_v2(0),d_v1(0) { }
+  gri_sos_iir () : d_s2(0),d_s1(0) { }
 
   ~gri_sos_iir () {}
 
@@ -91,16 +94,13 @@
 protected:
   std::vector<float>   d_fftaps;
   std::vector<float>   d_fbtaps;
-  int                  d_latest_n;
-  int                  d_latest_m;
-  double               d_v2;
-  double                d_v1;
-  std::vector<float>   d_prev_input;
+  double               d_s2;
+  double                d_s1;
 };
 
 
 //
-// general case.  We may want to specialize this
+//  Transposed Direct Form II SOS 
 //
 float
 gri_sos_iir<float, float, float>::filter (const float input)
@@ -113,26 +113,10 @@
     return (float) 0;
 
   
-  acc = d_fbtaps[0] * input - d_v1;
-  d_v1 = d_fbtaps[1] * input - d_fftaps[1] * 
-  acc -= d_fbtaps[1] * d_v1 + d_fbtaps[2]*d_v2;
-  acc += d_fftaps[0]*d_v0 + d_fftapbs[1]*d_v1 +d_fftaps[2]*d_v2;
+  acc  = (double)d_fftaps[0] * input + d_s1;
+  d_s1 = (double)d_fftaps[1] * input - acc * d_fbtaps[1] + d_s2;
+  d_s2 = (double)d_fftaps[2] * input - acc * d_fbtaps[2];
 
-  // 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 (float) acc;
 }