心電信號(hào)預(yù)處理(共15頁)

1、生物醫(yī)學(xué)信號(hào)處理(xn ho ch l)實(shí)習(xí)報(bào)告學(xué)生姓名：學(xué)號(hào)：實(shí)驗(yàn)室名稱：項(xiàng)目名稱：心電信號(hào)的預(yù)處理項(xiàng)目內(nèi)容：閱讀相關(guān)文獻(xiàn)資料，理解噪聲的相關(guān)知識(shí)，據(jù)此綜述抑制信號(hào)中干擾噪聲信號(hào)的基本理論和基本方法；設(shè)計(jì)兩種分別用于抑制不同噪聲的濾波器；基于W-H方程，設(shè)計(jì)最優(yōu)濾波器；運(yùn)用前面設(shè)計(jì)的幾種濾波器，對加有噪聲的模擬ECG信號(hào)進(jìn)行去噪；總結(jié)不同濾波器的去噪效果。原理（寫出具體的計(jì)算公式） 人體心電信號(hào)是一種弱電信號(hào),信噪比低一般正常的心電信號(hào)頻率范圍為0.05一100Hz,幅度為10產(chǎn)v一smV而這個(gè)心動(dòng)周期信號(hào)帶寬主要集中在058士1白Hz,P波的帶寬為0一55士19Hz,T波的帶寬為011士Z

2、Hz而90%的心電信號(hào)(ECG)頻譜能量集中在0.25/35Hz之間采集心電信號(hào)時(shí),一會(huì)受到各種噪聲的干擾,噪聲來源通常有下面幾種:(1)工頻干擾 50Hz工頻干擾是由人體的分布電容所引起,工頻干擾的模型由50Hz的正弦信號(hào)及其諧波組成由于本實(shí)驗(yàn)中所采用的是MIT一BIH數(shù)據(jù)庫的心電信號(hào),所有的工頻干擾都為60Hz幅值通常與ECG峰峰值相當(dāng)或更強(qiáng)。(2)電極接觸噪聲 電極接觸噪聲是瞬時(shí)干擾,來源于電極與肌膚的不良接觸,即病人與檢測系統(tǒng)的連接不好其連接不好可能是瞬時(shí)的,如病人的運(yùn)動(dòng)和振動(dòng)導(dǎo)致松動(dòng);也可能是檢測系統(tǒng)不斷的開關(guān)!放大器輸入端連接不好等電極接觸噪聲可抽象為快速!隨機(jī)變化的階躍信號(hào),它按

3、指數(shù)形式衰減到基線值,包含工頻成分這種瞬態(tài)過渡過程可發(fā)生一次或多次!其特征值包括初始瞬態(tài)的幅值和工頻成分的幅值!衰減的時(shí)間常數(shù);其持續(xù)時(shí)間一般為15左右,幅值可達(dá)記錄儀的最大值。(3)人為運(yùn)動(dòng) 人為運(yùn)動(dòng)是瞬時(shí)的(但非階躍)基線改變,由電極移動(dòng)中電極與皮膚阻抗改變所引起人為運(yùn)動(dòng)由病人的運(yùn)動(dòng)和振動(dòng)所引起,造成的基線干擾形狀可認(rèn)為類似周期正弦信號(hào),其峰值幅度和持續(xù)時(shí)間是變化的,幅值通常為幾十毫伏。(4)肌電干擾(EMG) 肌電干擾來自于人體的肌肉顫動(dòng),肌肉運(yùn)動(dòng)產(chǎn)生毫伏級電勢EMG基線通常在很小電壓范圍內(nèi)所以一般不明顯肌電干擾可視為瞬時(shí)發(fā)生的零均值帶限噪聲,主要能量集中在30一300Hz范圍內(nèi)。(5)

4、基線漂移和呼吸時(shí)ECG幅值的變化 基線漂移和呼吸時(shí)ECG幅值的變化一般由人體呼吸!電極移動(dòng)等低頻干擾所引起,頻率小于5Hz;其變化可視為一個(gè)加在心電信號(hào)上的與呼吸頻率同頻率的正弦分量,在0.015一0.3Hz處基線變化變化幅度的為ECG峰峰值的15%。 上面的電極接觸噪聲與人為運(yùn)動(dòng)所產(chǎn)生的噪聲是人為因素造成的，當(dāng)然也可以通過人為因素來避免。然而工頻干擾、肌電干擾(EMG)與基線漂移和呼吸時(shí)ECG幅值的變化就不是人為因素所能消除的了。 為了濾除掉上述三種噪聲，我按照實(shí)驗(yàn)要求設(shè)計(jì)了三種不同的濾波器。分別是巴特沃斯濾波器與切比雪夫?yàn)V波器。為了對比他們的濾波效果，又設(shè)計(jì)了一個(gè)維納濾波器。最后運(yùn)用SNR

5、指標(biāo)定量分析了不同濾波器的去噪能力。以下是3種濾波器的原理：巴特沃斯濾波器的設(shè)計(jì)原理其特點(diǎn)是通頻帶內(nèi)的頻率響應(yīng)曲線最大限度平坦，沒有起伏，而在阻頻帶則逐漸下降為零（對理想低通濾波的逼近：巴特沃思濾波器是以原點(diǎn)附近的最大平坦響應(yīng)來逼近理想低通濾波器）。而濾波器的幅頻特性是隨著濾波器的階次N的增加而變得越來越好，在截止頻率有：（1）衰減具有不變性。通帶、阻帶均具有單調(diào)下降的特性。低通巴特沃斯濾波器的幅頻響應(yīng)為： （2）其中，為濾波器的階數(shù)，取正整數(shù)； 為濾波器的截止頻率，為其的頻率； 為濾波器的通帶截止頻率；切比雪夫?yàn)V波器的設(shè)計(jì)原理切比雪夫?yàn)V波器有兩類，一類是在同帶內(nèi)幅度頻率響應(yīng)呈現(xiàn)等波紋特性，而

6、阻帶內(nèi)是單調(diào)的全極點(diǎn)濾波器；另一類是在通帶內(nèi)幅度頻率響應(yīng)呈現(xiàn)單調(diào)特性，而阻帶內(nèi)是等波紋特性，同時(shí)有零點(diǎn)和極點(diǎn)的濾波器，這類濾波器的零點(diǎn)位于s平面的虛軸上。第類切比雪夫?yàn)V波器的幅度平方函數(shù)定義為（3）其中，為切比雪夫多項(xiàng)式（），為濾波器階數(shù)；為通帶截止角頻率，此處是指被通帶波紋所限制的最高角頻率， ；為小于1的正數(shù)，表示通帶內(nèi)幅度波動(dòng)的程度，越小，幅度波動(dòng)越?。黄涮攸c(diǎn)為：當(dāng)時(shí)，在之間呈等波紋變化，越小，波動(dòng)幅度越??；無論為何值，所有的曲線都通過點(diǎn)，該點(diǎn)被定義為點(diǎn)；當(dāng)時(shí)，若為奇數(shù)，；若為偶數(shù)，（4）；當(dāng)時(shí)，曲線呈單調(diào)下降；通帶內(nèi)的起伏使對應(yīng)的相頻特性呈現(xiàn)非線性；維納濾波器的設(shè)計(jì)原理維納濾波器屬于現(xiàn)

7、代濾波器，傳統(tǒng)的濾波器只能濾除信號(hào)和干擾頻帶沒有重疊的情況，當(dāng)信號(hào)和干擾頻帶有重疊的時(shí)候傳統(tǒng)濾波器將無能為力，這時(shí)就需要用到現(xiàn)代濾波器，現(xiàn)代濾波器利用信號(hào)和干擾的統(tǒng)計(jì)特征（如自相關(guān)函數(shù)、功率譜等）導(dǎo)出一套最佳估計(jì)算法，然后用硬件或軟件予以實(shí)現(xiàn)。下面是對維納濾波器的設(shè)計(jì)仿真：維納濾波器是以均方誤差最?。↙MS）為準(zhǔn)則的，它根據(jù)過去觀測值和當(dāng)前觀測值來估計(jì)信號(hào)的當(dāng)前值，因此它的解形式是系統(tǒng)的傳遞函數(shù)或單位脈沖響應(yīng) 圖一設(shè)計(jì)維納濾波器的過程就是尋求在最小均方誤差下濾波器的單位脈沖響應(yīng)或傳遞函數(shù)的表達(dá)式，其實(shí)質(zhì)就是解維納-霍夫（Wiener-Hopf）方程。均方誤差為：(5)維納霍夫方程最小均方誤差是

8、下的解為：(6)其中，為濾波器的系數(shù)向量，為含有噪聲的混合信號(hào)的自相關(guān)矩陣，為混合信號(hào)和原始信號(hào)的互相關(guān)向量，為此先求出混合信號(hào)的自相關(guān)函數(shù)以及混合信號(hào)和原始信號(hào)的互相關(guān)函數(shù)，這兩個(gè)函數(shù)，我們需要有樣本得到。編寫的源程序：巴特沃斯濾波器(Butter.m) B,A=butter(3,2*pi*1000,s); %脈沖響應(yīng)不變法低通濾波器 num1,den1=impinvar(B,A,4000); h1,w=freqz(num1,den1); B,A=butter(3,2/0.00025,s); %雙線性變換法低通濾波器 num2,den2=bilinear(B,A,4000); h2,w=fr

9、eqz(num2,den2); f=w/pi*2000; plot(f,abs(h1),BLUE,f,abs(h2),RED); legend(脈沖響應(yīng)不變法,雙線性變換法); grid; xlabel(頻率/Hz ) ylabel(幅值/dB) title(巴特沃斯低通濾波器);圖二2.切比雪夫?yàn)V波器（qiebixuefu1.m） clc,clear %確定數(shù)字濾波器指標(biāo) Fs=10000; %采樣頻率 fp=280;fs=450; wp=2*pi*fp/Fs; %通帶截止頻率 Ap=0.1; %通帶內(nèi)的最大衰減 ws=2*pi*fs/Fs; %阻帶截止頻率 As=40; %阻帶內(nèi)的最小衰減

10、 %數(shù)字角頻率轉(zhuǎn)換成模擬角頻率 %雙線性變換法 wps=2*Fs*tan(wp/2); wss=2*Fs*tan(ws/2); %設(shè)計(jì)模擬濾波器 k1=sqrt(10(0.1*As)-1)/(10(0.1*Ap)-1); qsp=wss/wps; N=ceil(acosh(k1)/acosh(qsp); e=sqrt(10(0.1*Ap)-1); Ap=10*log(1+e2); %計(jì)算得到模擬低通H(s)分子和分母的系數(shù) z,p,k=cheb1ap(N,Ap); bp,ap=zp2tf(z,p,k); bs,as=lp2lp(bp,ap,Fs); %復(fù)變量映射s-z b,a=bilinear

11、(bs,as,Fs); %畫出切比雪夫1型低通濾波器的幅頻響應(yīng)曲線圖 H,w =freqs(bp,ap,2000);%按n指定的頻率點(diǎn) 給出頻率響應(yīng) magH2 = (abs(H).2; %給出傳遞函數(shù)的幅度平方 %輸出波形 plot(w,magH2);title(切比雪夫1型低通濾波器的幅頻響應(yīng)特性 ); xlabel (w/wc); ylabel(Chebyshev 1 | H(jw)|2);圖三3.維納濾波器（Winner.m） functionH,Emin=Winner(Rss,Rww,N) L1=(length(Rs)+1)/2; Rss=zeros(1,L1); for k=1:L

12、1 Rss(k)=Rs(k+L1-1); end L2=(length(Rw)+1)/2; Rww=zeros(1,L2); for k=1:L2 Rww(k)=Rw(k+L2-1); end Rx=zeros(1,N); for k=1:N Rx(k)=Rss(k)+Rww(k); end Rxx=zeros(N,N); for i=1:N for j=1:N if(i=j) Rxx(i,j)=Rx(j-i+1); else Rxx(i,j)=Rx(i-j+1); end end end H=pinv(Rxx)*Rss; Emin=Rss(1)-sum(H*Rss);測試程序：（Winner

13、Test.m） clc; t=-50:0.1:50; w=0.1*pi; s=cos(w.*t); w=randn(1,length(t); x=s+w; Rs=xcorr(s);% 信號(hào)的自相關(guān)函數(shù) Rw=xcorr(w); % 白噪聲的自相關(guān)函數(shù) h,E=Winner(Rs,Rw,length(t);% 維納濾波 se=filter(h,1,x); plot(t,s); hold on plot(t,x,g); hold on plot(t,se,r); e=s-se; plot(t,e,BLACK); legend(原始信號(hào),加噪后信號(hào),維納濾波后信號(hào),誤差);圖四4.分別運(yùn)用上面設(shè)計(jì)的

14、濾波器和最優(yōu)濾波器對含高斯噪聲心律失常ECG數(shù)據(jù)消噪，并運(yùn)用SNR指標(biāo)定量分析濾波器的去噪能力1）設(shè)計(jì)SNR（信噪比）的計(jì)算方法 function snr=SNR(I,In) %計(jì)算信噪比函數(shù) %I:original signal %In:noisy signal(ie.original signal+noise signal) snr=0; Ps=sum(abs(I).2);%signal power Pn=sum(abs(In-I).2);%noise power snr=10*log10(Ps/Pn); %其中I是純信號(hào)，In是帶噪聲的信號(hào)，snr是信噪比巴特沃茲低通濾波器(snrbut

15、ter.m)clc;clear;Test_Signal_1=load(118e00m.mat); % 導(dǎo)入數(shù)據(jù)Signal_1=Test_Signal_1.val(1,:);Fs=2*max(abs(Signal_1); %采樣頻率確定Signal_1=Signal_1/(Fs/2); %數(shù)據(jù)歸一化 for i=1:2000 Signal1(i)=Signal_1(1,i);endsubplot(4,1,1),plot(Signal1); title(原始信號(hào)) Signal_1_1 = awgn(Signal_1,8); %高斯白噪聲for i=1:2000 Signal1_1(i)=Sig

16、nal_1_1(1,i);endsnr_1_1=SNR(Signal1,Signal1_1)subplot(4,1,2),plot(Signal1_1);title(加入高斯噪聲以后的信號(hào)) fp = 100; %通帶截止頻率 fs = 105; %阻帶截止頻率Rp=1; Rs=40; %通帶波紋 阻帶衰減 Fs=360;%計(jì)算歸一化角頻率wp=fp/(Fs/2); ws=fs/(Fs/2); %數(shù)字轉(zhuǎn)模擬wp = tan(pi*wp/2);ws = tan(pi*ws/2);N,wc=buttord(wp,ws,Rp,Rs,s);B,A=butter(N,wc,s); b,a=bilinea

17、r(B,A,.5); %雙線性變換法h, w = freqz(b, a, 256, Fs);subplot(4,1,3)plot(w, abs(h);title(巴特沃茲濾波器) Signal_N=filter(b,a,Signal1_1);snr_N=SNR(Signal1,Signal_N)subplot(4,1,4),plot(Signal_N);title(濾波以后圖像)切比雪夫低通濾波器(snrqiebixuefu1.m)clc;clear;Test_Signal_1=load(118e00m.mat); %導(dǎo)入數(shù)據(jù)Signal_1=Test_Signal_1.val(1,:);Fs

18、=2*max(abs(Signal_1); %采樣頻率確定Signal_1=Signal_1/(Fs/2); %歸一化 for i=1:2000 Signal1(i)=Signal_1(1,i);endsubplot(4,1,1),plot(Signal1); title(原始信號(hào))Signal_1_1 = awgn(Signal_1,8); %高斯白噪聲for i=1:2000 Signal1_1(i)=Signal_1_1(1,i);endsnr_1_1=SNR(Signal1,Signal1_1)subplot(4,1,2),plot(Signal1_1);title(加入高斯噪聲以后的

19、信號(hào))% Filter Designfp = 100; %通帶截止頻率 fs = 105; %阻帶截止頻率Fs=360;Rp=1; Rs=40; %通帶波紋 阻帶衰減 %計(jì)算歸一化角頻率wp=fp/(Fs/2); ws=fs/(Fs/2); %數(shù)字轉(zhuǎn)模擬wp = tan(pi*wp/2);ws = tan(pi*ws/2);N,wc=cheb1ord(wp,ws,Rp,Rs,s); %切比雪夫一型B,A=cheby1(N,Rp,wc,low,s); b,a=bilinear(B,A,.5); %雙線性變換法h, w = freqz(b, a, Fs);subplot(4,1,3)plot(w,

20、 abs(h);title(切比雪夫?yàn)V波器)Signal_N=filter(b,a,Signal1_1);snr_N=SNR(Signal1,Signal_N)subplot(4,1,4),plot(Signal_N);title(去噪之后的信號(hào))維納濾波器(snrWinner.m)clear,clc;n = input(信號(hào)的長度 n = );Test_Signal_1=load(118e00m.mat); % Read data signals 1-2 minutesSignal_1=Test_Signal_1.val(1,:);Fs=2*max(abs(Signal_1); %采樣頻率確

21、定Signal_1=Signal_1/(Fs/2); %歸一化for i=1:2000 % Initialization data Signal1(i)=Signal_1(1,i);endsubplot(3,1,1),plot(Signal1); title(原始信號(hào))Signal_1_1 = awgn(Signal_1,8); % Gaussian noisefor i=1:2000 Signal1_1(i)=Signal_1_1(1,i);endsnr_1_1=SNR(Signal1,Signal1_1)subplot(3,1,2),plot(Signal1_1);title(加入高斯噪聲以后的信號(hào))%維納濾波Mlag=10