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

时间  2019-11-13
标签 iir 滤波器 实现 c++ 转载 栏目 C&C++ 繁體版
原文   原文链接

转载自 http://blog.csdn.net/liyuanbhu/article/details/38849897html

最近在写的一个程序须要用到IIR滤波器，并且IIR滤波器的系数须要动态调整。所以就花了点时间研究IIR 滤波器的实现。python

之前用到的IIR滤波器的参数都是事先肯定好的，有个网站，只要把滤波器的参数特性输进去，直接就能生成须要的C代码。数组

http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html函数

一直都偷懒直接用这个网站的结果，因此手上也没积累什么有用的代码。此次就须要本身从头作起。测试

我面临的问题有两个：网站

1. 根据滤波器的参数（滤波器的类型、截止频率、滤波器的阶数等），计算出滤波器对应的差分方程的系数。spa

2. 利用获得的差分方程的系数构造一个能够工做的滤波器。.net

其中第一个问题，对于不一样类型的滤波器，好比Butterworth型、Bessel型等，滤波器系数的计算方法都不一样。这部分工做我还没作彻底，等我把常见的几种滤波器类型的系数计算方法都实现后再来写一篇文章。htm

这里就先写写第二个问题。IIR 滤波器对应的差分方程为：blog

 

相应的系统函数为：

 

这里默认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. /** \brief 将滤波器的内部状态清零，滤波器的系数保留 
  2.  * \return 
  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. /** \brief 
  25.  * 
  26.  * \param num 分子多项式的系数，升序排列,num[0] 为常数项 
  27.  * \param m 分子多项式的阶数 
  28.  * \param den 分母多项式的系数，升序排列,den[0] 为常数项 
  29.  * \param m 分母多项式的阶数 
  30.  * \return 
  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. /** \brief 滤波函数，采用直接I型结构 
  57.  * 
  58.  * \param data 传入输入数据 
  59.  * \return 滤波后的结果 
  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];  
  80.     }  
  81.     return m_py[0];  
  82. }  
  83.   
  84.   
  85. /** \brief 滤波函数，采用直接I型结构 
  86.  * 
  87.  * \param data[] 传入输入数据，返回时给出滤波后的结果 
  88.  * \param len data[] 数组的长度 
  89.  * \return 
  90.  */  
  91. void IIR_I::filter(double data[], int len)  
  92. {  
  93.     int i, k;  
  94.     for(k = 0; k < len; k++)  
  95.     {  
  96.         m_px[0] = data[k];  
  97.         data[k] = 0.0;  
  98.         for(i = 0; i <= m_num_order; i++)  
  99.         {  
  100.             data[k] = data[k] + m_pNum[i] * m_px[i];  
  101.         }  
  102.         for(i = 1; i <= m_den_order; i++)  
  103.         {  
  104.             data[k] = data[k] - m_pDen[i] * m_py[i];  
  105.         }  
  106.         // we get the y value now  
  107.         m_py[0] = data[k];  
  108.         for(i = m_num_order; i >= 1; i--)  
  109.         {  
  110.             m_px[i] = m_px[i-1];  
  111.         }  
  112.         for(i = m_den_order; i >= 1; i--)  
  113.         {  
  114.             m_py[i] = m_py[i-1];  
  115.         }  
  116.     }  
  117. }  
  118. /** \brief 滤波函数，采用直接I型结构 
  119.  * 
  120.  * \param data_in[] 输入数据 
  121.  * \param data_out[] 保存滤波后的数据 
  122.  * \param len 数组的长度 
  123.  * \return 
  124.  */  
  125. void IIR_I::filter(double data_in[], double data_out[], int len)  
  126. {  
  127.     int i, k;  
  128.     for(k = 0; k < len; k++)  
  129.     {  
  130.         m_px[0] = data_in[k];  
  131.         m_py[0] = 0.0;  
  132.         for(i = 0; i <= m_num_order; i++)  
  133.         {  
  134.             m_py[0] = m_py[0] + m_pNum[i] * m_px[i];  
  135.         }  
  136.         for(i = 1; i <= m_den_order; i++)  
  137.         {  
  138.             m_py[0] = m_py[0] - m_pDen[i] * m_py[i];  
  139.         }  
  140.         for(i = m_num_order; i >= 1; i--)  
  141.         {  
  142.             m_px[i] = m_px[i-1];  
  143.         }  
  144.         for(i = m_den_order; i >= 1; i--)  
  145.         {  
  146.             m_py[i] = m_py[i-1];  
  147.         }  
  148.         data_out[k] = m_py[0];  
  149.     }  
  150. }  

除此以外，还提供了个resp函数，这个函数能够计算滤波器的时域响应。而且不要求data_in与data_out 的数组长度相同。默认data_in[0] 以前的数据全为0，data_in[M-1]以后的数字所有为data_in[M-1]。所以，用这个函数计算阶跃响应输入数据只须要提供一个数据点就好了。而且这个函数不改变滤波器的内部状态。

 

[cpp]  view plain  copy
 
  1. /** \brief 计算 IIR 滤波器的时域响应，不影响滤波器的内部状态 
  2.  * \param data_in 为滤波器的输入，0 时刻以前的输入默认为 0，data_in[M] 及以后的输入默认为data_in[M-1] 
  3.  * \param data_out 滤波器的输出 
  4.  * \param M 输入数据的长度 
  5.  * \param N 输出数据的长度 
  6.  * \return 
  7.  */  
  8. void IIR_I::resp(double data_in[], int M, double data_out[], int N)  
  9. {  
  10.     int i, k, il;  
  11.     for(k = 0; k < N; k++)  
  12.     {  
  13.         data_out[k] = 0.0;  
  14.         for(i = 0; i <= m_num_order; i++)  
  15.         {  
  16.             if( k - i >= 0)  
  17.             {  
  18.                 il = ((k - i) < M) ? (k - i) : (M - 1);  
  19.                 data_out[k] = data_out[k] + m_pNum[i] * data_in[il];  
  20.             }  
  21.         }  
  22.         for(i = 1; i <= m_den_order; i++)  
  23.         {  
  24.             if( k - i >= 0)  
  25.             {  
  26.                 data_out[k] = data_out[k] - m_pDen[i] * data_out[k - i];  
  27.             }  
  28.         }  
  29.     }  
  30. }  

直接II型

[cpp]  view plain  copy
 
  1. /**< IIR 滤波器直接II型实现 */  
  2. class IIR_II  
  3. {  
  4. public:  
  5.     IIR_II();  
  6.     void reset();  
  7.     void setPara(double num[], int num_order, double den[], int den_order);  
  8.     void resp(double data_in[], int m, double data_out[], int n);  
  9.     double filter(double data);  
  10.     void filter(double data[], int len);  
  11.     void filter(double data_in[], double data_out[], int len);  
  12. protected:  
  13.     private:  
  14.     double *m_pNum;  
  15.     double *m_pDen;  
  16.     double *m_pW;  
  17.     int m_num_order;  
  18.     int m_den_order;  
  19.     int m_N;  
  20. };  
  21.   
  22. class IIR_BODE  
  23. {  
  24. private:  
  25.     double *m_pNum;  
  26.     double *m_pDen;  
  27.     int m_num_order;  
  28.     int m_den_order;  
  29.     complex<double> poly_val(double p[], int order, double omega);  
  30. public:  
  31.     IIR_BODE();  
  32.     void setPara(double num[], int num_order, double den[], int den_order);  
  33.     complex<double> bode(double omega);  
  34.     void bode(double omega[], int n, complex<double> resp[]);  
  35. };  

 

直接II型实现中所需的存储单元少了不少。另外，这两种实现的外部接口是彻底相同的。

 

[cpp]  view plain  copy
 
  1. IIR_II::IIR_II()  
  2. {  
  3.     //ctor  
  4.     m_pNum = NULL;  
  5.     m_pDen = NULL;  
  6.     m_pW = NULL;  
  7.     m_num_order = -1;  
  8.     m_den_order = -1;  
  9.     m_N = 0;  
  10. };  
  11.   
  12. /** \brief 将滤波器的内部状态清零，滤波器的系数保留 
  13.  * \return 
  14.  */  
  15. void IIR_II::reset()  
  16. {  
  17.     for(int i = 0; i < m_N; i++)  
  18.     {  
  19.         m_pW[i] = 0.0;  
  20.     }  
  21. }  
  22. /** \brief 
  23.  * 
  24.  * \param num 分子多项式的系数，升序排列,num[0] 为常数项 
  25.  * \param m 分子多项式的阶数 
  26.  * \param den 分母多项式的系数，升序排列,den[0] 为常数项 
  27.  * \param m 分母多项式的阶数 
  28.  * \return 
  29.  */  
  30. void IIR_II::setPara(double num[], int num_order, double den[], int den_order)  
  31. {  
  32.     delete[] m_pNum;  
  33.     delete[] m_pDen;  
  34.     delete[] m_pW;  
  35.     m_num_order = num_order;  
  36.     m_den_order = den_order;  
  37.     m_N = max(num_order, den_order) + 1;  
  38.     m_pNum = new double[m_N];  
  39.     m_pDen = new double[m_N];  
  40.     m_pW = new double[m_N];  
  41.     for(int i = 0; i < m_N; i++)  
  42.     {  
  43.         m_pNum[i] = 0.0;  
  44.         m_pDen[i] = 0.0;  
  45.         m_pW[i] = 0.0;  
  46.     }  
  47.     for(int i = 0; i <= num_order; i++)  
  48.     {  
  49.          m_pNum[i] = num[i];  
  50.     }  
  51.     for(int i = 0; i <= den_order; i++)  
  52.     {  
  53.         m_pDen[i] = den[i];  
  54.     }  
  55. }  
  56. /** \brief 计算 IIR 滤波器的时域响应，不影响滤波器的内部状态 
  57.  * \param data_in 为滤波器的输入，0 时刻以前的输入默认为 0，data_in[M] 及以后的输入默认为data_in[M-1] 
  58.  * \param data_out 滤波器的输出 
  59.  * \param M 输入数据的长度 
  60.  * \param N 输出数据的长度 
  61.  * \return 
  62.  */  
  63. void IIR_II::resp(double data_in[], int M, double data_out[], int N)  
  64. {  
  65.     int i, k, il;  
  66.     for(k = 0; k < N; k++)  
  67.     {  
  68.         data_out[k] = 0.0;  
  69.         for(i = 0; i <= m_num_order; i++)  
  70.         {  
  71.             if( k - i >= 0)  
  72.             {  
  73.                 il = ((k - i) < M) ? (k - i) : (M - 1);  
  74.                 data_out[k] = data_out[k] + m_pNum[i] * data_in[il];  
  75.             }  
  76.         }  
  77.         for(i = 1; i <= m_den_order; i++)  
  78.         {  
  79.             if( k - i >= 0)  
  80.             {  
  81.                 data_out[k] = data_out[k] - m_pDen[i] * data_out[k - i];  
  82.             }  
  83.         }  
  84.     }  
  85. }  
  86. /** \brief 滤波函数，采用直接II型结构 
  87.  * 
  88.  * \param data 输入数据 
  89.  * \return 滤波后的结果 
  90.  */  
  91. double IIR_II::filter(double data)  
  92. {  
  93.     m_pW[0] = data;  
  94.     for(int i = 1; i <= m_den_order; i++) // 先更新 w[n] 的状态  
  95.     {  
  96.         m_pW[0] = m_pW[0] - m_pDen[i] * m_pW[i];  
  97.     }  
  98.     data = 0.0;  
  99.     for(int i = 0; i <= m_num_order; i++)  
  100.     {  
  101.         data = data + m_pNum[i] * m_pW[i];  
  102.     }  
  103.     for(int i = m_N - 1; i >= 1; i--)  
  104.     {  
  105.         m_pW[i] = m_pW[i-1];  
  106.     }  
  107.     return data;  
  108. }  
  109. /** \brief 滤波函数，采用直接II型结构 
  110.  * 
  111.  * \param data[] 传入输入数据，返回时给出滤波后的结果 
  112.  * \param len data[] 数组的长度 
  113.  * \return 
  114.  */  
  115. void IIR_II::filter(double data[], int len)  
  116. {  
  117.     int i, k;  
  118.     for(k = 0; k < len; k++)  
  119.     {  
  120.         m_pW[0] = data[k];  
  121.         for(i = 1; i <= m_den_order; i++) // 先更新 w[n] 的状态  
  122.         {  
  123.             m_pW[0] = m_pW[0] - m_pDen[i] * m_pW[i];  
  124.         }  
  125.         data[k] = 0.0;  
  126.         for(i = 0; i <= m_num_order; i++)  
  127.         {  
  128.             data[k] = data[k] + m_pNum[i] * m_pW[i];  
  129.         }  
  130.   
  131.         for(i = m_N - 1; i >= 1; i--)  
  132.         {  
  133.             m_pW[i] = m_pW[i-1];  
  134.         }  
  135.     }  
  136. }  
  137. /** \brief 滤波函数，采用直接II型结构 
  138.  * 
  139.  * \param data_in[] 输入数据 
  140.  * \param data_out[] 保存滤波后的数据 
  141.  * \param len 数组的长度 
  142.  * \return 
  143.  */  
  144. void IIR_II::filter(double data_in[], double data_out[], int len)  
  145. {  
  146.     int i, k;  
  147.     for(k = 0; k < len; k++)  
  148.     {  
  149.         m_pW[0] = data_in[k];  
  150.         for(i = 1; i <= m_den_order; i++) // 先更新 w[n] 的状态  
  151.         {  
  152.             m_pW[0] = m_pW[0] - m_pDen[i] * m_pW[i];  
  153.         }  
  154.         data_out[k] = 0.0;  
  155.         for(i = 0; i <= m_num_order; i++)  
  156.         {  
  157.             data_out[k] = data_out[k] + m_pNum[i] * m_pW[i];  
  158.         }  
  159.   
  160.         for(i = m_N - 1; i >= 1; i--)  
  161.         {  
  162.             m_pW[i] = m_pW[i-1];  
  163.         }  
  164.     }  
  165. }  

 

测试

 

下面是几个测试例子，首先计算一个4阶切比雪夫低通滤波器的阶跃响应。

 

[python]  view plain  copy
 
  1. void resp_test(void)  
  2. {  
  3.     double b[5] = {0.001836, 0.007344, 0.011016, 0.007344, 0.001836};  
  4.     double a[5] = {1.0, -3.0544, 3.8291, -2.2925, 0.55075};  
  5.     double x[2] = {1.0, 1.0};  
  6.     double y[100];  
  7.   
  8.     IIR_II filter;  
  9.     filter.setPara(b, 4, a, 4);  
  10.     filter.resp(x, 2, y, 100);  
  11.   
  12.     for(int i= 0; i< 100; i++)  
  13.     {  
  14.         cout << y[i] << endl;  
  15.     }  
  16. }  


获得的结果以下：

 

 

一样是这个滤波器，计算输入信号为 delta 函数时的结果。

 

[cpp]  view plain  copy
 
  1. void filter_test(void)  
  2. {  
  3.     double b[5] = {0.001836, 0.007344, 0.011016, 0.007344, 0.001836};  
  4.     double a[5] = {1.0, -3.0544, 3.8291, -2.2925, 0.55075};  
  5.     double x[100], y[100];  
  6.     int len = 100;  
  7.     IIR_I filter;  
  8.     //IIR_II filter;  
  9.     filter.setPara(b, 4, a, 4);  
  10.   
  11.     for (int i = 0; i < len; i++)  
  12.     {  
  13.         x[i] = 0.0;  
  14.         y[i] = 0.0;  
  15.     }  
  16.     x[0] = 1.0;  
  17.     filter.filter(x, y, 100);  
  18.     filter.reset();  
  19.     filter.filter(x, 100);  
  20.   
  21.     for (int i = 0; i < len; i++)  
  22.     {  
  23.         cout << x[i] <<", " << y[i]<<  endl;  
  24.     }  
  25. }  



获得的结果以下：

 

相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程