IIR數字濾波器C語言

1.類比濾波器的設計       1.1巴特沃斯濾波器的次數         根據給定的參數設計類比濾波器，然後進行變數變換，求取數字濾波器的方法，稱為濾波器的間接設計。做為數字濾波器的設計基礎的類比濾波器，稱之為原型濾波器。這裡，我們首先介紹的是最簡單最基礎的原型濾波器，巴特沃斯低通濾波器。由於IIR濾波器不具有線性相位特性，因此不必考慮相位特性，直接考慮其振幅特性。
       在這裡，N是濾波器的次數，Ωc是截止頻率。從上式的振幅特性可以看出，這個是單調遞減的函數，其振幅特性是不存在紋波的。設計的時候，一般需要先計算跟所需要設計參數相符合的次數N。首先，就需要先由阻帶頻率，計算出阻帶衰減 將巴特沃斯低通濾波器的振幅特性，直接帶入上式，則有
最後，可以解得次數N為
當然，這裡的N只能為正數，因此，若結果為小數，則捨棄小數，向上取整。       1.2巴特沃斯濾波器的傳遞函數          巴特沃斯低通濾波器的傳遞函數，可由其振幅特性的分母多項式求得。其分母多項式

根據S解開，可以得到極點。這裡，為了方便處理，我們分為兩種情況去解這個方程。當N為偶數的時候，

這裡，使用了歐拉公式。同樣的，當N為奇數的時候，

同樣的，這裡也使用了歐拉公式。歸納以上，極點的解為

上式所求得的極點，是在s平面內，在半徑為Ωc的圓上等間距的點，其數量為2N個。為了使得其IIR濾波器穩定，那麼，只能選取極點在S平面左半平面的點。選定了穩定的極點之後，其類比濾波器的傳遞函數就可由下式求得。

      1.3巴特沃斯濾波器的實現（C語言）           首先，是次數的計算。次數的計算，我們可以由下式求得。           其對應的C語言程式為 [cpp]  view plain copy N = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /                log10 (Stopband/Cotoff) ));  

         然後是極點的選擇，這裡由於涉及到複數的操作，我們就聲明一個複數結構體就可以了。最重要的是，極點的計算含有自然指數函數，這點對於電腦來講，不是太方便，所以，我們將其替換為三角函數，

這樣的話，實部與虛部就還可以分開來計算。其代碼實現為 [cpp]  view plain copy typedef struct    {       double Real_part;       double Imag_Part;   } COMPLEX;         COMPLEX poles[N];      for(k = 0;k <= ((2*N)-1) ; k++)   {       if(Cotoff*cos((k+dk)*(pi/N)) < 0)       {           poles[count].Real_part = -Cotoff*cos((k+dk)*(pi/N));       　　poles[count].Imag_Part= -Cotoff*sin((k+dk)*(pi/N));                 count++;           if (count == N) break;       }   }   
       計算出穩定的極點之後，就可以進行傳遞函數的計算了。傳遞的函數的計算，就像下式一樣

這裡，為了得到類比濾波器的係數，需要將分母乘開。很顯然，這裡的極點不一定是整數，或者來說，這裡的乘開需要做複數運算。其複數的乘法代碼如下， [cpp]  view plain copy int Complex_Multiple(COMPLEX a,COMPLEX b,                    double *Res_Real,double *Res_Imag)          {          *(Res_Real) =  (a.Real_part)*(b.Real_part) - (a.Imag_Part)*(b.Imag_Part);          *(Res_Imag)=  (a.Imag_Part)*(b.Real_part) + (a.Real_part)*(b.Imag_Part);            return (int)1;    }   有了乘法代碼之後，我們現在簡單的情況下，看看其如何計算其濾波器係數。我們做如下假設 這個時候，其傳遞函數為

將其乘開，其大致的關係就像下圖所示一樣。

計算的關係一目瞭然，這樣的話，實現就簡單多了。高階的情況下也一樣，重複這種計算就可以了。其代碼為 [cpp]  view plain copy