數(shù)字濾波器設計及在心電信號濾波中應用

1、課 程 設 計 報 告 課程名稱 數(shù)字信號處理課程設計 課題名稱 數(shù)字濾波器設計及在心電信號濾波中的應用 專 業(yè) 通信工程 班 級 學 號 姓 名 指導教師 彭 禎 張鏖烽 郭 芳 2012年 9 月 10日湖南工程學院課 程 設 計 任 務 書課程名稱 數(shù)字信號處理課程設計 課題名稱 數(shù)字濾波器設計及在心電信號濾波中的應用 專 業(yè) 通信工程 班 級 學 號 姓 名 指導教師 彭 禎 張鏖烽 郭 芳 審 批 任務書下達日期 2012 年 9月 1 日任務完成日期 2012 年 9月 10日一、 課程設計的性質與目的數(shù)字信號處理課程是通信專業(yè)的一門重要專業(yè)基礎課，是信息的數(shù)字化處理、存儲和應用的

2、基礎。通過該課程的課程設計實踐，使學生對信號與信息的采集、處理、傳輸、顯示、存儲、分析和應用等有一個系統(tǒng)的掌握和理解；鞏固和運用在數(shù)字信號處理課程中所學的理論知識和實驗技能，掌握數(shù)字信號處理的基礎理論和處理方法，提高分析和解決信號與信息處理相關問題的能力，為以后的工作和學習打下基礎。數(shù)字濾波器是一種用來過濾時間離散信號的數(shù)字系統(tǒng)，通過對抽樣數(shù)據(jù)進行數(shù)學處理來達到頻域濾波的目的。根據(jù)其單位沖激響應函數(shù)的時域特性可分為兩類：無限沖激響應（IIR）濾波器和有限沖激響應（FIR）濾波器。二、課程設計題目數(shù)字濾波器設計及在心電信號濾波中的應用，其內容如下：1、心電信號采集心電信號作為心臟電活動在人體體表

3、的表現(xiàn)，信號一般比較微弱，幅度在10V5mV，頻率為0.05100Hz。在心電信號的采集、放大、檢測及記錄過程中，有來自外界的各種干擾。記錄一段時間內的人體心電信號波形，要求長度不小于10秒，并對記錄的信號進行數(shù)字化，保存為數(shù)據(jù)文件；這里，請同學們使用美國的MIT/BIH心電原始數(shù)據(jù)，由實驗老師給出一定長度的的心電原始數(shù)據(jù)，數(shù)據(jù)保存在文件“a01.txta10.txt”中，在MATLAB中通過如下語句讀?。?從當前路徑下的a01.txt文件讀取心電原始數(shù)據(jù)到變量a01中，a01為二維數(shù)據(jù)，第一列%為心電信號時間，第二列為心電信號幅度。2、心電信號分析使用MATLAB繪出數(shù)字化后的心電信號的時域

4、波形和頻譜圖。根據(jù)頻譜圖求出其帶寬，并說明心電信號的基本特征。3、含噪心電信號合成在MATLAB軟件平臺下，給原始的心電信號疊加上噪聲或干擾，干擾類型分為如下幾種：（1）白噪聲；（2）工頻干擾（50Hz）；（3）諧波干擾（二次、三次諧波為主，分別為100Hz、150Hz）；（4）其它干擾，可設置為低頻、高頻、帶限噪聲，或沖激干擾。繪出疊加噪聲后的心電信號時域和頻譜圖，在視覺上與原始心電信號圖形對比，繪出其時域波形差，分析頻域基本特征變化。4、數(shù)字濾波器設計及濾波給定濾波器的規(guī)一化性能指標（參考指標,實際中依據(jù)每個同學所疊加噪聲情況而定）例如：通帶截止頻率wp=0.25*pi, 阻通帶截止頻率w

5、s=0.3*pi; 通帶最大衰減Rp=1 dB; 阻帶最小衰減Rs=15 dB，每個題目至少設計出5個用不同方法的不同類型濾波器。采用雙線性變換法與脈沖響應不變法，分別利用不同的原型低通濾波器（Butterworth型與切比雪夫I型）來設計各型IIR濾波器（低通、高通、帶通、帶阻中的至少3種類型），繪出濾波器的頻域響應；并用這些數(shù)字濾波器對含噪心電信號分別進行濾波處理，比較不同方法下設計出來的數(shù)字濾波器的濾波效果，并從理論上進行分析（或解釋）。5、心電信號波形觀察、頻譜觀察對濾波后的心電信號觀察其時域、頻域特征變化。繪出濾波后、濾波前、加噪后三個心電信號的差值波形，觀察相互間的差異性；同時，分

6、析頻譜變化。學生也可選用持續(xù)時間更長的心電原始數(shù)據(jù)，加上干擾后按上述要求設計濾波器。三、課程設計要求1、在一周內學生須上機16小時以上，程序調試完后，須由指導老師在機器上檢查運行結果，經教師認可后的源程序可通過打印機輸出，并請教師在程序清單上簽字。2、課程設計報告內容和格式：設計題目，設計的詳細步驟，設計過程中的結果、圖形等，設計總結。3、每組每人必須獨立完成，成績的考核按設計結果、答辯成績及課程設計報告來綜合評定。成績分為優(yōu)、良、中、及格、不及格五級分評定。4、指導教師：彭禎，張鏖烽，郭芳。四、設計進度安排通信工程1001/1002:1周周一上午，E-412, 任務講解與布置，學生分組選題，

7、查找相關資料，準備課程設計，學生上機，按任務要求進行課程設計；分組選題；1周周二上午，E-412,學生上機，按任務要求進行課程設計；分組選題；1周周三上午，E-412,學生上機，按任務要求進行課程設計；分組選題；1周周四上午，E-412,學生任務完成，答辯并提交課程設計報告。附：課程設計報告裝訂順序：封面、任務書、目錄、正文、評分、附件（A4大小的圖紙及程序清單）。 正文的格式:一級標題用3號黑體,二級標題用四號宋體加粗,正文用小四號宋體;行距為22。正文的內容:一、課題的主要功能；二、課題的功能模塊的劃分；三、主要功能的實現(xiàn)；四、程序調試；五、總結；六、附件（所有程序的原代碼，要求對程序寫出

8、必要的注釋）；七、評分表。目 錄1 課程設計的目的12 課程設計的原理12.1 Butterworth低通數(shù)字濾波器的設計12.2切比雪夫I型數(shù)字低通濾波器2 2.3 IIR數(shù)字濾波器的性質.23 課程設計設計步驟及結果分析43.1心電數(shù)據(jù)的導入43.2繪出心電信號的時域圖和頻譜圖43.3加入噪聲干擾53.3.1白噪聲.53.3.2工頻干擾（50Hz）.73.3.3帶限chirp噪聲.83.4濾波器的設計93.4.1 Butterwort型低通數(shù)字濾波器.93.4.2切比雪夫I型數(shù)字低通濾波器.123.4.3帶阻濾波器的設計.154 心得體會171 課程設計的目的數(shù)字濾波器是指輸入,輸出均為數(shù)

9、字信號，通過數(shù)值運算處理改變輸入信號所含頻率成分的相對比例，或者濾除某些頻率成分的數(shù)字器件或程序。因此，數(shù)字濾波的概念和模擬濾波相同，只是信號的形式和實現(xiàn)濾波方法不同。正因為數(shù)字濾波通過數(shù)值運算實現(xiàn)濾波，所以數(shù)字濾波器處理精度高，穩(wěn)定，體積小，重量輕，靈活，不存在阻抗匹配問題，可以實現(xiàn)模擬濾波器無法實現(xiàn)的特殊濾波功能。希望學生運用數(shù)字信號處理課程中所學的理論知識和實驗技能，基本掌握數(shù)字信號處理的基礎理論和處理方法，提高分析和解決信號與信息處理相關問題的能力，為以后的工作和學習打下基礎。2 課程設計的原理2.1 Butterworth低通數(shù)字濾波器的設計巴特沃斯低通濾波器的平方幅度響應為其中，n

10、為濾波器的階數(shù)，為低通濾波器的截止頻率。 該濾波器具有 一些特殊的性質： 對所有的n,都有當時, ； 對所有的n,都有當時, ；是的單調遞減函數(shù),即不會出現(xiàn)幅度響應的起伏； 當時,巴特沃斯濾波器趨向于理想的低通濾波器； 在處平方幅度響應的各級導數(shù)均存在且等于0，因此在該點上取得最大值，且具有最大平坦特性。 圖1展示了2階、4階、8階巴特沃斯低通濾波器的幅頻特性。可見階數(shù)n越高，其幅頻特性越好，低頻檢測信號保真度越高，過渡帶變窄，即衰減加劇，但半功率點不變。圖1 巴特沃斯低通濾波器的幅頻特性2.2 切比雪夫I型數(shù)字低通濾波器（1）確定數(shù)字低通濾波器的技術指標：通帶截止頻率p、通帶衰減ap、阻帶截

11、止頻率s、阻帶衰減as切比雪夫濾波器的振幅平方特性如圖2所示：圖2 切比雪夫濾波器的振幅平方特性(2)將數(shù)字低通濾波器的技術指標轉換成模擬低通濾波器的技術指標。如果采用脈沖響不變法，邊界頻率的轉換關系為:如果采用雙線性變換法，邊界頻率的轉換關系為(3) 按照模擬低通濾波器的技術指標設計模擬低通濾波器。(4) 利用雙線性變換法將模擬濾波器Ha(s)，從s平面轉換到z平面，得到數(shù)字低通濾波器系統(tǒng)函數(shù)H(z)。（5）數(shù)字低通技術指標為：p=0.4rad, a p=1dB; s=0.5rad, a s=40Db（6）模擬低通的技術指標為：歸一化截止角頻率wp=2pi*Fs/Ft; ws=2pi*Fs/

12、Ft(7)利用模擬切比雪夫濾波器設計數(shù)字濾波器。通帶截止頻率為：wp=0.4*pi; 阻帶截止頻率為：ws=0.5*pi;通帶最大衰減為：Rp=1;阻帶最大衰減為：As=15;設定周期為1s；模擬低通濾波器的生成：b,a=cheby1(n,1,Wn,'low','s');滿足設計指標的最小階數(shù)和截止頻率：Wnn,Wn=cheb1ord(OmegaP,OmegaS,1,40,'s')。最后實現(xiàn)輸入輸出、幅頻特性、相頻特性的圖形。2.3 IIR數(shù)字濾波器的性質無限長沖激響應( IIR) 數(shù)字濾波器是數(shù)字濾器的一種，數(shù)字濾波器根據(jù)其沖激響應函數(shù)的時域特

13、性,還包括有限長沖激響應( FIR)數(shù)字濾波器。IIR數(shù)字濾波器的特征是：具有無限持續(xù)時間沖激響應；需要用遞歸模型來實現(xiàn)，這可以從其差分方程得出，也可以從其系數(shù)函數(shù)為：得出。 數(shù)字巴特沃思濾波器屬于IIR濾波器，該類濾波器具有特定的性質和設計方法。目前比較成熟的IIR數(shù)字濾波器設計方法有兩種：1）直接法目前所用的方法主要是：零極點累試法、頻域幅度平方誤差最小法和時域單位脈沖響應逼近法。直接法的最大優(yōu)點在于可以設計任意幅頻特征的濾波器。2)間接法，目前所用的方法主要是：沖激響應不變法、階躍響應不變法和雙線性法。它們都是借助于 已經成熟的現(xiàn)有低通濾波器原型進行設計，即對數(shù)字低通數(shù)字濾波器，先將數(shù)字

14、低通濾波器的技術指標按希望的設計方法轉換為模擬低通濾波器的技術指標，再按指定的模擬低通濾波器的類型設計模擬濾波器H(s)，然后，將模擬濾波器的系統(tǒng)函數(shù)H(s)從s平面轉換到z平面，得到數(shù)字低通濾波器的系統(tǒng)函數(shù)H(z)；如所設計的數(shù)字濾波器為高通、帶通或帶阻濾波器，則可借助模擬濾波器的頻帶變換轉換為低通模擬濾波器。 由于直接法設計巴特沃思濾波器相對復雜，在不需要任意幅頻特征的情況下，一般采用問接法，同時由于沖激響應不變法和脈沖響應不變法，從s平面轉換到z平面的映射為多值映射，容易造成頻譜混疊，故而本文采用不會產生頻譜混疊的雙線性變換法。 3 課程設計步驟及結果分析3.1心電數(shù)據(jù)的導入將老師給的心

15、電信號原始數(shù)據(jù)存于桌面，然后我選用桌面上的第一組數(shù)據(jù)a=load('C:Documents and Settingshnie桌面心電信號數(shù)據(jù)a01.txt');3.2 繪出心電信號的時域圖和頻譜圖將導入的兩行數(shù)據(jù)分別用t,b來替換，然后通過調用plot函數(shù)來畫出時域圖，然后通過對1000個心電數(shù)據(jù)的幅值進行FFT運算，再次調用plot函數(shù)來繪出頻域圖，具體設計如下：figure(1);subplot(2,1,1);t=a(1:1000,1);b=a(1:1000,2);plot(t,b);title('原始波形圖');xlabel('時間（s）'

16、);ylabel('幅值（A）');y1=fft(a(:,2),1000);f1=100*(0:999)/1000;subplot(2,1,2);plot(f1,abs(y1);title('原始頻譜圖');xlabel('頻率（Hz）');ylabel('幅度(dB)');圖圖1.1 原始時域和頻譜圖通過導入的心電信號數(shù)據(jù)發(fā)現(xiàn)在其頻譜圖上的020Hz和80100Hz之間的幅值比較大，而在3070Hz之間的幅值相對較小。3.3 加入噪聲干擾這里我加入的噪聲是：白噪聲，50Hz的工頻噪聲，帶限chirp噪聲。3.3.1 白噪聲通過用

17、s來代表加入白噪聲后的信號，并進行數(shù)字濾波器的頻率響應，對s中的1000個頻率點調用plot函數(shù)畫出加入白噪聲后的時域圖，再對ws/pi,abshs調用plot函數(shù)畫出加入白噪聲后的頻譜圖，具體操作如下：q=0.8*rand(1000,1);s=a(:,2)+q;hs,ws=freqz(s,1,1024);abshs=abs(hs);figure(2);subplot(2,1,1);plot(s(1:1000);title('加入白噪聲后的時域圖');xlabel(時間（s）');ylabel('幅值（A）');subplot(2,1,2);plot(w

18、s/pi,abshs);title('加入白噪聲后的頻譜圖');xlabel('幅度');ylabel('Hz');圖1.2 加入白噪聲后的時域圖和頻譜圖通過觀察加入白噪聲后的時域圖和頻域圖，將它與未加入白噪聲進行比較，可以發(fā)現(xiàn)頻譜圖在0Hz時的幅度增加的很大，而且又在沒有譜線的頻率上竟然出現(xiàn)了頻譜，這是由于白噪聲在所有頻率上都有頻率造成的。3.3.2 工頻干擾（50Hz）用x表示加入工頻干擾后的信號，再對工頻信號的1000個頻率點進行FFT運算，在進行相關的運算后，通過調用plot函數(shù)直接繪出加入工頻干擾后的時域圖和頻譜圖。具體步驟如下：x2=

19、sin(2*pi*50*t); t=0:0.001:0.001*(1000-1);x1=a(:,2);x=x1+x2;y2=fft(x2,1000);f2=100*(0:999)/1000;figure(3);subplot(2,1,1);plot(t,x);title('加入工頻干擾后的時域圖');xlabel('時間（s）');ylabel('幅值（A）');subplot(2,1,2);plot(f2,abs(y2);title('加入工頻干擾后的頻譜圖');xlabel('幅度');ylabel('

20、Hz');圖1.3 加入工頻干擾后的時域圖和頻譜圖比較圖1.1和圖1.3，發(fā)現(xiàn)加噪聲后的幅值有比較微小的變化，而在頻譜圖在38Hz和63Hz附近增值很大，這是由于50Hz的工頻噪聲造成的，而在其他頻率范圍內也有比較明顯的變化。3.3.3 帶限chirp噪聲 同樣用k表示加入帶限chirp噪聲后的信號，再用freqz進行濾波器的頻率響應，調用plot對k畫出加入帶限噪聲后的時域和頻域圖，具體描述如下：p=0.5*chirp(a(:,1),0,a(1000,1),200);k=x1+p;hc,wc=freqz(k,1,1024);abshc=abs(hc);figure(4);subplo

21、t(2,1,1);plot(k(1:1000);title('加入帶限噪聲后的時域圖');xlabel('時間（s）');ylabel('幅值（A）');subplot(2,1,2);plot(wc/pi,abshc);title('加入帶限噪聲后的頻譜圖');xlabel('幅度');ylabel('Hz');圖1.4 加入帶限噪聲后的時域圖和頻譜圖加入帶限chirp噪聲干擾后，發(fā)現(xiàn)加噪之前時域圖是均勻分布的，而加噪后，則變?yōu)榍懊嫦∈韬竺婷芗那闆r了，對于頻譜圖而言，加噪前是中間凹兩邊凸，加噪后是

22、前凸后平，在個別點上幅值增加很多，這是由于加入chirp噪聲的結果。3.4 濾波器的設計3.4.1 Butterworth型低通數(shù)字濾波器用wp和ws表示分別將通帶，阻帶截止頻率的角頻率表示，在分別計算階數(shù)n1和截止頻率Wn,再設計低通Butterworth型模擬濾波器，然后采用雙線性法將模擬濾波器系數(shù)變?yōu)閿?shù)字濾波器系數(shù)，畫出濾波器頻譜圖，調用filter實現(xiàn)對工頻干擾的濾波，用plot函數(shù)畫出濾除工頻干擾后的時域圖和濾除白噪聲后的頻譜圖。具體操作如下：figure(5);fs=100;f1=5;f2=10;wp=(f1/fs)*2*pi;ws=(f2/fs)*2*pi;Omegap=2*fs

23、*tan(wp/2);Omegas=2*fs*tan(ws/2);n1,Wn=buttord(Omegap,Omegas,1,50,'s');b,a=butter(n1,Wn,'s');bz,az=bilinear(b,a,fs);freqz(bz,az,512,fs);y=filter(bz,az,x2);figure(6);subplot(2,1,1)plot(y);title('濾除工頻干擾后的時域圖');xlabel('時間（s）');ylabel('幅值（A）');y3=fft(y,1000);f1=10

24、0*(0:999)/1000;subplot(2,1,2);plot(f1,abs(y3);title('濾除白噪聲后的頻譜圖');xlabel('頻率(Hz)');ylabel('幅值(dB)');圖1.5 Butterworth型低通數(shù)字濾波器頻率響應圖圖1.6 濾除單頻正弦波后的時域圖和濾除噪聲后的頻譜圖由于信號處于頻段的低頻部分，而工頻信號的頻譜在整個上是呈對稱分布的，與源信號的頻段分布是相似的，采用低通濾波器將噪聲信號的高頻部分濾掉，由圖1.3和1.6會發(fā)現(xiàn)加入工頻干擾后的源信號，濾波后時域圖的幅值比原來變小了，濾波后的頻譜圖頻率在01

25、0Hz和90100Hz之間的幅值變化較大，而在1090Hz的區(qū)間，頻譜圖基本趨于直線。3.4.2切比雪夫I型數(shù)字低通濾波器用Wp1,Wp2,Ws1,Ws2表示分別將通帶，阻帶截止頻率的角頻率表示,算出頻帶寬帶，計算階數(shù)n1和截止頻率WN，再設計切比雪夫I型模擬濾波器，采用雙線性法將模擬濾波器系數(shù)變?yōu)閿?shù)字濾波器系數(shù)，畫出切比雪夫I型數(shù)字濾波器的頻率響應，調用filter實現(xiàn)對白噪聲的濾波，再最后調用plot函數(shù)畫出濾除白噪聲后的時域圖和頻域圖。具體過程如下：figure(7);fs=100;f11=10;f12=25;f21=5;f22=30;Wp1=(f11/fs)*2*pi;Ws1=(f21

26、/fs)*2*pi;Wp2=(f12/fs)*2*pi;Ws2=(f22/fs)*2*pi;Omegap1=2*fs*tan(Wp1/2);Omegap2=2*fs*tan(Wp2/2);Omegas1=2*fs*tan(Ws1/2);Omegas2=2*fs*tan(Wp2/2);BW=Omegap2-Omegap1;W0=Omegap1*Omegap2;W00=sqrt(W0);WP=1;WS=WP*(W02-Ws12)/(Ws1*BW);n1,WN=buttord(WP,WS,1,50,'s');B,A=cheby1(n1,1,WN,'s');BT,AT=

27、lp2bp(B,A,W00,BW);num,den=bilinear(BT,AT,0.5);freqz(num,den,64);y=filter(num,den,s);figure(8);subplot(2,1,1);plot(y);title('濾除白噪聲后的時域圖');xlabel('時間（s）');ylabel('幅值（A）');s3=fft(y,1000);f1=100*(0:999)/1000;subplot(2,1,2);plot(f1,abs(s3);title('濾除白噪聲后的頻譜圖');xlabel('頻率(Hz)');ylabel('幅值(dB)');圖1.7 切比雪夫1型數(shù)字濾波器的頻率響應圖圖1.8 濾除白噪聲后的時域圖和頻域圖由于白噪聲在整個頻段上都存在，對與源信號共存的低頻信號用選頻濾波器是無法濾除的，故而采用低通濾波器將噪聲的高頻部分去掉。比較圖1.2和圖1.7，它們在時域圖上差不多是一致的，但在頻譜圖中的圖1.2在0Hz上的幅值是很搞的，在濾波后的就明顯變小了。3.4.3 帶阻濾波器的設計 設計帶阻濾波器