濾波器設計及在心電信號濾波中的應用

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

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

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

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

5、與信息處理相關問題的能力，為以后的工作和學習打下基礎。二 設計原理 2.1 Butterworth低通數(shù)字濾波器的設計巴特沃斯低通濾波器的平方幅度響應為其中，n為濾波器的階數(shù)，為低通濾波器的截止頻率。 該濾波器具有 一些特殊的性質： 對所有的n,都有當時, ； 對所有的n,都有當時, ；是的單調遞減函數(shù),即不會出現(xiàn)幅度響應的起伏； 當時,巴特沃斯濾波器趨向于理想的低通濾波器； 在處平方幅度響應的各級導數(shù)均存在且等于0，因此在該點上取得最大值，且具有最大平坦特性。 圖1展示了2階、4階、8階巴特沃斯低通濾波器的幅頻特性?？梢婋A數(shù)n越高，其幅頻特性越好，低頻檢測信號保真度越高，過渡帶變窄，即衰減加

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

7、0.4rad，ap=1dB； s=0.5rad，as=40dB（6）模擬低通的技術指標為： 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)利用模擬切比雪夫濾波器設計數(shù)字濾波器。通帶截止頻率為：wp=0.4*pi; 阻帶截止頻率為：ws=0.5*pi;通帶最大衰減為：Rp=1;阻帶最大衰減為：As=15;設定周期為1s；模擬低通濾波器的生成：b,a=cheby1(n,1,Wn,'low','s');滿足設計

8、指標的最小階數(shù)和截止頻率：Wnn,Wn=cheb1ord(OmegaP,OmegaS,1,40,'s')。最后實現(xiàn)輸入輸出、幅頻特性、相頻特性的圖形。三 詳細設計步驟3.1心電數(shù)據(jù)的導入將老師給的心電信號原始數(shù)據(jù)存于桌面文件夾，導入數(shù)據(jù)代碼為：a=load('C:UsersAdministratorDesktop題目七數(shù)據(jù)文件題目七數(shù)據(jù)文件ecg.txt');3.2 繪出心電信號的時域圖和頻譜圖將導入的數(shù)據(jù)分別用t來替換，通過調用plot函數(shù)來畫出時域圖，然后通過對4000個心電數(shù)據(jù)的幅值進行FFT運算，再次調用plot函數(shù)來繪出頻域圖，具體設計如下：figur

9、e(1); %新建圖像subplot(2,1,1); %將2個圖畫到一個平面t=a(1:4000,1); %用t替換導入的4000個數(shù)據(jù)plot(t); %繪出t的圖形title('原始波形圖');xlabel('時間（s）');ylabel('幅值（A）');y1=fft(a(:,1),4000); %行取全部f1=100*(0:3999)/4000; %先生成一個0,1,2,.,3999的整數(shù)向量，然后對對這個向量的每一項乘以100除以4000.subplot(2,1,2); %將2個圖畫到一個平面plot(f1,abs(y1); %f1為橫

10、坐標，abs(y1)為縱坐標作圖title('原始頻譜圖');xlabel('頻率（Hz）');ylabel('幅度(dB)');3.3 加入干擾這里我加入的干擾是：白噪聲高頻干擾，50Hz的電源線干擾。3.3.1 白噪聲通過用s來代表加入白噪聲后的信號，并進行數(shù)字濾波器的頻率響應，對s中4000個頻率點調用plot函數(shù)畫出加入白噪聲后的時域圖，再對ws/pi,abshs調用plot函數(shù)畫出加入白噪聲后的頻譜圖，具體操作如下：q=50*rand(4000,1); %產(chǎn)生4000行1列的位于（0，1）區(qū)間的隨機數(shù)s=a(:,1)+q; %產(chǎn)生的隨機

11、數(shù)與原始信號疊加賦給shs,ws=freqz(s,1,4096); %求離散系統(tǒng)頻響特性的函數(shù)freqz()abshs=abs(hs); %hs取絕對值賦給abshsfigure(2); %新建圖像2subplot(2,1,1); plot(s(1:4000); %繪出從1取到4000時s的圖形title('加入高頻干擾信號后的時域圖');xlabel('時間（s）');ylabel('幅值（A）');subplot(2,1,2);plot(ws/pi,abshs); % ws/pi為橫坐標abshs為縱坐標作圖title('加入高頻干擾

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

13、25*(4000-1); %從0開始間隔0.00025取值到 0.00025*(4000-1)x1=a(:,1); %將矩陣a用x1表示x=x1+x2; %x1和x2疊加賦給xy2=fft(x,4000); %對x信號做快速傅里葉變換f2=100*(0:3999)/4000; %先生成一個0,1,2,.,3999的整數(shù)向量，然后對對這個向量的每一項乘以100除以4000.figure(3); %新建圖像3subplot(2,1,1);plot(t,x); %t為橫坐標，x為縱坐標作圖title('加入電源線干擾后的時域圖');xlabel('時間（s）');yl

14、abel('幅值（A）');subplot(2,1,2);plot(f2,abs(y2); %f2為橫坐標，abs(y2)為縱坐標作圖title('加入電源線干擾后的頻譜圖');xlabel('幅');ylabel('Hz');3.4 濾波器的設計3.4.1切比雪夫I型數(shù)字低通濾波器用Wp1,Wp2,Ws1,Ws2表示分別用通帶和阻帶截止頻率的角頻率算出頻帶寬帶，計算階數(shù)n1和截止頻率WN，再設計切比雪夫I型模擬濾波器，采用雙線性法將模擬濾波器系數(shù)變?yōu)閿?shù)字濾波器系數(shù)，畫出切比雪夫I型數(shù)字濾波器的頻率響應，調用filter實現(xiàn)對白噪

15、聲的濾波，再最后調用plot函數(shù)畫出濾除白噪聲后的時域圖和頻域圖。具體過程如下：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ù)字轉化為模擬Omegap2=2*fs*tan(Wp2/2); %數(shù)字轉化為模擬BW=Omegap2-Omegap1; %頻帶寬度W0=Omegap1*Omegap2; W00=sqrt(W0);WP=1;n1,WN=buttord(WP,WS,1,50

16、,'s'); %參數(shù)WP和WS分別是通帶邊界頻率和阻帶邊界頻率，通帶最大衰減為1阻帶最小衰減50，返回的參數(shù)n1和WN分別為濾波器的階數(shù)和3dB截止頻率B,A=cheby1(n1,1,WN,'s'); %計算出階數(shù)為n1、截止頻率為WN、通帶波紋最大衰減為1的數(shù)字低通濾波器，它的返回值a、b分別表示數(shù)字低通濾波器的系統(tǒng)函數(shù)的分子和分母的多項式系數(shù)BT,AT=lp2bp(B,A,W00,BW); num,den=bilinear(BT,AT,0.5); %實現(xiàn)模數(shù)的映射freqz(num,den,64); y=filter(num,den,s);figure(5)

17、;subplot(2,1,1);plot(y);title('濾除高頻信號后的時域圖');xlabel('時間（s）');ylabel('幅值（A）');s3=fft(y,4000);f1=100*(0:3999)/4000; subplot(2,1,2);plot(f1,abs(s3);title('濾除高頻信號后的頻譜圖');xlabel('頻率(Hz)');ylabel('幅值(dB)');3.4.2 高通濾波器的設計 高通濾波器的設計可通過模擬低通濾波器再經(jīng)頻率變換而實現(xiàn)。其中模擬低通濾波器

18、可根據(jù)已經(jīng)存在的典型濾波器，如巴特沃斯濾波器等逼近實現(xiàn)，而由低通到高通轉換理論依據(jù)在緒論部分已經(jīng)進行了詳細的論證，又本設計基于MATLAB仿真軟件實現(xiàn)，可利用MATLAB信號處理工具箱提供的各種函數(shù)模型實現(xiàn)，可忽略其中的一些復雜的函數(shù)變換，從而簡化理論設計和論證，具體如下：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); %調用Butterworth模擬低通濾波器原型函數(shù)b,a=zp2tf(

19、z,p,k); %將零點極點增益形式轉換為傳遞函數(shù)形式H,w=freqs(b,a,n); %按n指定的頻率點給出頻率響應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') %對不同曲線做標記text(1.3,0.08,'n=5')text(1.16,0.08,'n=10

20、9;)text(0.93,0.98,'n=20')grid on; %模擬高通濾波器設計m=0:0.01:2;for ii=1:4 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); %由低通原型濾波器轉換為截止頻率為0.25Hz的高通濾波器Ht,w=freqs(bt,at,m); subplot(2,1,2);hold on; %hold on 是保存axes內(nèi)圖像用的,如果在新

21、畫圖像之后不想覆蓋原圖像就要加上hold on這句話plot(w,abs(Ht);endtitle('模擬高通濾波器');xlabel('w/pi');ylabel('|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

22、,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;title(' 模擬高通濾波器幅頻特性');xlabel('w/wc');ylabel('幅度/dB');3.4.3 帶阻濾波器的設計 設計帶阻濾波器，再畫出濾波器的頻率響應圖，然后用filter實現(xiàn)對工頻

23、干擾信號的濾波，調用plot函數(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('濾除電源線干擾后的時域圖');xlabel('時間（s）');ylabel('幅值（A）');subplot(2,1,2);y5=fft(y4,4000);f1=100*(0:3999)/4000;plot(f1,abs(

24、y5);title('濾除電源線干擾后的頻譜圖');xlabel('頻率(Hz)');ylabel('幅值(dB)');四 設計結果及分析4.1 心電數(shù)據(jù)的導入，心電信號的時域圖和頻譜圖如下：圖3 原始時域和頻譜圖通過導入的心電信號數(shù)據(jù)發(fā)現(xiàn)，在其頻譜圖上的020Hz和80100Hz之間的幅值比較大，而在3070Hz之間的幅值相對較小。4.2 加入高頻干擾信號后的波形圖通過觀察加入高頻信號后的時域圖和頻域圖，將它與未加入前進行比較，可以發(fā)現(xiàn)頻譜圖在0Hz時的幅度增加的很大，而且又在沒有譜線的頻率上出現(xiàn)了頻譜，這是由于白噪聲在所有頻率上都有頻率造成的。圖4 加入高頻干擾后的時域和