濾波器設(shè)計(jì)及在心電信號(hào)濾波中應(yīng)用

1、西南科技大學(xué)課 程 設(shè) 計(jì) 報(bào) 告課程名稱： 通信系統(tǒng)課程設(shè)計(jì) 設(shè)計(jì)名稱：濾波器設(shè)計(jì)及在心電信號(hào)濾波中的應(yīng)用姓 名： 學(xué) 號(hào): 班 級(jí)： 指導(dǎo)教師： 起止日期： 2013.6.24-2013.7.6 西南科技大學(xué)信息工程學(xué)院制課 程 設(shè) 計(jì) 任 務(wù) 書(shū)學(xué)生班級(jí)： 學(xué)生姓名： 學(xué)號(hào)： 設(shè)計(jì)名稱： 濾波器設(shè)計(jì)及在心電信號(hào)濾波中的應(yīng)用 起止日期： 2013.6.24-2013.7.6 指導(dǎo)教師： 設(shè)計(jì)要求：一般正常人的心電信號(hào)頻率在0.7100HZ范圍內(nèi)，幅度為10v（胎兒）5mv（成人）。人體心電信號(hào)微弱，信噪比小，因此，在采集心電信號(hào)時(shí)，易受到儀器、人體活動(dòng)等因素的影響，而且所采集的心電信號(hào)常伴

2、有干擾。 采集心電數(shù)據(jù)時(shí)，由于人的說(shuō)話呼吸，常常會(huì)混有約為0.1Hz到0.25Hz頻段的干擾，對(duì)于這些低頻干擾，可以讓信號(hào)通過(guò)一個(gè)高頻濾波器，低截止頻率設(shè)置為0.25，來(lái)濾波低頻信號(hào)； 對(duì)于高頻信號(hào)干擾，可以讓信號(hào)再通過(guò)一個(gè)低頻濾波器，其中截止頻率設(shè)置為99Hz；由于電子設(shè)備采集到的信號(hào)經(jīng)常會(huì)混有電源線干擾，電源線干擾是以50 Hz為中心的窄帶噪聲，帶寬小于1Hz。對(duì)于電源線干擾，可以讓信號(hào)通過(guò)一個(gè)帶阻濾波器。 某一個(gè)實(shí)際心電圖信號(hào)的采樣序列存于文件ecg.txt中，共4000個(gè)采樣點(diǎn)。要求根據(jù)以上指標(biāo)進(jìn)行如下設(shè)計(jì)： 1、設(shè)計(jì)一個(gè)低通濾波器濾除心電信號(hào)中的高頻信號(hào)干擾，對(duì)濾波前后的心電信號(hào)的頻

3、譜進(jìn)行分析比較。 2、設(shè)計(jì)一個(gè)高通濾波器濾除心電信號(hào)中的低頻干擾，對(duì)濾波前后的心電信號(hào)的頻譜進(jìn)行分析比較。 3、設(shè)計(jì)一個(gè)帶阻濾波器（50Hz陷波器）濾除心電信號(hào)中的電源線干擾，對(duì)濾波前后的心電信號(hào)的頻譜進(jìn)行分析比較。說(shuō)明：1、使用MATLAB或其它編程軟件編程實(shí)現(xiàn)。2、如果濾波效果不明顯，則需變動(dòng)濾波器參數(shù)指標(biāo)，重新設(shè)計(jì)濾波器，通過(guò)頻譜分析，多次試驗(yàn)確定最合適的濾波器。課 程 設(shè) 計(jì) 學(xué) 生 日 志時(shí)間設(shè)計(jì)內(nèi)容2013.6.24-2013.6.30查閱資料，確定初步方案2013.6.31-2013.7.2編寫代碼，調(diào)試改錯(cuò)2013.7.3-2013.7.5完善代碼并撰寫課程設(shè)計(jì)報(bào)告 2013.

4、7.6答辯課 程 設(shè) 計(jì) 評(píng) 語(yǔ) 表指導(dǎo)教師評(píng)語(yǔ)： 成績(jī)： 指導(dǎo)教師： 年 月 日濾波器設(shè)計(jì)及在心電信號(hào)濾波中的應(yīng)用一 設(shè)計(jì)目的和意義數(shù)字濾波器是指輸入,輸出均為數(shù)字信號(hào)，通過(guò)數(shù)值運(yùn)算處理改變輸入信號(hào)所含頻率成分的相對(duì)比例，或者濾除某些頻率成分的數(shù)字器件或程序。因此，數(shù)字濾波的概念和模擬濾波相同，只是信號(hào)的形式和實(shí)現(xiàn)濾波方法不同。正因?yàn)閿?shù)字濾波通過(guò)數(shù)值運(yùn)算實(shí)現(xiàn)濾波，所以數(shù)字濾波器處理精度高，穩(wěn)定，體積小，重量輕，靈活，不存在阻抗匹配問(wèn)題，可以實(shí)現(xiàn)模擬濾波器無(wú)法實(shí)現(xiàn)的特殊濾波功能。希望學(xué)生運(yùn)用數(shù)字信號(hào)處理課程中所學(xué)的理論知識(shí)和實(shí)驗(yàn)技能，基本掌握數(shù)字信號(hào)處理的基礎(chǔ)理論和處理方法，提高分析和解決信號(hào)

5、與信息處理相關(guān)問(wèn)題的能力，為以后的工作和學(xué)習(xí)打下基礎(chǔ)。二 設(shè)計(jì)原理 2.1 Butterworth低通數(shù)字濾波器的設(shè)計(jì)巴特沃斯低通濾波器的平方幅度響應(yīng)為其中，n為濾波器的階數(shù)，為低通濾波器的截止頻率。 該濾波器具有 一些特殊的性質(zhì)： 對(duì)所有的n,都有當(dāng)時(shí), ； 對(duì)所有的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階巴特沃斯低通濾波器的幅頻特性?？梢?jiàn)階數(shù)n越高，其幅頻特性越好，低頻檢測(cè)信號(hào)保真度越高，過(guò)渡帶變窄，即衰減加

6、劇，但半功率點(diǎn)不變。圖1 巴特沃斯低通濾波器的幅頻特性2.2 切比雪夫I型數(shù)字低通濾波器（1）確定數(shù)字低通濾波器的技術(shù)指標(biāo)：通帶截止頻率p、通帶衰減ap、阻帶截止頻率s、阻帶衰減as切比雪夫?yàn)V波器的振幅平方特性如圖2所示：圖2 切比雪夫?yàn)V波器的振幅平方特性(2)將數(shù)字低通濾波器的技術(shù)指標(biāo)轉(zhuǎn)換成模擬低通濾波器的技術(shù)指標(biāo)。如果采用脈沖響不變法，邊界頻率的轉(zhuǎn)換關(guān)系為:如果采用雙線性變換法，邊界頻率的轉(zhuǎn)換關(guān)系為(3) 按照模擬低通濾波器的技術(shù)指標(biāo)設(shè)計(jì)模擬低通濾波器。(4) 利用雙線性變換法將模擬濾波器Ha(s)，從s平面轉(zhuǎn)換到z平面，得到數(shù)字低通濾波器系統(tǒng)函數(shù)H(z)。（5）數(shù)字低通技術(shù)指標(biāo)為： p=

7、0.4rad，ap=1dB； s=0.5rad，as=40dB（6）模擬低通的技術(shù)指標(biāo)為： p=2Ttan12p，T=1 p=2tan0.1=0.65rad/s，p=1dB s=2tan0.15=1.019/s, s=15d歸一化截止角頻率wp=2pi*Fs/Ft; ws=2pi*Fs/Ft(7)利用模擬切比雪夫?yàn)V波器設(shè)計(jì)數(shù)字濾波器。通帶截止頻率為：wp=0.4*pi; 阻帶截止頻率為：ws=0.5*pi;通帶最大衰減為：Rp=1;阻帶最大衰減為：As=15;設(shè)定周期為1s；模擬低通濾波器的生成：b,a=cheby1(n,1,Wn,low,s);滿足設(shè)計(jì)指標(biāo)的最小階數(shù)和截止頻率：Wnn,Wn=

8、cheb1ord(OmegaP,OmegaS,1,40,s)。最后實(shí)現(xiàn)輸入輸出、幅頻特性、相頻特性的圖形。三 詳細(xì)設(shè)計(jì)步驟3.1心電數(shù)據(jù)的導(dǎo)入將老師給的心電信號(hào)原始數(shù)據(jù)存于桌面文件夾，導(dǎo)入數(shù)據(jù)代碼為：a=load(C:UsersAdministratorDesktop題目七數(shù)據(jù)文件題目七數(shù)據(jù)文件ecg.txt);3.2 繪出心電信號(hào)的時(shí)域圖和頻譜圖將導(dǎo)入的數(shù)據(jù)分別用t來(lái)替換，通過(guò)調(diào)用plot函數(shù)來(lái)畫(huà)出時(shí)域圖，然后通過(guò)對(duì)4000個(gè)心電數(shù)據(jù)的幅值進(jìn)行FFT運(yùn)算，再次調(diào)用plot函數(shù)來(lái)繪出頻域圖，具體設(shè)計(jì)如下：figure(1); %新建圖像subplot(2,1,1); %將2個(gè)圖畫(huà)到一個(gè)平面t=

9、a(1:4000,1); %用t替換導(dǎo)入的4000個(gè)數(shù)據(jù)plot(t); %繪出t的圖形title(原始波形圖);xlabel(時(shí)間（s）);ylabel(幅值（A）);y1=fft(a(:,1),4000); %行取全部f1=100*(0:3999)/4000; %先生成一個(gè)0,1,2,.,3999的整數(shù)向量，然后對(duì)對(duì)這個(gè)向量的每一項(xiàng)乘以100除以4000.subplot(2,1,2); %將2個(gè)圖畫(huà)到一個(gè)平面plot(f1,abs(y1); %f1為橫坐標(biāo)，abs(y1)為縱坐標(biāo)作圖title(原始頻譜圖);xlabel(頻率（Hz）);ylabel(幅度(dB);3.3 加入干擾這里我加

10、入的干擾是：白噪聲高頻干擾，50Hz的電源線干擾。3.3.1 白噪聲通過(guò)用s來(lái)代表加入白噪聲后的信號(hào)，并進(jìn)行數(shù)字濾波器的頻率響應(yīng)，對(duì)s中4000個(gè)頻率點(diǎn)調(diào)用plot函數(shù)畫(huà)出加入白噪聲后的時(shí)域圖，再對(duì)ws/pi,abshs調(diào)用plot函數(shù)畫(huà)出加入白噪聲后的頻譜圖，具體操作如下：q=50*rand(4000,1); %產(chǎn)生4000行1列的位于（0，1）區(qū)間的隨機(jī)數(shù)s=a(:,1)+q; %產(chǎn)生的隨機(jī)數(shù)與原始信號(hào)疊加賦給shs,ws=freqz(s,1,4096); %求離散系統(tǒng)頻響特性的函數(shù)freqz()abshs=abs(hs); %hs取絕對(duì)值賦給abshsfigure(2); %新建圖像2s

11、ubplot(2,1,1); plot(s(1:4000); %繪出從1取到4000時(shí)s的圖形title(加入高頻干擾信號(hào)后的時(shí)域圖);xlabel(時(shí)間（s）);ylabel(幅值（A）);subplot(2,1,2);plot(ws/pi,abshs); % ws/pi為橫坐標(biāo)abshs為縱坐標(biāo)作圖title(加入高頻干擾信號(hào)后的頻譜圖);xlabel(Hz);ylabel(幅度);通過(guò)觀察加入白噪聲后的時(shí)域圖和頻域圖，將它與未加入白噪聲進(jìn)行比較可以發(fā)現(xiàn)頻譜圖在0Hz時(shí)的幅度增加的很大，而且又在沒(méi)有譜線的頻率上竟然出現(xiàn)了頻譜，這是由于白噪聲在所有頻率上都有頻率造成的。3.3.2 電源線干擾

12、（50Hz）用x表示加入電源線干擾后的信號(hào)，再對(duì)電源線信號(hào)的4000個(gè)頻率點(diǎn)進(jìn)行FFT運(yùn)算，在進(jìn)行相關(guān)的運(yùn)算后，通過(guò)調(diào)用plot函數(shù)直接繪出加入電源線干擾后的時(shí)域圖和頻譜圖。具體步驟如下：x2=sin(2*pi*50*t); %x2表示正弦信號(hào)t=0:0.00025:0.00025*(4000-1); %從0開(kāi)始間隔0.00025取值到 0.00025*(4000-1)x1=a(:,1); %將矩陣a用x1表示x=x1+x2; %x1和x2疊加賦給xy2=fft(x,4000); %對(duì)x信號(hào)做快速傅里葉變換f2=100*(0:3999)/4000; %先生成一個(gè)0,1,2,.,3999的整數(shù)向

13、量，然后對(duì)對(duì)這個(gè)向量的每一項(xiàng)乘以100除以4000.figure(3); %新建圖像3subplot(2,1,1);plot(t,x); %t為橫坐標(biāo)，x為縱坐標(biāo)作圖title(加入電源線干擾后的時(shí)域圖);xlabel(時(shí)間（s）);ylabel(幅值（A）);subplot(2,1,2);plot(f2,abs(y2); %f2為橫坐標(biāo)，abs(y2)為縱坐標(biāo)作圖title(加入電源線干擾后的頻譜圖);xlabel(幅);ylabel(Hz);3.4 濾波器的設(shè)計(jì)3.4.1切比雪夫I型數(shù)字低通濾波器用Wp1,Wp2,Ws1,Ws2表示分別用通帶和阻帶截止頻率的角頻率算出頻帶寬帶，計(jì)算階數(shù)n1

14、和截止頻率WN，再設(shè)計(jì)切比雪夫I型模擬濾波器，采用雙線性法將模擬濾波器系數(shù)變?yōu)閿?shù)字濾波器系數(shù)，畫(huà)出切比雪夫I型數(shù)字濾波器的頻率響應(yīng)，調(diào)用filter實(shí)現(xiàn)對(duì)白噪聲的濾波，再最后調(diào)用plot函數(shù)畫(huà)出濾除白噪聲后的時(shí)域圖和頻域圖。具體過(guò)程如下：figure(4); %新建圖像4fs=1000; %采樣頻率f11=10; %通帶頻率f12=25; %阻帶頻率Wp1=(f11/fs)*2*pi; %通帶角頻率Wp2=(f12/fs)*2*pi; %阻帶角頻率Omegap1=2*fs*tan(Wp1/2); %數(shù)字轉(zhuǎn)化為模擬Omegap2=2*fs*tan(Wp2/2); %數(shù)字轉(zhuǎn)化為模擬BW=Omega

15、p2-Omegap1; %頻帶寬度W0=Omegap1*Omegap2; W00=sqrt(W0);WP=1;n1,WN=buttord(WP,WS,1,50,s); %參數(shù)WP和WS分別是通帶邊界頻率和阻帶邊界頻率，通帶最大衰減為1阻帶最小衰減50，返回的參數(shù)n1和WN分別為濾波器的階數(shù)和3dB截止頻率B,A=cheby1(n1,1,WN,s); %計(jì)算出階數(shù)為n1、截止頻率為WN、通帶波紋最大衰減為1的數(shù)字低通濾波器，它的返回值a、b分別表示數(shù)字低通濾波器的系統(tǒng)函數(shù)的分子和分母的多項(xiàng)式系數(shù)BT,AT=lp2bp(B,A,W00,BW); num,den=bilinear(BT,AT,0.5

16、); %實(shí)現(xiàn)模數(shù)的映射freqz(num,den,64); y=filter(num,den,s);figure(5);subplot(2,1,1);plot(y);title(濾除高頻信號(hào)后的時(shí)域圖);xlabel(時(shí)間（s）);ylabel(幅值（A）);s3=fft(y,4000);f1=100*(0:3999)/4000; subplot(2,1,2);plot(f1,abs(s3);title(濾除高頻信號(hào)后的頻譜圖);xlabel(頻率(Hz);ylabel(幅值(dB);3.4.2 高通濾波器的設(shè)計(jì) 高通濾波器的設(shè)計(jì)可通過(guò)模擬低通濾波器再經(jīng)頻率變換而實(shí)現(xiàn)。其中模擬低通濾波器可根據(jù)

17、已經(jīng)存在的典型濾波器，如巴特沃斯濾波器等逼近實(shí)現(xiàn)，而由低通到高通轉(zhuǎn)換理論依據(jù)在緒論部分已經(jīng)進(jìn)行了詳細(xì)的論證，又本設(shè)計(jì)基于MATLAB仿真軟件實(shí)現(xiàn)，可利用MATLAB信號(hào)處理工具箱提供的各種函數(shù)模型實(shí)現(xiàn)，可忽略其中的一些復(fù)雜的函數(shù)變換，從而簡(jiǎn)化理論設(shè)計(jì)和論證，具體如下：figure(6)n=0:0.01:2; for ii=1:4 %定義循環(huán)，產(chǎn)生不同階數(shù)的曲線 switch ii case 1,N=2; case 2,N=5; case 3,N=10; case 4,N=30; endz,p,k=buttap(N); %調(diào)用Butterworth模擬低通濾波器原型函數(shù)b,a=zp2tf(z,p

18、,k); %將零點(diǎn)極點(diǎn)增益形式轉(zhuǎn)換為傳遞函數(shù)形式H,w=freqs(b,a,n); %按n指定的頻率點(diǎn)給出頻率響應(yīng)magH2=(abs(H).2; subplot(2,1,1);hold on;plot(w,magH2); endxlabel(w/wc);ylabel(|H(jw)|2);title(Butterworth 模擬低通濾波器原型);text(1.5,0.18,n=2) %對(duì)不同曲線做標(biāo)記text(1.3,0.08,n=5)text(1.16,0.08,n=10)text(0.93,0.98,n=20)grid on; %模擬高通濾波器設(shè)計(jì)m=0:0.01:2;for ii=1:4

19、 switch ii case 1,N=2; case 2,N=5; case 3,N=10; case 4,N=30; endz,p,k=buttap(N); b,a=zp2tf(z,p,k); bt,at=lp2hp(b,a,0.25*2*pi); %由低通原型濾波器轉(zhuǎn)換為截止頻率為0.25Hz的高通濾波器Ht,w=freqs(bt,at,m); subplot(2,1,2);hold on; %hold on 是保存axes內(nèi)圖像用的,如果在新畫(huà)圖像之后不想覆蓋原圖像就要加上hold on這句話plot(w,abs(Ht);endtitle(模擬高通濾波器);xlabel(w/pi);y

20、label(|H(jw)|2); text(0.5,0.28,n=2) text(1.0,0.12,n=4)text(1.5,0.28,n=10)text(2.0,0.10,n=30)grid on; %模擬高通濾波器性能pha=angle(Ht); %輸出系統(tǒng)的相頻特性figure(7);subplot(2,1,1);plot(w,20.*log10(pha);grid;title(模擬高通濾波器相頻特性); xlabel(w/wc);ylabel(相位/dB); mag=abs(Ht); %輸出系統(tǒng)的幅頻特性subplot(2,1,2);plot(w,20*log10(mag);grid;

21、title( 模擬高通濾波器幅頻特性);xlabel(w/wc);ylabel(幅度/dB);3.4.3 帶阻濾波器的設(shè)計(jì) 設(shè)計(jì)帶阻濾波器，再畫(huà)出濾波器的頻率響應(yīng)圖，然后用filter實(shí)現(xiàn)對(duì)工頻干擾信號(hào)的濾波，調(diào)用plot函數(shù)分別畫(huà)出濾除工頻干擾后的時(shí)域圖和頻譜圖。具體操作如下：figure(8);d=fir1(4000,0.157 0.17,stop);%freqz(d,512);y4=filter(d,1,x2);figure(9);subplot(2,1,1);plot(y4);title(濾除電源線干擾后的時(shí)域圖);xlabel(時(shí)間（s）);ylabel(幅值（A）);subplot

22、(2,1,2);y5=fft(y4,4000);f1=100*(0:3999)/4000;plot(f1,abs(y5);title(濾除電源線干擾后的頻譜圖);xlabel(頻率(Hz);ylabel(幅值(dB);四 設(shè)計(jì)結(jié)果及分析4.1 心電數(shù)據(jù)的導(dǎo)入，心電信號(hào)的時(shí)域圖和頻譜圖如下：圖3 原始時(shí)域和頻譜圖通過(guò)導(dǎo)入的心電信號(hào)數(shù)據(jù)發(fā)現(xiàn)，在其頻譜圖上的020Hz和80100Hz之間的幅值比較大，而在3070Hz之間的幅值相對(duì)較小。4.2 加入高頻干擾信號(hào)后的波形圖通過(guò)觀察加入高頻信號(hào)后的時(shí)域圖和頻域圖，將它與未加入前進(jìn)行比較，可以發(fā)現(xiàn)頻譜圖在0Hz時(shí)的幅度增加的很大，而且又在沒(méi)有譜線的頻率上出現(xiàn)了頻譜，這是由于白噪聲在所有頻率上都有頻率造成的。圖4 加入高頻干擾后的時(shí)域和頻譜圖圖5 加入電源線干擾后的時(shí)域和頻譜圖圖6 切比雪夫1型數(shù)字濾波器的