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教程