網(wǎng)站內(nèi)容注意事項(xiàng)福州關(guān)鍵詞排名優(yōu)化
? ? ? ? 我們常說的經(jīng)典濾波器是根據(jù)傅里葉分析和變換設(shè)計(jì)出來的，只允許一定頻率范圍內(nèi)的信號成分正常通過，而阻止另一部分頻率成分通過。按照最佳逼近特性或者濾波通帶特性分類，主要為巴特沃斯濾波器（Butterworth）、切比雪夫?yàn)V波器（Chebyshev）、貝塞爾濾波器（Bessel）和橢圓濾波器（Elliptic）四種。每種MATLAB都有相應(yīng)的函數(shù)，用起來也比較方便，但是卻缺少C/C++的程序，于是自己仔細(xì)研究了每種濾波器的特性和原理，并且部分濾波器實(shí)現(xiàn)了C語言的代碼化，接下來的時(shí)間會對這些濾波器的原理和C語言的實(shí)現(xiàn)進(jìn)行介紹。
該系列均以低通濾波器為原型來介紹，其他類型的濾波器可以由低通濾波器通過頻率變換轉(zhuǎn)換得到，這里不過多介紹。低通濾波器的主要性能指標(biāo)有以下幾個：通帶截止頻率fp、阻帶截止頻率fs、通帶衰減 ( Ap)、阻帶衰減 ( As) 以及歸一化頻率時(shí)需要用到的-3dB的轉(zhuǎn)折頻率fc。
1.?Butterworth濾波器原理
? ? ?Butterworth濾波器因其在通帶內(nèi)的幅值特性具有最大平坦的特性而聞名，是四種經(jīng)典濾波器中最簡單的，巴特沃斯濾波器只需要兩個參數(shù)表征，濾波器的階數(shù)N和-3dB處的截止頻率
。其幅度平方函數(shù)為：
??????????????????????????????????????????????????????????????????????????? 
?????????????????????????????? 
N是濾波器的階數(shù)，從幅度平方函數(shù)可以看出，N階濾波器有2N個極點(diǎn)，而且這2N個極點(diǎn)均布在一個圓上，圓的半徑為，稱之為Butterworth圓，Butterworth濾波器系統(tǒng)是一個線性系統(tǒng)，要使其穩(wěn)定，其極點(diǎn)必須位于S平面的左半平面，所以取左半平面內(nèi)的N個極點(diǎn)作為濾波器的極點(diǎn)，濾波器就是穩(wěn)定的了，求出極點(diǎn)之后，計(jì)算模擬濾波器的系數(shù)as、bs，然后通過雙線性變換（不懂得自行查書）由模擬域到數(shù)字域，求出系數(shù)az和bz?。最后通過差分方程就可以計(jì)算濾波結(jié)果了。
????????????????????????????? ?? ???????????? 
2.?C語言實(shí)現(xiàn)
A．求階數(shù)
公式為：
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
代碼：
N = ceil(0.5*( log10 (( pow (10, Stopband_attenuation/10) - 1)/( pow (10, Passband_attenuation/10) - 1)) / log10 (Stopband/Passband) ));
B．求極點(diǎn)
公式：
代碼：
for(k = 0;k <= ((2*N)-1) ; k++){if(Cutoff*cos((2*k+N-1)*(pi/(2*N))) < 0.0){ poles[count].x = -Cutoff*cos((2*k+N-1)*(pi/(2*N)));poles[count].y = -Cutoff*sin((2*k+N-1)*(pi/(2*N))); count++;if (count == N) break;}} C．計(jì)算模擬濾波器系數(shù)
公式：
代碼（這部分需要自己推導(dǎo)，多算幾步就找到規(guī)律了）：
Res[0].x = poles[0].x; Res[0].y = poles[0].y;Res[1].x = 1; Res[1].y= 0;for(count_1 = 0;count_1 < N-1;count_1++)//N個極點(diǎn)相乘次數(shù){for(count = 0;count <= count_1 + 2;count++){if(0 == count){Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] ); }else if((count_1 + 2) == count){Res_Save[count].x += Res[count - 1].x;Res_Save[count].y += Res[count - 1].y; } else {Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] ); Res_Save[count].x += Res[count - 1].x;Res_Save[count].y += Res[count - 1].y; }}for(count = 0;count <= N;count++)//Res[i]=a(i),i越大次數(shù)越高{Res[count].x = Res_Save[count].x; Res[count].y = Res_Save[count].y;*(a + N - count) = Res[count].x;} }*(b+N) = *(a+N);D．雙線性變換
? ? ? ? 用下式進(jìn)行替換
?中的
變量，得到
，然后類似上面的計(jì)算方法計(jì)算Z域的系數(shù)az和bz，其中T為采樣周期，但是因?yàn)樵谟?jì)算中會被約去，所以簡化計(jì)算，這里取1。
公式：
代碼：
int Count = 0,Count_1 = 0,Count_2 = 0,Count_Z = 0;double *Res, *Res_Save;Res = new double[N+1]();Res_Save = new double[N+1](); memset(Res, 0, sizeof(double)*(N+1));memset(Res_Save, 0, sizeof(double)*(N+1));for(Count_Z = 0;Count_Z <= N;Count_Z++){*(az+Count_Z) = 0;*(bz+Count_Z) = 0;}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++)//計(jì)算（1-Z^-1）^N-Count的系數(shù),{ //Res_Save[]=Z^-1多項(xiàng)式的系數(shù)，從常數(shù)項(xiàng)開始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++)//計(jì)算(1-Z^-1)^N-Count*（1+Z^-1）^Count的系數(shù),{ //Res_Save[]=Z^-1多項(xiàng)式的系數(shù)，從常數(shù)項(xiàng)開始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]; } }//最外層for循環(huán)for(Count_Z = N;Count_Z >= 0;Count_Z--){*(bz+Count_Z) = (*(bz+Count_Z))/(*(az+0));*(az+Count_Z) = (*(az+Count_Z))/(*(az+0));}E．由由差分方程計(jì)算濾波結(jié)果
采用濾波器直接II型結(jié)構(gòu)，可以減少一半的中間緩存內(nèi)存，具體代碼如下：
double FiltButter(double *pdAz, //濾波器參數(shù)表1double *pdBz, //濾波器參數(shù)表2int nABLen, //參數(shù)序列的長度double dDataIn,//輸入數(shù)據(jù)double *pdBuf) //數(shù)據(jù)緩沖區(qū)
{int i;int nALen;int nBLen;int nBufLen;double dOut;if( nABLen<1 )return 0.0;//根據(jù)參數(shù),自動求取序列有效長度nALen = nABLen;for( i=nABLen-1; i; --i ){if( *(pdAz+i) != 0.0 )//從最后一個系數(shù)判斷是否為0{nALen = i+1; break;}}//printf("%lf ", nALen);if( i==0 ) nALen = 0;nBLen = nABLen;for( i=nABLen-1; i; --i ){if( *(pdBz+i) != 0.0 ){nBLen = i+1;break;}}//printf("%lf ", nBLen);if( i==0 ) nBLen = 0;//計(jì)算緩沖區(qū)有效長度nBufLen = nALen;if( nALen < nBLen)nBufLen = nBLen;//濾波: 與系數(shù)a卷乘dOut = ( *pdAz ) * dDataIn; // a(0) * x(i) for( i=1; i<nALen; i++) // a(i) * w(n-i),i=1toN{dOut -= *(pdAz+i) * *(pdBuf + (nBufLen - 1) - i);} //卷乘結(jié)果保存為緩沖序列的最后一個*(pdBuf + nBufLen - 1) = dOut;//濾波: 與系數(shù)b卷乘dOut = 0.0;for( i=0; i<nBLen; i++) // b(i) * w(n-i){ dOut += *(pdBz+i) * *(pdBuf + (nBufLen - 1) - i);}//丟棄緩沖序列中最早的一個數(shù), 最后一個數(shù)清零for( i=0; i<nBufLen-1; i++){*(pdBuf + i) = *(pdBuf + i + 1);}*(pdBuf + nBufLen - 1) = 0;//返回輸出值return dOut;
}VC6.0開發(fā)環(huán)境濾波結(jié)果如下：
至此就完成了Butterworth濾波器的設(shè)計(jì)過程，完整代碼請自行下載。
?
參考資料：http://blog.csdn.net/zhoufan900428/article/details/9069475
源碼地址：http://download.csdn.net/detail/zhwzhaowei/9830114
?
?