高階IIR濾波的設計與實現

剛離開學校時發現在學校學的很多東西在工作中都用不上。隨著時間的推移，當自己需要用這些東西時卻發現當初所學的知識都已經還給了老師。 對我來說，數字濾波器的設計就是例子。 

 

前段時間花了大概有一兩個月的時間來重溫當初在學校學習過的數字訊號處理相關內容，從歐拉公式到傳遞函數等，一點一滴慢慢的撿起來。以下就是我這一兩個月的工作成果^_^

 

由類比濾波器設計數字IIR濾波器： 脈衝響應不變法、雙線性變換法

 

我主要關注雙線性變換法。 

 

在選擇類比濾波器原型時，有如下幾種選擇： Butterworth/Chebyshev I/Chebyshev II/Elliptic/... 

 

N階IIR LPF： 

 

首先，我們要找到類比濾波器原型的幅度平方函數。這個在wiki上很容易找到。不要問我怎麼來的，我也不知道。 

 

 以Chebyshev I為例： 

     *                             1
     *      |H(jw)|^2 = --------------------------
     *                   sqrt(1 + e^2*Tn^2(w/wc))

由於有s = jw 

 

     *                             1
     *      H(s)*H(-s) = -----------------------
     *                    1 + e^2*Tn^2(-j*s/wc)

 

由上面的公式即可得到類比濾波器的極點，別問我怎麼解出來的，wiki上都有

 

     *      pk = (+-sinh(A)*sin(B) + jcosh(A)*cos(B))*wc
     *
     *         = Rp + j*Ip
     *
     *            1          1
     * where A = ---*arsinh(---), 
     *            n          e
     *
     *            PI   2*k-1
     *       B =  ---*-------, k = 1, 2, ... n/2
     *             2     n 

 

當然，上式中只顯示的n/2個極點均處於實軸之上，別一半極點位於實軸之下，且二者成共軛。為了在代碼中寫注釋方便，共軛我就寫成了pk~ 

 

 到這裡，所有與類比濾波器類型相關的部分都處理完了。接下來的處理對各種類比濾波器原型都適用： 

 

代碼裡面實現數字濾波器的關鍵是要把高階IIR濾波器轉換成1階或2階濾波器，這裡可以使用共軛複數的加法和乘法可以消去虛部的特性讓n階IIR濾波器轉換成n/2個2階結。 

 

繼續算啦

 

     *             n/2          1      
     *      H(s) = MUL -----------------------
     *             k=1 (1 - s/pk)*(1 - s/pk~)
     *
     *             n/2                     1
     *           = MUL ------------------------------------------
     *             k=1  1 - ((pk+pk~)/(pk*pk~))*s + s^2/(pk*pk~)
     *
     *             n/2                     1        
     *           = MUL ---------------------------------------------
     *             k=1  1 - ((2*Rp)/(Rp^2+Ip^2))*s + s^2/(Rp^2+Ip^2)
     *
     *             n/2          1
     *           = MUL -------------------
     *             k=1  1 + A1*s + A2*s^2

 

 到這裡就已經很簡單了，接下來就是雙線性變將類比濾波器轉換成數字濾波器：

     *             n/2             1 + 2*z^-1 + z^-2
     *      H(z) = MUL --------------------------------------------
     *             k=1  (1+A1+A2) + (2-2*A2)*z^-1 + (1-A1+A2)*z^-2

 

 接下來用代碼實現就easy啦。

 

LPF OK了，HPF用頻率轉換方法也很easy。BPF有些麻煩，需要用4階結。 

 

 ======================

 

接下來是一些概念：

 

 自然頻率： Hz

規一化的數字頻率：   自然頻率/採樣率

數字頻率： 2*PI*規一化的數字頻率 

類比頻率： tan(PI*規一化的數字頻率)

 

 給大家猜猜： 上面的wc是什麼頻率 

 

  =============================

最後是一點點小體會： 

 

在做公式推導時，多用變數代替某一階段的結果，這樣寫代碼時清晰可見。而且當發現公式推匯出錯時，不會造成整個代碼重寫。當然，在做效能最佳化時又另當別論啦。