趙曉群，張 潔

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

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

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

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

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

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

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

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

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

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

2 基于FFT的巴特沃斯低通濾波器設(shè)計(jì)方法

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

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

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

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

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

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

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

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

設(shè)計(jì)巴特沃斯濾波器的程序?qū)崿F(xiàn)如下，Butter函數(shù)可在給定濾波器性能的情況下，選擇巴特沃斯濾波器的階數(shù)n和截止頻率ωc，從而可利用butter函數(shù)設(shè)計(jì)巴特沃斯濾波器的傳遞函數(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')可設(shè)計(jì)截止頻率為ωc的n階低通模擬巴特沃斯濾波器.利用buttord函數(shù)、butter函數(shù)編制設(shè)計(jì)巴特沃斯低通濾波器的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ǔn)確性將直接影響解碼端合成語音的質(zhì)量，MELP編碼器的基音提取中，首先要用截止頻率為1 kHz的6階巴特沃斯低通濾波器對語音信號進(jìn)行濾s(n)波，消除語音幀中高頻成分對基音周期估算的影響［4］。在這里我們用Matlab實(shí)現(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ù)對信號實(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語言實(shí)現(xiàn)

利用重復(fù)計(jì)算由差分方程得出的遞推公式來實(shí)現(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對應(yīng)帶入差分方程，在C語言下變成得到函數(shù)文件如下:

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

5 運(yùn)行時空開銷的比較

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

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

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

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

6 查看變量使用空間

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

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

7 結(jié)論

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

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

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

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

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

［4］韓笑蕾.(甚)低速率水聲語音編碼算法研究［D］.上海:同濟(jì)大學(xué)，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.