IIR 滤波器的实现（C++）

最近在写的一个程序需要用到IIR滤波器，而且IIR滤波器的系数需要动态调整。因此就花了点时间研究IIR 滤波器的实现。 以前用到的IIR滤波器的参数都是事先确定好的，有个网站，只要把滤波器的参数特性输进去，直接就能生成需要的C代码。
http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html 一直都偷懒直接用这个网站的结果，所以手上也没积累什么有用的代码。这次就需要自己从头做起。 我面临的问题有两个： 1. 根据滤波器的参数（滤波器的类型、截止频率、滤波器的阶数等），计算出滤波器对应的差分方程的系数。 2. 利用得到的差分方程的系数构造一个可以工作的滤波器。 其中第一个问题，对于不同类型的滤波器，比如Butterworth型、Bessel型等，滤波器系数的计算方法都不同。这部分工作我还没做完全，等我把常见的几种滤波器类型的系数计算方法都实现后再来写一篇文章。 这里就先写写第二个问题。IIR 滤波器对应的差分方程为：
相应的系统函数为： 这里默认a[0] = 1。实际上，总可以通过调整a[k] 与 b[k] 的值使得a[0] = 1，所以这个条件时总能满足的。 按照奥本海默写的《离散时间信号处理》上面的介绍，IIR 滤波器有两种基本的实现形式，分别成为直接I型和直接II型。我分别写了两个类，实现这两种形式。

直接I型

[cpp] view plain copy

1. class IIR_I  
2. {  
3. private:  
4.     double *m_pNum;  
5.     double *m_pDen;  
6.     double *m_px;  
7.     double *m_py;  
8.     int m_num_order;  
9.     int m_den_order;  
10. public:  
11.     IIR_I();  
12.     void reset();  
13.     void setPara(double num[], int num_order, double den[], int den_order);  
14.     void resp(double data_in[], int m, double data_out[], int n);  
15.     double filter(double data);  
16.     void filter(double data[], int len);  
17.     void filter(double data_in[], double data_out[], int len);  
18. };  

其中 m_px 存放x[n-k] 的值（m_px[0]存放x[n-0]、 m_px[1] 存放x[n-1]，以此类推），m_py存放y[n-k] 的值（m_py[0]存放y[n-0]、 m_py[1] 存放y[n-1]，以此类推）。 三个filter函数用来做实际的滤波操作。在这之前，需要用setPara函数初始化滤波器的系数。 下面是实现代码： [cpp] view plain copy

1. /** rief 将滤波器的内部状态清零，滤波器的系数保留 
2.  *  eturn 
3.  */  
4. void IIR_I::reset()  
5. {  
6.     for(int i = 0; i <= m_num_order; i++)  
7.     {  
8.         m_pNum[i] = 0.0;  
9.     }  
10.     for(int i = 0; i <= m_den_order; i++)  
11.     {  
12.         m_pDen[i] = 0.0;  
13.     }  
14. }  
15. IIR_I::IIR_I()  
16. {  
17.     m_pNum = NULL;  
18.     m_pDen = NULL;  
19.     m_px = NULL;  
20.     m_py = NULL;  
21.     m_num_order = -1;  
22.     m_den_order = -1;  
23. };  
24. /** rief 
25.  * 
26.  * param num 分子多项式的系数，升序排列,num[0] 为常数项 
27.  * param m 分子多项式的阶数 
28.  * param den 分母多项式的系数，升序排列,den[0] 为常数项 
29.  * param m 分母多项式的阶数 
30.  *  eturn 
31.  */  
32. void IIR_I::setPara(double num[], int num_order, double den[], int den_order)  
33. {  
34.     delete[] m_pNum;  
35.     delete[] m_pDen;  
36.     delete[] m_px;  
37.     delete[] m_py;  
38.     m_pNum = new double[num_order + 1];  
39.     m_pDen = new double[den_order + 1];  
40.     m_num_order = num_order;  
41.     m_den_order = den_order;  
42.     m_px = new double[num_order + 1];  
43.     m_py = new double[den_order + 1];  
44.     for(int i = 0; i <= m_num_order; i++)  
45.     {  
46.         m_pNum[i] = num[i];  
47.         m_px[i] = 0.0;  
48.     }  
49.     for(int i = 0; i <= m_den_order; i++)  
50.     {  
51.         m_pDen[i] = den[i];  
52.         m_py[i] = 0.0;  
53.     }  
54. }  
55.   
56. /** rief 滤波函数，采用直接I型结构 
57.  * 
58.  * param data 传入输入数据 
59.  *  eturn 滤波后的结果 
60.  */  
61. double IIR_I::filter(double data)  
62. {  
63.     m_py[0] = 0.0; // 存放滤波后的结果  
64.     m_px[0] = data;  
65.     for(int i = 0; i <= m_num_order; i++)  
66.     {  
67.         m_py[0] = m_py[0] + m_pNum[i] * m_px[i];  
68.     }  
69.     for(int i = 1; i <= m_den_order; i++)  
70.     {  
71.         m_py[0] = m_py[0] - m_pDen[i] * m_py[i];  
72.     }  
73.     for(int i = m_num_order; i >= 1; i--)  
74.     {  
75.         m_px[i] = m_px[i-1];  
76.     }  
77.     for(int i = m_den_order; i >= 1; i--)  
78.     {  
79.         m_py[i] = m_py[i-1];