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  Res[0].Real_part = poles[0].Real_part;     Res[0].Imag_Part= poles[0].Imag_Part;    Res[1].Real_part = 1;     Res[1].Imag_Part= 0;      for(count_1 = 0;count_1 < N-1;count_1++)   {     for(count = 0;count <= count_1 + 2;count++)     {         if(0 == count)    　{                 Complex_Multiple(Res[count], poles[count_1+1],                      &(Res_Save[count].Real_part),                      &(Res_Save[count].Imag_Part));         }         else if((count_1 + 2) == count)         {               Res_Save[count].Real_part  += Res[count - 1].Real_part;      　Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;         }            　else      　{                 Complex_Multiple(Res[count], poles[count_1+1],                      &(Res_Save[count].Real_part),                      &(Res_Save[count].Imag_Part));                  1   　　Res_Save[count].Real_part  += Res[count - 1].Real_part;     　　 Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;    　}     }    　　*(b+N) = *(a+N);   到此，我們就可以得到一個類比濾波器巴特沃斯低通濾波器了。

2.雙1次z變換       2.1雙1次z變換的原理         我們為了將類比濾波器轉換為數字濾波器的，可以用的方法很多。這裡著重說說雙1次z變換。我們希望通過雙1次z變換，建立一個s平面到z平面的映射關係，將類比濾波器轉換為數字濾波器。         和之前的例子一樣，我們假設有如下類比濾波器的傳遞函數。
將其做拉普拉斯逆變換，可得到其時間域內的連續微分方程式，
其中，x(t)表示輸入，y(t)表示輸出。然後我們需要將其離散化，假設其採樣周期是T，用差分方程去近似的替代微分方程，可以得到下面結果
然後使用z變換，再將其化簡。可得到如下結果
從而，我們可以得到了s平面到z平面的映射關係，即
由於所有的高階系統都可以視為一階系統的並聯，所以，這個映射關係在高階系統中，也是成立的。
然後，將關係式
帶入上式，可得
到這裡，我們可以就可以得到Ω與ω的對應關係了。

         這裡的Ω與ω的對應關係很重要。我們最終的目的設計的是數字濾波器，所以，設計時候給的參數必定是數字濾波器的指標。而我們通過間接設計設計IIR濾波器時候，首先是要設計類比濾波器，再通過變換，得到數字濾波器。那麼，我們首先需要做的，就是將數字濾波器的指標，轉換為類比濾波器的指標，基於這個指標去設計類比濾波器。另外，這裡的採樣時間T的取值很隨意，為了方便計算，一般取1s就可以。        2.2雙1次z變換的實現（C語言）          我們設計好的巴特沃斯低通濾波器的傳遞函數如下所示。       我們將其進行雙1次z變換，我們可以得到如下式子
可以看出，我們還是需要將式子乘開，進行合并同類項，這個跟之前說的演算法相差不大。其代碼為。 [cpp] view plain copy for(Count = 0;Count<=N;Count++)       {                    　for(Count_Z = 0;Count_Z <= N;Count_Z++)               {                    Res[Count_Z] = 0;                Res_Save[Count_Z] = 0;                 }                   Res_Save [0] = 1;            　 for(Count_1 = 0; Count_1 < N-Count;Count_1++)               {               　　for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)                   {                       if(Count_2 == 0)  Res[Count_2] += Res_Save[Count_2];       　　　　　　　　　        　else if((Count_2 == (Count_1+1))&&(Count_1 != 0))                            　    Res[Count_2] += -Res_Save[Count_2 - 1];    　　　　　　　　　　        else  Res[Count_2] += Res_Save[Count_2] - Res_Save[Count_2 - 1];   　　           　for(Count_Z = 0;Count_Z<= N;Count_Z++)                 　　　{                       　 Res_Save[Count_Z]  =  Res[Count_Z] ;                      　　Res[Count_Z]  = 0;                 　　　}                         }           for(Count_1 = (N-Count); Count_1 < N;Count_1++)               {                           for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)                   {                        if(Count_2 == 0) 　Res[Count_2] += Res_Save[Count_2];                       else if((Count_2 == (Count_1+1))&&(Count_1 != 0))                           　　　　　　　 Res[Count_2] += Res_Save[Count_2 - 1];                    else                            Res[Count_2] += Res_Save[Count_2] + Res_Save[Count_2 - 1];                  　}                     for(Count_Z = 0;Count_Z<= N;Count_Z++)                     {                          Res_Save[Count_Z]  =  Res[Count_Z] ;                      Res[Count_Z]  = 0;                     }               }               for(Count_Z = 0;Count_Z<= N;Count_Z++)           {                       *(az+Count_Z) += 　pow(2,N-Count) * (*(as+Count)) *   　　　　　　　　　　　      　　　　 　Res_Save[Count_Z];                   *(bz+Count_Z) +=  (*(bs+Count)) * Res_Save[Count_Z];                         }         }   到此，我們就已經實現了一個數字濾波器。

3.IIR濾波器的間接設計代碼（C語言） [cpp] view plain copy #include <stdio.h>   #include <math.h>   #include <malloc.h>   #include <string.h>         #define     pi     ((double)3.1415926)         struct DESIGN_SPECIFICATION   {       double Cotoff;          double Stopband;       double Stopband_attenuation;   };      typedef struct    {       double Real_part;       double Imag_Part;   } COMPLEX;            int Ceil(double input)   {        if(input != (int)input) return ((int)input) +1;        else return ((int)input);    }