基于Matlab的IIR Butterworth低通數(shù)字濾波器設(shè)計(jì).doc

問題：信號(hào)濾波是一個(gè)非常重要的信號(hào)處理手段，大量的信號(hào)處理系統(tǒng)中，信號(hào)濾波的質(zhì)量將直接影響到系統(tǒng)整體性能。現(xiàn)在有家醫(yī)院向你提出一個(gè)問題，她們的心電圖信號(hào)記錄儀器由于受到了比較嚴(yán)重的電源干擾(50Hz附近)，記錄到的心電圖ECG波形有較大的畸變，請問如何設(shè)計(jì)一個(gè)干擾濾除系統(tǒng)，讓記錄下的ECG數(shù)據(jù)繼續(xù)可以進(jìn)行疾病診斷。(ECG的波形信息非常重要，設(shè)計(jì)濾波器需要注意這一點(diǎn))基于Matlab的IIR Butterworth低通數(shù)字濾波器設(shè)計(jì)0 引言心電信號(hào)是心臟電活動(dòng)在人體體表的表現(xiàn)，它一般比較微弱，其頻率為0.05Hz 100Hz（能量集中在0. 05Hz 44Hz），幅度為幾百微伏至幾毫伏。因此，在心電信號(hào)的采集過程中，極易受到內(nèi)、外界環(huán)境的干擾，其中50Hz的工頻干擾尤為突出。心電圖信號(hào)的干擾對心電圖數(shù)據(jù)分析和壓縮有一定影響。如何消除50Hz工頻干擾，成為處理心電信號(hào)的首要任務(wù)。為了消除ECG數(shù)據(jù)中的主要干擾，保留有用信號(hào)或者從干擾中“抽取”有用信號(hào)，人們提出了許多方法，諸如，平滑濾波，中心頻率固定的帶阻濾波器，自適應(yīng)濾波，低通數(shù)字濾波等等。這些信號(hào)處理方法各有優(yōu)缺點(diǎn)，平滑濾波算法簡單，處理速度快，濾波效果較好，但存在一定的削峰作用；中心頻率固定的帶阻濾波器原理簡單，能較大程度地抑制50Hz工頻干擾，但存在“群延時(shí)”現(xiàn)象；自適應(yīng)濾波器的中心頻率能夠跟隨工頻信號(hào)的頻率幅度變化而自動(dòng)調(diào)節(jié)并抵消工頻干擾，但需要附加參考信號(hào)通道，算法相對復(fù)雜，難以用于實(shí)時(shí)處理?？紤]到工頻干擾是50Hz的低頻信號(hào)，假設(shè)要處理的ECG信號(hào)為0.05HZ44HZ段信號(hào)，且Matlab軟件中有一個(gè)數(shù)字濾波器軟件包，可以直接進(jìn)行調(diào)用處理信號(hào)，所以本文設(shè)計(jì)IIR巴特沃斯(Butterworth)低通數(shù)字濾波器來實(shí)現(xiàn)人體ECG信號(hào)的濾波。1 數(shù)字濾波器介紹數(shù)字濾波器是一種用來過濾時(shí)間離散信號(hào)的數(shù)字系統(tǒng)，它可以用軟件(計(jì)算機(jī)程序)或用硬件來實(shí)現(xiàn)，而且在兩種情況下都可以用來過濾實(shí)時(shí)信號(hào)或非實(shí)時(shí)信號(hào)(記錄信號(hào))。盡管數(shù)字濾波器這個(gè)名稱一直到六十年代中期才出現(xiàn)。但是隨著科學(xué)技術(shù)的發(fā)展及計(jì)算機(jī)的更新普及，數(shù)字濾波器有著很好的發(fā)展前景。有些時(shí)候，它可以完全取代模擬濾波器。 數(shù)字濾波器具有如下的優(yōu)點(diǎn)：它沒有漂移，它能處理低頻信號(hào)，頻率響應(yīng)特征可做成極接近于理想的特性，它可以做成沒有插人損耗和有線性相位特性，可相當(dāng)簡單地獲得自適應(yīng)濾波，濾波器的設(shè)計(jì)者可以控制數(shù)字字長，因而可以精確地控制濾波的精度。其中，它最重要的優(yōu)點(diǎn)是：隨著濾波器參數(shù)的改變，很容易改變?yōu)V波器的性能。這一特點(diǎn)允許我們只用一種程序?yàn)V波器完成多重濾波任務(wù)。 本文采用雙線性變換法設(shè)計(jì)IIR巴特沃斯低通數(shù)字濾波器實(shí)現(xiàn)了對工頻干擾信號(hào)的有效濾波。2 巴特沃斯（Butterworth）低通數(shù)字濾波器的特性所謂低通濾波器，就是只能讓低頻分量通過。在低通濾波器的設(shè)計(jì)中，有許多不同的逼近方法：巴特沃斯逼近，切比雪夫逼近，橢圓逼近方法等。其中，最簡單的逼近為巴特沃斯逼近。巴特沃斯低通濾波器的平方幅度響應(yīng)為其中，n為濾波器的階數(shù)，為低通濾波器的截止頻率。 該濾波器具有 一些特殊的性質(zhì)： 對所有的n,都有當(dāng)時(shí), ； 對所有的n,都有當(dāng)時(shí), ；是的單調(diào)遞減函數(shù),即不會(huì)出現(xiàn)幅度響應(yīng)的起伏； 當(dāng)時(shí),巴特沃斯濾波器趨向于理想的低通濾波器； 在處平方幅度響應(yīng)的各級(jí)導(dǎo)數(shù)均存在且等于0，因此在該點(diǎn)上取得最大值，且具有最大平坦特性。 圖1展示了2階、4階、8階巴特沃斯低通濾波器的幅頻特性?？梢婋A數(shù)n越高，其幅頻特性越好，低頻檢測信號(hào)保真度越高，過渡帶變窄，即衰減加劇，但半功率點(diǎn)不變。圖1 巴特沃斯低通濾波器的幅頻特性3 IIR數(shù)字濾波器的性質(zhì)無限長沖激響應(yīng)( IIR) 數(shù)字濾波器是數(shù)字濾器的一種，數(shù)字濾波器根據(jù)其沖激響應(yīng)函數(shù)的時(shí)域特性,還包括有限長沖激響應(yīng)( FIR)數(shù)字濾波器。IIR數(shù)字濾波器的特征是：具有無限持續(xù)時(shí)間沖激響應(yīng)；需要用遞歸模型來實(shí)現(xiàn)，這可以從其差分方程得出，也可以從其系數(shù)函數(shù)為：得出。 數(shù)字巴特沃思濾波器屬于IIR濾波器，該類濾波器具有特定的性質(zhì)和設(shè)計(jì)方法。目前比較成熟的IIR數(shù)字濾波器設(shè)計(jì)方法有兩種：1）直接法目前所用的方法主要是：零極點(diǎn)累試法、頻域幅度平方誤差最小法和時(shí)域單位脈沖響應(yīng)逼近法。直接法的最大優(yōu)點(diǎn)在于可以設(shè)計(jì)任意幅頻特征的濾波器。2)間接法，目前所用的方法主要是：沖激響應(yīng)不變法、階躍響應(yīng)不變法和雙線性法。它們都是借助于 已經(jīng)成熟的現(xiàn)有低通濾波器原型進(jìn)行設(shè)計(jì)，即對數(shù)字低通數(shù)字濾波器，先將數(shù)字低通濾波器的技術(shù)指標(biāo)按希望的設(shè)計(jì)方法轉(zhuǎn)換為模擬低通濾波器的技術(shù)指標(biāo)，再按指定的模擬低通濾波器的類型設(shè)計(jì)模擬濾波器H(s)，然后，將模擬濾波器的系統(tǒng)函數(shù)H(s)從s平面轉(zhuǎn)換到z平面，得到數(shù)字低通濾波器的系統(tǒng)函數(shù)H(z)；如所設(shè)計(jì)的數(shù)字濾波器為高通、帶通或帶阻濾波器，則可借助模擬濾波器的頻帶變換轉(zhuǎn)換為低通模擬濾波器。 由于直接法設(shè)計(jì)巴特沃思濾波器相對復(fù)雜，在不需要任意幅頻特征的情況下，一般采用問接法，同時(shí)由于沖激響應(yīng)不變法和階躍響應(yīng)不變法，從s平面轉(zhuǎn)換到z平面的映射為多值映射，容易造成頻譜混疊，故而本文采用不會(huì)產(chǎn)生頻譜混疊的雙線性變換法。 4 數(shù)字巴特沃思濾波器的雙線性變換法設(shè)計(jì) 設(shè)計(jì)IIR濾波器的任務(wù)就是尋求一個(gè)因果、物理上可實(shí)現(xiàn)的系統(tǒng)函數(shù)H(z)使其頻率響應(yīng)滿足所希望得到的頻域指標(biāo)，即符合給定的通帶截止頻率、 阻帶截止、通帶衰減和阻帶衰減等指標(biāo)。 雙線性變換法的基本步驟如下：(1)根據(jù)任務(wù)，確定性能指標(biāo)：在設(shè)計(jì)之前，首根據(jù)工程實(shí)際的需要確定巴特沃思濾波器的技術(shù)指標(biāo)：邊界頻率：，；通帶最大衰減和阻帶最小衰減； (2)將數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換成模擬濾波器指標(biāo)。雙線性變換法的轉(zhuǎn)換關(guān)系為：，為數(shù)字角頻率，為模擬角頻率，為取樣周期； (3)用模擬濾波器設(shè)計(jì)方法得到模擬巴特沃斯濾波器的傳輸函數(shù)H(s)，這可借助巴特沃思濾波器相關(guān)的嚴(yán)格的設(shè)計(jì)公式、現(xiàn)成的曲線和圖表得到；(4)映射實(shí)現(xiàn)：利用雙線性變換法關(guān)系式 將模擬巴特沃思濾波器H(s)轉(zhuǎn)換成數(shù)字巴特沃思濾波器H(z)； (5)用有限精度算法實(shí)現(xiàn)這個(gè)系統(tǒng)函數(shù)H(z)，同時(shí)要確定該濾波器的結(jié)構(gòu)， 因?yàn)槿绻谟布蠈?shí)現(xiàn)該濾波器必須考慮計(jì)算的復(fù)雜度、存儲(chǔ)容量限制等因素， 而用Matlab實(shí)現(xiàn)無此影響； (6)用合適的軟、硬件技術(shù)實(shí)現(xiàn)該濾波器。 5 數(shù)字巴特沃思濾波器的Matlab實(shí)現(xiàn)按照傳統(tǒng)的方法設(shè)計(jì)巴特沃思數(shù)字濾波器，不論采用采用何種方法都面臨著繁雜的數(shù)學(xué)計(jì)算或反復(fù)的累加試驗(yàn)，以驗(yàn)證設(shè)計(jì)結(jié)果是否滿足指標(biāo)要求，工作量巨大、且效率極低。而現(xiàn)代分析信號(hào)時(shí)，使用軟件Matlab是一種高效、簡便的設(shè)計(jì)方法。 Matlab是一款高性能數(shù)值計(jì)算和可視化軟件，應(yīng)用及其廣泛。Matlab軟件的誕生，尤其是數(shù)字信號(hào)處理工具箱(Signal Processing Toolbox) 和濾波器設(shè)計(jì)工具箱(Filter Design Toolbox)推出后，使數(shù)字信號(hào)處理系統(tǒng)的分析與設(shè)計(jì)問題變得簡單了,它為數(shù)字信號(hào)處理系統(tǒng)的設(shè)計(jì)與仿真提供了一個(gè)強(qiáng)有力的工具。在設(shè)計(jì)數(shù)字濾波器時(shí)，通常采用Matlab來進(jìn)行輔助設(shè)計(jì)和仿真。可方便地分析顯示濾波器系統(tǒng)的幅頻、相頻特性,以及對給定信號(hào)經(jīng)所設(shè)計(jì)的濾波器進(jìn)行濾波前和濾波后的幅頻特性。 由于沖擊響應(yīng)不變法和階躍響應(yīng)不變法的不足，下面以雙線性變換法設(shè)計(jì)一個(gè)巴特沃思低通濾波器，下面以雙線性變換法設(shè)計(jì)一個(gè)巴特沃思低通濾波器。這可以嚴(yán)格按照上述步驟進(jìn)行，并且有很多文獻(xiàn)。正是按照這個(gè)步驟進(jìn)行，但通過仔細(xì)研究Matlab軟件自帶的信號(hào)處理工具箱下Matlabtoolboxsignalsignalbutter.m和Matlabtoolboxsignalsignalbuttord.m 等文件，可以清楚看到：Matlab默認(rèn)采用雙線性變換法設(shè)計(jì)，因而無需完全嚴(yán)格按照上述步驟進(jìn)行，完全可以采用更簡潔的方法。所以設(shè)計(jì)濾波器可以運(yùn)用Matalb中的程序?qū)崿F(xiàn)也可以通過它中的Smulink仿真模塊來搭建系統(tǒng)實(shí)現(xiàn)，具體的做法是：運(yùn)行Matalb軟件，在路徑SmulinkSignal Processing Blockset Filtering Filter Designs Analog Filter Design，可以設(shè)置highpass,lowpass,bandpass等。如下圖所示。接著就可以看到濾波器設(shè)計(jì)方法，濾波器類型，階數(shù)及邊界頻率的參數(shù)，如下圖所示。本文采用的是數(shù)字濾波器的設(shè)計(jì)，所以選擇路徑SmulinkSignal Processing Blockset Filtering Filter Designs Digtal Filter Design或者直接在Matlab命令窗口中直接敲入FDATOOL運(yùn)行即可，即通過FDATool界面設(shè)計(jì)，則可以相應(yīng)的設(shè)置濾波器的類型及參數(shù)。本文采用FDATool界面設(shè)計(jì)濾波器。下面詳細(xì)介紹一下FDATool界面。如下圖所示。FDATool（Filter Design & Analysis Tool）是Matlab信號(hào)處理工具箱里專用的濾波器設(shè)計(jì)分析工具，Matlab6.0以上的版本還專門增加了濾波器設(shè)計(jì)工具箱（Filter Design Toolbox）。FDATool可以設(shè)計(jì)幾乎所有的基本的常規(guī)濾波器，包括FIR和IIR的各種設(shè)計(jì)方法。它操作簡單，方便靈活。 Filter Type（濾波器類型）選項(xiàng)，包括Lowpass（低通）、Highpass（高通）、Bandpass（帶通）、Bandstop（帶阻）和特殊的FIR濾波器。Design Method（設(shè)計(jì)方法）選項(xiàng)，包括IIR濾波器的Butterworth（巴特沃思）法、Chebyshev Type I（切比雪夫I型）法、Chebyshev Type II（切比雪夫II型） 法、Elliptic（橢圓濾波器）法和FIR濾波器的Equiripple法、Least-Squares（最小乘方）法、Window（窗函數(shù)）法。Filter Order（濾波器階數(shù)）選項(xiàng)，定義濾波器的階數(shù)，包括Specify Order（指定階數(shù)）和Minimum Order（最小階數(shù)）。在Specify Order中填入所要設(shè)計(jì)的濾波器的階數(shù)（N階濾波器，Specify OrderN-1），如果選擇Minimum Order則Matlab根據(jù)所選擇的濾波器類型自動(dòng)使用最小階數(shù)。Frenquency Specifications選項(xiàng)，可以詳細(xì)定義頻帶的各參數(shù)，包括采樣頻率Fs和頻帶的截止頻率。它的具體選項(xiàng)由Filter Type選項(xiàng)和Design Method選項(xiàng)決定，例如Bandpass（帶通）濾波器需要定義Fstop1（下阻帶截止頻率）、Fpass1（通帶下限截止頻率）、Fpass2（通帶上限截止頻率）、Fstop2（上阻帶截止頻率），而Lowpass（低通）濾波器只需要定義Fstop1、Fpass1。采用窗函數(shù)設(shè)計(jì)濾波器時(shí)，由于過渡帶是由窗函數(shù)的類型和階數(shù)所決定的，所以只需要定義通帶截止頻率，而不必定義阻帶參數(shù)。Magnitude Specifications選項(xiàng)，可以定義幅值衰減的情況。例如設(shè)計(jì)帶通濾波器時(shí)，可以定義Wstop1（頻率Fstop1處的幅值衰減）、Wpass（通帶范圍內(nèi)的幅值衰減）、Wstop2（頻率Fstop2處的幅值衰減）。當(dāng)采用窗函數(shù)設(shè)計(jì)時(shí)，通帶截止頻率處的幅值衰減固定為6dB，所以不必定義。Window Specifications選項(xiàng)，當(dāng)選取采用窗函數(shù)設(shè)計(jì)時(shí)，該選項(xiàng)可定義，它包含了各種窗函數(shù)。在設(shè)計(jì)過程中，可以對比濾波器幅頻相頻特性和設(shè)計(jì)要求，隨時(shí)調(diào)整參數(shù)和濾波器類型，以便得到最佳效果。最后設(shè)置參數(shù)，設(shè)計(jì)本文所需濾波器，參數(shù)選擇如下圖所示。如果問題中提供了被干擾的ECG信號(hào)數(shù)據(jù)，就可以用設(shè)計(jì)的濾波器通過Simulink搭建模塊來檢驗(yàn)它濾波的特性。6 總結(jié)利用Matlab的強(qiáng)大運(yùn)算功能，基于Matlab信號(hào)處理工具箱（Signal Processing Toolbox）的數(shù)字濾波器設(shè)計(jì)法可以快速有效的設(shè)計(jì)由軟件組成的常規(guī)數(shù)字濾波