巴特沃斯低通濾波器的實現(xiàn)方法研究

趙曉群，張 潔

(同濟大學電子與信息工程學院，上海201804)

巴特沃斯(Butterworth)濾波器是一種具有最大平坦幅度響應的低通濾波器，它在通信領域里已有廣泛應用。與貝塞爾(bessl)、契比雪夫(chebyshev)濾波器相比，巴特沃斯濾波器在線性相位、衰減斜率和加載特性3個方面具有特性均衡的優(yōu)點，因此在實際使用中，巴特沃斯濾波器已被列為首選［1］。

Matlab語言是一種面向科學與工程計算的語言，它的測試程序手段豐富，擴展能力強，內(nèi)涵豐富。信號處理工具箱(Signal Processing Toolbox)提供了設計巴特沃斯濾波器的函數(shù)，本文利用這些函數(shù)，進行了巴特沃斯濾波器的程序設計，并將其作為函數(shù)文件保存，可方便地進行調(diào)用。雖然Matlab擁有強大的數(shù)值運算功能，但編程運算效率并不高，對運行內(nèi)存要求巨大，通過本文可以證明，在相同條件下，Matlab對程序的處理時間比C增加了200倍。因此在工程應用中，常常將C與Matlab結(jié)合起來運用。

1 巴特沃斯低通濾波器的特性

巴特沃斯低通濾波器的幅頻特性模平方為

其中，N為濾波器的階數(shù)，Ωc為低通濾波器的截止頻率，當 Ω =Ωc時是濾波器的電壓-3dB點或稱半功率點。不同階次N的巴特沃斯濾波器特性如圖1，這一幅頻特性具有如下特點:

(1)最大平坦性:可以證明在Ω=0點，它的前2N-1階導數(shù)都等于零，這表示巴特沃斯濾波器在Ω=0附近一段范圍內(nèi)是非常平直的，它以原點的最大平坦性來逼近理想低通濾波器，“最平響應”即由此而得名。

(2)通帶、阻帶下降的單調(diào)性:這種濾波器具有良好的相頻特性。

(3)3dB的不變性:隨著N的增加，頻帶邊緣下降越陡峭，越接近理想特性。但不管N是多少，幅頻特性都通過-3dB點。當Ω=Ωc時，特性以20N dB/dec 速度下降［2］。

圖1 巴特沃斯濾波器幅頻特性

2 基于FFT的巴特沃斯低通濾波器設計方法

在連續(xù)時間系統(tǒng)中，可以利用卷積的方法求系統(tǒng)的零狀態(tài)響應，這時，首先把激勵信號分解為沖擊函數(shù)序列，然后令每一沖激函數(shù)單獨作用于系統(tǒng)求其沖激響應，最后把這些響應疊加即可得到系統(tǒng)對此激勵信號的零狀態(tài)響應。這個疊加的過程表現(xiàn)為求卷積積分。在離散時間系統(tǒng)中，可以采用大體相同的方法分析，由于離散信號本身就是一個不連續(xù)的序列，因此，激勵信號分解為脈沖序列的工作就很容易完成，對應每個樣值的激勵，系統(tǒng)得到對此樣值的響應，每一響應也是一個離散時間序列，把這些序列疊加即得到零狀態(tài)響應。

若長度為N1的序列x(n)與長度為N2的序列h(n)做線卷積得到y(tǒng)(n)為有線長序列，其長度為N1+N2-1。共需要進行N1N2次乘法運算，在N1=N2=N的情況下，需N2次乘法運算。

如果把求線性卷積改為求圓周卷積，并借助FFT技術(shù)，有可能減少卷積所需的運算工作量。如圖2給出直接卷積與快速卷積兩種方案原理圖。由圖2(b)可見，快速卷積的過程中，共需要兩次FFT，一次IFFT計算，相當于三次FFT的運算量。巴特沃斯低通濾波器的H(k)是知道的，數(shù)據(jù)已置于存儲器中，實際只需要兩個FFT運算量，如果假定N1=N2=N，經(jīng)補零后點數(shù)為N1+N2-1≈2N，因而需要2×(Nlog22N)次復數(shù)乘法運算。此外，為完成X(k)與H(k)兩序列相乘，還需要做2N次復乘。復乘運算全部次數(shù)為

顯然，隨著N值增大，式(2)的數(shù)字要比N2顯著減少。

圖2 直接卷積與快速卷積原理方框圖

當x(n)長度很長時，即N1?N2，通常不允許等x(n)全部采集齊后再進行卷積，否則使輸出相對于輸入有較長的延時，另外，若N1+N2-1太大，h(n)要補上太多的零點，很不經(jīng)濟，且FFT的計算時間也要很長。為此，采用分段卷積的方法，即把x(n)分成長度與h(n)相仿的一段段，分別求出每段卷積的結(jié)果，然后用相應的方式把它們結(jié)合起來，便是總的輸出。但是，N1可能很長，以至于趨于無窮大，以語音信號為例，如果不采用分段卷積的方法將遲遲不能給出結(jié)果，也無法找到那樣大的儲備來滿足N1的需要［2］。但是即使采用分段卷積，也要根據(jù)語音的分幀，一幀幀輸入語音，再對語音信號卷積，當前幀未輸入完，不能進行卷積，因此也是有著相當大的時延，對于實時性要求較高的語音通信系統(tǒng)，這種方法沒有本文給出的C語言實現(xiàn)的方法適用。

3 巴特沃斯低通濾波器的MATLAB實現(xiàn)

Matlab的信號處理工具箱提供了有關(guān)巴特沃斯濾波器的函數(shù) buttap，buttord，butter［3］。

設計巴特沃斯濾波器的程序?qū)崿F(xiàn)如下，Butter函數(shù)可在給定濾波器性能的情況下，選擇巴特沃斯濾波器的階數(shù)n和截止頻率ωc，從而可利用butter函數(shù)設計巴特沃斯濾波器的傳遞函數(shù)［n，ωc］=buttord(ωp，ωs，Rp，Rs，'s')可得到滿足性能的模擬巴特沃斯濾波器的最小階數(shù)n及截止頻率ωc，其中ωp為通帶的拐角頻率，ωs為阻帶的拐角頻率，ωp和ωs的單位均為rad/s;Rs為通帶區(qū)的最大波動系數(shù)，Rp為Rs阻帶區(qū)的最小衰減系數(shù)，Rp和 Rs的單位都為 dB。［b，a］=butter(n，ωc，'s')可設計截止頻率為ωc的n階低通模擬巴特沃斯濾波器.利用buttord函數(shù)、butter函數(shù)編制設計巴特沃斯低通濾波器的MATLAB函數(shù)文件butterdesign.m，其清單如下:

Function ［N，Wc，b，a］=butterdesign(Wp，Rp，Ws，As)

［N，Wc］=buttord(Wp，Ws，Rp，As，’s’);

［b，a］=butter(N，Wc，’s’);

［h，W］=freqs(b，a);

subplot(2，1，1);plot(W，abs(h));

subplot(2，1，2);plot(W，angle(h));

在低速率語音編碼中，基音周期的判別準確性將直接影響解碼端合成語音的質(zhì)量，MELP編碼器的基音提取中，首先要用截止頻率為1 kHz的6階巴特沃斯低通濾波器對語音信號進行濾s(n)波，消除語音幀中高頻成分對基音周期估算的影響［4］。在這里我們用Matlab實現(xiàn)用截止頻率為1 kHz的6階巴特沃斯低通濾波器對語音信號s(n)來濾波。采樣的一組語音信號:e=［5.405 1，64.149 0，80.667 6，115.356 8，48.432 1，-35.842 0，-90.292 7，-74.395 0，-97.388 3，-72.897 4 ］。調(diào)用函數(shù)［b，a］=butterπ(N，ωc)，N 為濾波器階數(shù) 6，ωc是歸一化截止頻率0.25。得到系數(shù)b=［0.001 1，0.006 3，0.015 8，0.021 0，0.015 8，0.0063，0.001 1］，a=［1.000 0，-2.978 5，4.136 1，-3.259 8，1.517 3，-0.391 1，0.043 4］。根據(jù)系數(shù)就可以調(diào)用濾波器函數(shù)對信號實行濾波，y=filter(b，a，e)，得到結(jié)果y=［0.005 7，0.118 5，0.904 3，3.977 9，11.971 4，26.804 9，46.517 8，63.225 6，65.475 1，45.087 3］。y為濾波后的信號。

4 巴特沃斯濾波器的C語言實現(xiàn)

利用重復計算由差分方程得出的遞推公式來實現(xiàn)線性時不變離散時間系統(tǒng)，根據(jù)線性時不變系統(tǒng)的直接I型和直接II型或規(guī)范型結(jié)構(gòu)，系統(tǒng)的輸入和輸出滿足如下差分方程［5］:

并具有如下系統(tǒng)函數(shù):

M=N=6，將之前得到的系數(shù)a和b對應帶入差分方程，在C語言下變成得到函數(shù)文件如下:

從結(jié)果可以看出，本文給出的基于C語言的設計方法與Matlab設計方法運行結(jié)果一致。

5 運行時空開銷的比較

可以用tic，toc函數(shù)計算下在Matlab里面運行的時間為0.000 303 s。用clock函數(shù)對C函數(shù)計算得到的運行時間為0.000 000 333 s。

將同樣數(shù)據(jù)的10個信號擴展到30個，Matlab運行時間0.000 394 s，C函數(shù)為0.000 001 33 s。將同樣的10個信號擴展到60個，Matlab運行時間為0.000 516 s，C 函數(shù)為0.000 002 33 s。

將同樣的信號擴展到150個，Matlab運行時間為0.001 090 s，C函數(shù)運行時間為0.000 005 33 s。

C語言運行時間效率要比用Matlab高得多，對于數(shù)字語音信號這樣的運算，一般信號數(shù)字運算量是相當巨大的，當數(shù)字信號運算量越大，C語言進行計算的效率優(yōu)勢就越明顯。

6 查看變量使用空間

在Matlab中命令窗口，用whos和whos global來查看所有局部和全局變量占用的內(nèi)存大?。?］，可以看到a和b是7為雙精度數(shù)值各56字節(jié)，e和y為輸入和輸出，是全部存在里面的，10個輸入的話就有80字節(jié)，輸出也是一樣的。

在C函數(shù)中，除了輸入輸出共設立了4個數(shù)組，p數(shù)組只是為了顯示用的，真正運行函數(shù)是不一定要用到的，a數(shù)組有6位占用48字節(jié)，b，x，y有7位雙精度值有56字節(jié)，但是C函數(shù)有最大的優(yōu)勢，就是它的輸入輸出占用的空間是可以變化的，這與Matlab必須先存在程序中不同，對于大量數(shù)值信號的計算來說，這樣是很費空間的，所以當計算數(shù)值數(shù)量越大時，C函數(shù)優(yōu)勢越明顯。

7 結(jié)論

本文給出了一種用C語言實現(xiàn)巴特沃斯低通濾波器的方法，并且通過與Matlab巴特沃斯濾波器一般方法以及另一種快速圓周卷積設計方法的比較中，得出本文用的C語言實現(xiàn)的方法在時間開銷以及運行時延等方面都比Matlab的方法要好，基于C語言的設計方法的處理時間是基于Matlab的設計方法的1/200。說明本文給出的算法在對于高實時性并且具有大輸入信號量的通信系統(tǒng)中是具有很大優(yōu)勢的。不足之處由于大輸入信號量的通信系統(tǒng)對于信號輸入量要求很大，也只能分段傳輸處理，所以還是有一定的延時。

本課題由水聲通信與海洋信息技術(shù)教育部重點實驗室(廈門大學)資助。

［1］李鐘慎.基于MATLAB設計巴特沃斯低通濾波器［J］.信息技術(shù)，2003，27(3):49-50.

［2］鄭君里，應啟行，楊為理.信號與系統(tǒng)［M］.北京:高等教育出版社，2000.

［3］張雪英.數(shù)字語音處理及MATLAB仿真［M］.北京:電子工業(yè)出版社，2010.

［4］韓笑蕾.(甚)低速率水聲語音編碼算法研究［D］.上海:同濟大學，2012:57-58.

［5］OPPENHEIM A V，SCHAFER R W，BUCK J R.Discrete-Time signal processing［D］.Cambridge:Massachusetts Institute of Technology，1999:284-285.