IIR數字濾波器的實現（C語言）

• 經典濾波器和數字濾波器

　　一般濾波器可以分為經典濾波器和數字濾波器。

1. 經典濾波器：假定輸入信號中的有用成分和希望去除的成分各自占有不同的頻帶。如果信號和噪聲的頻譜相互重迭，經典濾波器無能為力。比如 FIR 和 IIR 濾波器等。　　
2. 現代濾波器：從含有噪聲的時間序列中估計出信號的某些特征或信號本身。現代濾波器將信號和噪聲都視為隨機信號。包括 Wiener Filter、Kalman Filter、線性預測器、自適應濾波器等

• Z變換和差分方程

　　在連續系統中采用拉普拉斯變換求解微分方程，並直接定義了傳遞函數，成為研究系統的基本工具。在采樣系統中，連續變量變成了離散量，將Laplace變換用於離散量中，就得到了Z變換。和拉氏變換一樣，Z變換可用來求解差分方程，定義Z傳遞函數成為分析研究采樣系統的基本工具。

　　對於一般常用的信號序列，可以直接查表找出其Z變換。相應地，也可由信號序列的Z變換查出原信號序列，從而使求取信號序列的Z變換較為簡便易行。

　　Z變換有許多重要的性質和定理：

• 線性定理

　　設a，a₁，a₂為任意常數，連續時間函數f(t)，f₁(t)，f₂(t)的Z變換分別為F(z)，F₁(z)，F₂(z)，則有$$\mathbf{Z}[af(t)]=aF(z)$$ $$ \mathbf{Z}[a_1f_1(t)+a_2f_2(t)]=a_1F_1(z)+a_2F_2(z)$$

• 滯后定理

　　設連續時間函數在t<0時，f(t)=0，且f(t)的Z變換為F(z)，則有$$\mathbf{Z}[f(t-kT)]=z^{-k}F(z)$$

 

　　應用Z變換求解差分方程的一個例子：已知系統的差分方程表達式為$y(n)-0.9y(n-1)=0.05u(n)$，若邊界條件$y(-1)=1$，求系統的完全響應。

　　解：方程兩端取z變換$$Y(z)-0.9[z^{-1}Y(z)+y(-1)]=0.05\frac{z}{z-1}\\Y(z)=\frac{0.05z^2}{(z-1)(z-0.9)}+\frac{0.9y(-1)z}{z-0.9}$$

　　可得$$\frac{Y(z)}{z}=\frac{A_1z}{z-1}+\frac{A_2z}{z-0.9}$$

　　其中A₁=0.5，A₂=0.45，於是$y(n)=0.5+0.45 \times(0.9)^n \quad(n\geq0)$

•  IIR數字濾波器的差分方程和系統函數

　　 IIR數字濾波器是一類遞歸型的線性時不變因果系統，其差分方程可以寫為：$$y(n)=\sum_{i=0}^{M}a_ix(n-i)+\sum_{i=1}^{N}b_iy(n-i)$$

　　進行Z變換，可得：$$Y(z)=\sum_{i=0}^{M}a_iz^{-i}X(z)+\sum_{i=1}^{N}b_iz^{-i}Y(z)$$

　　於是得到IIR數字濾波器的系統函數：$$H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{i=0}^{M}a_iz^{-i}}{1-\sum_{i=1}^{N}b_iz^{-i}}=a_0\frac{\prod_{i=1}^{M}(1-c_iz^{-1})}{\prod_{i=1}^{N}(1-d_iz^{-1})}$$

　　其中c_i為零點而d_i為極點。H(z)的設計就是要確定系數、或者零極點，以使濾波器滿足給定的性能指標。

•  IIR數字濾波器結構

 　　數字濾波器的功能本質上是將一組輸入數字序列通過一定的運算后轉變為另一組輸出數字序列。濾波器系統函數可以表達為多種不同的形式，每一種對應着不同的算法，也就對應着不同的實現結構。例如：$$H(z)=\frac{1}{1-0.3z^{-1}-0.4z^{-2}}$$

　　可以分解為：$$H(z)=\frac{1}{1-0.8z^{-1}}\cdot \frac{1}{1+0.5z^{-1}}$$

　　或$$H(z)=\frac{0.6154}{1-0.8z^{-1}}+\frac{0.3846}{1+0.5z^{-1}}$$

　　上述同一系統的三種不同描述形式就對應了不同的實現結構，或者說不同的濾波器結構可以實現相同的傳遞函數。IIR濾波器常見的結構形式有直接Ⅰ型、直接Ⅱ型（典范型）、級聯型、並聯型。通過差分方程能夠畫出包含反饋結構的數字網絡稱為直接型。

　　直接Ⅰ型濾波器的網絡結構可以根據差分方程很直觀地畫出（The Direct-Form I structure implements the feed-forward terms first followed by the feedback terms）：

 　　可以看出直接Ⅰ型濾波器需要N+M個延時單元（N≥M）。直接Ⅱ型結構是對直接Ⅰ型的變型，將上面網絡的零點與極點的級聯次序互換，並將延時單元合並。實現N階濾波器只需要N個延時單元（The Direct-Form II structure uses fewer delay blocks than Direct-Form I），故稱為典范型。

　　直接Ⅱ型看上去不那么直觀，可以通過下圖進行理解。我們可以將整個濾波器系統看成A、B兩個子系統串聯而成，由於為線性系統因此交換順序不影響最終輸出結果，傳遞函數可寫為：$$H(z)=B(z)\cdot \frac{1}{A(z)}=\frac{1}{A(z)}\cdot B(z)$$

　　直接型優點：直接型結構簡單，用的延遲器較少（N和M中較大者的個數）；缺點：系數a_k，b_k對濾波器性能的控制關系不直接，因此調整不方便；具體實現濾波器時a_k，b_k的量化誤差將使濾波器的頻率響應產生很大的改變，甚至影響系統的穩定性。直接型結構一般用於實現低階系統。

　　級聯型：將系統傳遞函數H(z)因式分解為多個二階子系統，系統函數就可以表示為這些二階子系統傳遞函數的乘積。實現時將每個二階子系統用直接型實現，整個系統函數用二階環節的級聯實現。

　　 並聯型：與級聯型類似，用部分分式展開法將系統函數表示為二階子系統傳遞函數的和。每個二階子系統仍然用直接型實現，整個系統函數用二階環節的並聯實現。

 

---

　　在IIR濾波器設計過程中，通常利用模擬濾波器來設計數字濾波器，要先根據濾波器的性能指標設計出相應的模擬濾波器的系統函數H(s)，然后由H(s)經變換得到所需要的數字濾波器的系統函數H(z)。常用的變換方法有沖激響應不變法和雙線性變換法。下面使用MATLAB等工具直接生成數字濾波器系數：

　　在MATLAB命令行中輸入fdatool打開濾波器設計工具箱，為了便於分析，我們先從設計一個簡單的2階低通濾波器。Design Method用於選擇IIR濾波器還是FIR濾波器，這里我們選擇IIR濾波器，類型選擇Butterworth，當然也可以選擇其他類型，不同類型的頻率響應不同，選擇后默認的濾波器結構是直接II型。ResponseType用於選擇低通、高通、帶通、帶阻等類型，選擇低通濾波“Lowpass”。Frequency Specifications用於設置采樣頻率以及截止頻率，這里填入200以及20，即采樣率為200Hz，20Hz以上的頻率將被濾除掉。Fiter Order 選擇濾波器階數，為了簡單起見，先選擇一個2階濾波器做實驗。

 　　參數設置好后點擊Design filter按鈕，將按要求設計濾波器。默認生成的IIR濾波器類型是Direct-Form II,Second-Order Sections（直接Ⅱ型，每個Section是一個二階濾波器），在工具欄上點擊Filter Coefficients圖標或菜單欄上選擇Analysis→Filter Coefficients可以查看生成的濾波器系數。

　　MATLAB中二階濾波器差分方程公式如下（注意反饋項符號為負號）：$$y[n]=b_0 \cdot x[n]+b_1 \cdot x[n-1] + b_2 \cdot x[n-2] -a_1 \cdot y[n-1] -a_2 \cdot y[n-2]$$

　　高階IIR濾波器的實現是采用二階濾波器級聯的方式來實現的。默認情況下，Filter Coefficients把結果分成多個2階Section顯示，其中還有增益。增益的目的是為了保證計算的精度和系統的穩定性。選擇[edit]→[Convert to Single Section]，這時候系數變成我們熟悉的形式：

　　按照上面的公式，濾波器差分方程為：y[n] = 0.06745527*x[n] + 0.134910547*x[n-1] + 0.06745527*x[n-2] - (-1.1429805025)*y[n-1] - (0.412801596)*y[n-2]

　　濾波器設計完成后還可以生成Simulink模型進行仿真：按照下圖中數字標號進行，第一步點擊左邊Realize Model圖標，第二步勾選“Build model using basic elements”這一項，右邊四個灰色的項將自動打鈎，最后點擊“Realize Model”，matlab將自動生成濾波器模型，在彈出的窗口中雙擊模型可以觀察該模型的內部結構。

　　下面是按照設計要求生成的2階濾波器直接Ⅰ型的結構：

　　Direct-Form I

　　下面是直接Ⅱ型的內部結構：

Direct-Form II

　　使用生成的濾波器搭建一個簡答的測試模型：將兩個幅值為1，頻率分別為10Hz、50Hz的正弦波疊加，輸入濾波器后觀察濾波前后的波形。仿真時間設為1s，仿真參數中求解器類型設為固定步長，求解器選擇discrete（它適用於離散無連續狀態的系統），步長設為0.005s（200Hz）

　　點擊Run按鈕開始進行仿真：

 　　打開示波器結果如下圖所示：上面一欄是不同頻率疊加的波形，下面是10Hz正弦波和濾波后得到的波形的對比。由於50Hz正弦波頻率高於濾波器截止頻率20Hz，因此被濾除，同時濾波也產生了一定的滯后和失真。

　　知道了差分方程的形式並通過MATLAB得到濾波器系數后很容易寫出相應的代碼來實現數字濾波，另外還有一個網站能根據設計指標直接生成C代碼：http://www-users.cs.york.ac.uk/~fisher/mkfilter/

　　根據前面的設計指標，在網頁上填入相應參數后提交，會得到下面的C語言代碼。簡單修改后就可以使用：

#define NZEROS 2
#define NPOLES 2
#define GAIN   1.482463775e+01

static float xv[NZEROS+1], yv[NPOLES+1];

static void filterloop()
{ 
　　for (;;)
   {
　　　　 xv[0] = xv[1]; xv[1] = xv[2]; 
        xv[2] = next input value / GAIN;
        yv[0] = yv[1]; yv[1] = yv[2]; 
        yv[2] =   (xv[0] + xv[2]) + 2 * xv[1] + ( -0.4128015981 * yv[0]) + (  1.1429805025 * yv[1]);
        next output value = yv[2];
    }
}

　　

　　在LabVIEW中為了自己實現IIR濾波器可以使用反饋結點來存儲數據，下面的程序框圖實現了與MATLAB模型相同的功能：

　　前面板波形圖如下圖所示：

 

 

參考：

Butterworth / Bessel / Chebyshev Filters

ARM官方DSP庫IIR濾波器的實現(STM32)

ARM官方DSP庫濾波器基礎知識

手把手教你用matlab生成STM32官方IIR濾波器的系數(二)

手把手教你用matlab生成STM32官方IIR濾波器的系數(三)

A Collection of Useful C++ Classes for Digital Signal Processing

數字濾波器

IIR Filters

簡單常用濾波算法C語言實現