基于MATLAB的心電信號的分析與處理設(shè)計(jì)

河南科技大學(xué)程設(shè)計(jì)說明書課程名稱 醫(yī)學(xué)信號處理題 目基于MATLAB的心電信號的分析與處理設(shè)計(jì) （2）院 系 醫(yī)學(xué)技術(shù)與工程學(xué)院班 級 醫(yī)療器械工程 111班學(xué)生姓名指導(dǎo)教師 侯海燕宋衛(wèi)東＿日期2014年9月11號課程設(shè)計(jì)任務(wù)書（指導(dǎo)教師填寫）課程設(shè)計(jì)名稱醫(yī)學(xué)信號處理 學(xué)生姓名 專業(yè)班級醫(yī)療器械工程 111班設(shè)計(jì)題目 基于 MATLA的心電信號的分析與處理設(shè)計(jì)（B 2）一、課程設(shè)計(jì)目的1．熟練掌握使用 MATLA程序設(shè)計(jì)方法B2．掌握數(shù)字信號處理的基本概念、理論、方法3．掌握序列離散傅里葉變換的 MATLA實(shí)現(xiàn)，并進(jìn)行頻譜分析B4．熟練掌握使用 MATLA設(shè)計(jì)BIIR或FIR數(shù)字濾波器5．學(xué)會用 MATLA對信號進(jìn)行分析和處理B二、設(shè)計(jì)內(nèi)容、技術(shù)條件和要求一）設(shè)計(jì)內(nèi)容與技術(shù)條件1．根據(jù)給定的一段 MIT-BIH心電信號（ 101號），畫出心電信號的時域波形和頻譜圖（幅頻和相頻）；2．根據(jù)心電信號頻率范圍及其噪聲的頻率范圍設(shè)計(jì) 2個濾波器（一個 IIR，一個FIR）實(shí)現(xiàn)對心電信號濾波。濾波器的種類（高通，低通，帶通，帶阻），濾波器性能指標(biāo)（通阻帶截止頻率，衰減系數(shù)），濾波器的設(shè)計(jì)方法（ IIR有沖擊響應(yīng)不變法和雙線性變換法， FIR有窗函數(shù)法及頻率抽樣法）等自行設(shè)計(jì)。要求輸出所設(shè)計(jì)的濾波器的系統(tǒng)函數(shù)，畫出濾波器的頻率響應(yīng)（幅頻響應(yīng)和相頻響應(yīng)）曲線；3．用該濾波器對心電信號進(jìn)行濾波，畫出濾波以后心電信號的時域波形和頻譜（幅頻）；分析信號濾波前后心電信號的時域和頻域的變化；4．兩個濾波器濾波效果異同分析；5．運(yùn)用 GUI設(shè)計(jì)一個心電信號處理系統(tǒng)界面。（選作）二）設(shè)計(jì)要求1．根據(jù)濾波器的性能指標(biāo)要求，設(shè)計(jì)數(shù)字濾波器；2．程序中按照 IIR濾波器的步驟一步步完成設(shè)計(jì)； 盡可能的少調(diào)用 MATLA自帶B的函數(shù)文件；

3．設(shè)計(jì)程序要通用性好，整齊易懂，并要求主要語句有注釋；4．設(shè)計(jì)結(jié)果中的 圖示要美觀，整齊，有標(biāo)題，有縱橫坐標(biāo)標(biāo)示；5．課程設(shè)計(jì)報(bào)告要 有理論依據(jù)、設(shè)計(jì)過程，結(jié)果分析 。報(bào)告要求實(shí)事求是、文理通順、字跡端正。三、時間進(jìn)度安排9月1日：查閱資料， 復(fù)習(xí)醫(yī)學(xué)信號處理課程的相關(guān)內(nèi)容， 熟悉 MATLA的工作環(huán)境，B掌握運(yùn)算編程方法。9月2日~5日：按照課程設(shè)計(jì)要求，完成課程設(shè)計(jì)的程序設(shè)計(jì)部分。9月8日~10日：調(diào)試并完善程序，撰寫課程設(shè)計(jì)說明書。9月11~12日：答辯。四、主要參考文獻(xiàn)參考資料（1）數(shù)字信號處理——理論、算法與實(shí)現(xiàn)胡廣書清華大學(xué)出版社（2）數(shù)字信號處理及 MATLA實(shí)現(xiàn)余成波清華大學(xué)出版社B（3）MATLAB7.在數(shù)字信號處理中的應(yīng)用羅軍輝機(jī)械工業(yè)出版社0（4）數(shù)字信號處理原理及其 MATLA實(shí)現(xiàn)從玉良電子工業(yè)出版社B（5）MATLA幫助文件B相關(guān)網(wǎng)站MATLA中國論壇B MATLA學(xué)習(xí)網(wǎng)B MATLA中文論壇B 指導(dǎo)教師簽字：侯海燕2014年指導(dǎo)教師簽字：侯海燕2014年8月25日1、設(shè)計(jì)目的意義 ····································12、設(shè)計(jì)內(nèi)容 ········································1設(shè)計(jì)原理 ········································1要求 ············································1內(nèi)容 ············································13、計(jì)過程及分析 ····································1原始心電信號分析 ································1設(shè)計(jì)濾波器 ·····································34、結(jié)果分析 ········································75、參考文獻(xiàn) ·······································7、設(shè)計(jì)目的意義本次課程主要達(dá)到以下幾個目的：熟練掌握使用 MATLAB程序設(shè)計(jì)方法掌握序列離散傅里葉變換的 MATLAB實(shí)現(xiàn)，并進(jìn)行頻譜分析熟練掌握使用 MATLAB設(shè)計(jì) IIR或FIR數(shù)字濾波器2、設(shè)計(jì)內(nèi)容設(shè)計(jì)原理頻譜分析（傅里葉變換）（2）濾波器設(shè)計(jì)IIR（沖激響應(yīng)不變法，雙線性變換法） 詳解見附錄FIR（窗函數(shù)法，頻率抽樣法） 詳解見附錄（3）信號濾波（時域卷積，頻域相乘）要求要求設(shè)計(jì)出心電數(shù)據(jù)處理的處理與分析程序。處理對象：心電數(shù)據(jù)；內(nèi)容：心電數(shù)據(jù)仿真，心電數(shù)據(jù)處理（仿真數(shù)據(jù)，真實(shí)數(shù)據(jù)）；結(jié)果：得到處理結(jié)果。內(nèi)容始心電信號的時域波形和頻譜圖（幅頻和相頻）。IIR濾波器的頻率響應(yīng)（幅頻響應(yīng)和相頻響應(yīng)）曲線；原始心電信號經(jīng)IIR濾波器濾波后的心電信號的時域波形和頻譜圖（幅頻和相頻）；FIR濾波器的頻率響應(yīng)（幅頻響應(yīng)和相頻響應(yīng)）曲線；原始心電信號經(jīng)FIR濾波器濾波后的心電信號的時域波形和頻譜圖（幅頻和相頻）；3、設(shè)計(jì)過程及分析原始心電信號分析用load函數(shù)將原心電信號導(dǎo)入 課程設(shè)計(jì) \101ecg.txt'），并畫出心電信號的時域波形和頻譜圖（幅頻和相頻） ，如圖1所示 :原始心電信號的時域波形圖2原始心電信號的時域波形圖2-2TOC\o"1-5"\h\z0 1 2 3 4 5 6 7 8 9 1015010050050原始心電信號的頻譜始心電信號的頻譜圖0100 150 200 250 300 3500頻率/HZ原始心電信號的相頻圖55d/ra0角相-50 50 100 150 200 250 300 350頻率/Hz圖1原始心電信號的時域波形圖及頻譜圖心電信號由于受到人體諸多因素的影響，因而有著一般信號所沒有的特點(diǎn)：信號弱。心電信號是體表的電生理信號，一般比較微弱，幅度在 10pV～5mV，頻率為 0．05～100Hz。例如從母體腹部收取到的胎兒心電信號僅 10／zV～50／IV。噪聲強(qiáng)。由于人體自身信號弱，加之人體又是一個復(fù)雜的系統(tǒng)，因此信隨機(jī)性強(qiáng)。 心電信號不僅是隨機(jī)的，而且是非平穩(wěn)的。同時，在心電圖從而使圖像質(zhì)量變差， 使均勻和連續(xù)變在心電圖上形成一些毛刺。 使原本很微弱的信號很難和工頻干擾50Hz(MIT-BIH數(shù)60Hz)的基波及其各次諧波，其幅值成分ECG峰一峰值的 0—50％范圍內(nèi)變化。引起基線漂移的干擾心電信號有時候會出現(xiàn)信號基線起伏不平的現(xiàn)象， 造成這樣的現(xiàn)象有很多原因，主要的有：呼吸運(yùn)動人體呼吸時胸腔內(nèi)器官和組織會發(fā)生一定程度的變化，會對在體表記錄到的心電圖波形的幅度和形態(tài)有所影響， 表現(xiàn)為基線隨呼吸產(chǎn)生周期性或非周期性漂移， 從而導(dǎo)致心電波形的幅度隨呼氣和吸氣而分別上抬和下移。 呼吸運(yùn)動是引起心電基線漂移的主要原因。運(yùn)動偽跡運(yùn)動偽跡是由于人體輕微運(yùn)動造成電極與入體的接觸電阻發(fā)生變化而引入的一種干擾， 它的產(chǎn)生原因僅僅是接觸電阻的變化， 而不是接觸的斷續(xù)。這種干擾同樣導(dǎo)致信號基線的變化， 但不是基線的躍變。 ③信號記錄和處理中電子設(shè)備引起的干擾這種干擾對信號影響很大， 嚴(yán)重時可完全淹沒心電信號或使得基線劇烈漂移， 其中導(dǎo)聯(lián)開路和放大器的熱移是主要因素。 這種干擾往往無法通過心電分析算法來校正。 由于心電波形已經(jīng)完全畸變， 此時對這些數(shù)據(jù)分析已無太大意義。所以一般跳過此段數(shù)據(jù)。高頻噪聲心電信號中的高頻噪聲主要是肌電噪聲。 肌肉收縮會產(chǎn)生 mV級的肌電干擾，表現(xiàn)為心電圖上不規(guī)則的細(xì)小波紋， 使心電圖模糊不清或產(chǎn)生失真。 肌電噪聲的特點(diǎn)是頻率范圍較廣，頻譜分布非常復(fù)雜。3.2設(shè)計(jì)濾波器IIR濾波器的設(shè)計(jì))IIR濾波器的設(shè)計(jì)過程：按照技術(shù)要求設(shè)計(jì)一個模擬濾波器， 得到模擬低通濾波器的傳輸函數(shù) H(s)，再按一定的轉(zhuǎn)換關(guān)系將 H(s)轉(zhuǎn)換成數(shù)字低通濾波器的系數(shù)函數(shù) H(z)。這樣設(shè)計(jì)的關(guān)鍵問題就是找到這樣的轉(zhuǎn)換關(guān)系，將 s平面上的 H(s)轉(zhuǎn)換成z平面上的H(z)。2)巴特沃斯濾波器分母多項(xiàng)式的因式表示，如表 1所示：表1巴特沃斯濾波器分母多項(xiàng)式的因式表示

33)巴特沃斯低通濾波器的階數(shù)公式N=log10((10^(As/10)-1)/(10^(Rp/10)-1))/(2*log10(ws/wp)))4)巴特沃斯低通濾波器函數(shù)由巴特沃斯低通濾波器的階數(shù)公式和巴特沃斯濾波器分母多項(xiàng)式的因式表示求出歸一化巴特沃斯低通濾波器 Has(s)N=7則Has(s)=1/((s+1)*(s^2+0.4450s+1)*(s^2+1.247s+1)*(s^2+1.8022s+1))5)巴特沃斯低通濾波器的頻域特性， ，如圖 2所示：低通濾波器幅頻特性1.51度幅0.50 20 40 60 80 100 120 140 160 180頻率低通濾波器幅頻特性1.51度幅0.50 20 40 60 80 100 120 140 160 180頻率/Hz圖2巴特沃斯低通濾波器的相頻和幅頻特性06)經(jīng)過巴特沃斯低通濾波器器后心電信號的時域波形和頻譜圖，如圖 3所示：2濾波后時域波形-20 1 2 3 4 5 6 7 8 9 10時間t/s圖3經(jīng)過巴特沃斯低通濾波器器后心電信號的時域波形圖和頻譜圖對比原始信號的時域波形圖和頻譜圖可得通過低通濾波器后的心電信號波形圖可以明顯看出波形變得平滑，由工頻干擾產(chǎn)生的毛刺被低通濾波器成功濾除。FIR濾波器的設(shè)計(jì))FIR濾波器的設(shè)計(jì)過程：給定理想的頻率響應(yīng)函數(shù) Hd(e^jw)及技術(shù)指標(biāo) δ,Δw;求出理想的單位抽樣響應(yīng) hd(n);根據(jù)阻帶衰減選擇窗函數(shù) w(n);根據(jù)過渡帶寬度確定 N值N=A/Δw;求所設(shè)計(jì)的 FIR濾波器的單位脈沖響應(yīng) h(n)=hd(n)*w(n);計(jì)算頻率響應(yīng) Hd(e^jw),驗(yàn)算指標(biāo)是否滿足要求。2)布拉克曼窗低通濾波器的頻域特性，如圖 4所示：100)Bd(edutingaM0000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1NormalizedFrequency( 100)Bd(edutingaM0000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1NormalizedFrequency( rad/sample)4布拉克曼窗低通濾波器的相頻和幅頻特性3）經(jīng)過布拉克曼窗低通濾波器器后心電信號的時域波形和頻譜圖，如圖 5所示/A值幅/A值幅時間t/s濾波后心電信號的頻譜圖b150/d100值50幅00 50 100 150 200 250 300 350頻率/Hz濾波后心電信號的相頻特性對比原始信號的時域波形圖和頻譜圖可得通過低通濾波器后的心電信號波頻率/Hz5對比原始信號的時域波形圖和頻譜圖可得通過低通濾波器后的心電信號波頻率/Hz5經(jīng)過布拉克曼窗低通濾波器器后心電信號的時域波形和頻譜圖4、結(jié)果分析通過用matlab設(shè)計(jì)的IIR低通濾波器和 FIR低通濾波器對帶有噪聲的心電信號進(jìn)行處理和分析，消除了原心電信號中由工頻干擾產(chǎn)生的毛刺現(xiàn)象，雖然仍存在一些誤差，但基本上達(dá)到了預(yù)期效果。通過做課程設(shè)計(jì)熟悉并掌握了相關(guān)的matlab操作環(huán)境；進(jìn)一步加深了對濾波器設(shè)計(jì)相關(guān)知識的理解，熟練了濾波器設(shè)計(jì)的方法和過程。5、參考文獻(xiàn)：TOC\o"1-5"\h\z陳天華.數(shù)字圖像處理 [M].北京.清華大學(xué)出版社， 2009劉衛(wèi)國.MALTAB程序設(shè)計(jì)與應(yīng)用 [M].北京.高等教育出版社， 2008程正興.小波分析算法與應(yīng)用 .西安.西安交通大學(xué)出版社， 1998程佩青.數(shù)字信號處理 .北京.清華大學(xué)出版社， 2007蘇金明，王永利 .MALTAB應(yīng)用指南 [M].上冊.北京電子工業(yè)出版社， 2004夏良正.數(shù)字圖像處理（修訂版） [M].南京.東南大學(xué)出版社， 1999霍紅濤，林小竹，何薇 .數(shù)字圖像處理 [M].北京.北京理工大學(xué)出版社， 2003張開滋，郭繼鴻，劉海洋 .臨床心電信息學(xué) [M].長沙.湖南科技出版社， 2002黃寶晨，朱怡然 .心電圖基本知識 [J].中國鄉(xiāng)村醫(yī)藥雜志， 2004（第7頁）許原.心電圖解讀心電圖如何解讀和診斷 [J].中國臨床醫(yī)生， 2004（第5頁）附錄：巴特沃斯數(shù)字低通濾波器：巴特沃斯濾波器是電子濾波器的一種。 巴特沃斯濾波器的特點(diǎn)是通頻帶的頻率響應(yīng)曲線最平滑。 巴特沃斯濾波器的特點(diǎn)是通頻帶內(nèi)的頻率響應(yīng)曲線最大限度平坦，沒有起伏，而在阻頻帶則逐漸下降為零。 在振幅的對數(shù)對角頻率的波得圖上,從某一邊界角頻率開始 ,振幅隨著角頻率的增加而逐步減少 ,趨向負(fù)無窮大。巴特沃斯濾波器的振幅對角頻率單調(diào)下降， 并且也是唯一的無論階數(shù)， 振幅對角頻率曲線都保持同樣的形狀的濾波器。 只不過濾波器階數(shù)越高， 在阻頻帶振幅衰減速度越快。 其他濾波器高階的振幅對角頻率圖和低級數(shù)的振幅對角頻率有不同的形狀。窗函數(shù)法：窗函數(shù)法是設(shè)計(jì) FIR濾波器的最主要方法之一， 實(shí)際中遇到的離散時間信號總是有限長的， 因此不可避免的要遇到數(shù)據(jù)截短的問題， 在信號處理中， 對離散序列的截短是通過序列與窗函數(shù)相乘來實(shí)現(xiàn)的。 在信號處理中， 窗函數(shù)是一種除在給定區(qū)間之外取值均為 0的實(shí)函數(shù)。 譬如：在給定區(qū)間內(nèi)為常數(shù)而在區(qū)間外為0的窗函數(shù)被形象地稱為矩形窗。 任何函數(shù)與窗函數(shù)之積仍為窗函數(shù)， 所以相乘的結(jié)果就像透過窗口“看”其他函數(shù)一樣。窗函數(shù)在光譜分析、濾波器設(shè)計(jì)以及音頻數(shù)據(jù)壓縮等方面有廣泛的應(yīng)用。程序流程圖：clearall;closeall;%時域波形圖以及頻譜圖課程設(shè)計(jì) \101ecg.txt');t=b(:,1);%時間c=b(:,2);%幅值figure(1);subplot(311);plot(t,c);title('原始心電信號的時域波形圖 ');xlabel('時間t/s');ylabel('幅值/A');n=3600;m=(fft(c,n));fs=360;%采樣頻率f=fs/n*(0:n-1);%定位坐標(biāo)subplot(312);plot(f,abs(m));title('原始心電信號的頻譜圖 ');xlabel('頻率/HZ');ylabel('幅值/db');axis([0,360,0,150]);subplot(313);plot(f,angle(m));title('原始心電信號的相頻圖 ');xlabel('頻率/Hz');ylabel('相角/rad');axis([0,360,-5,5]);%低通IIR濾波器;wp=2*pi*10/fs; %通帶數(shù)字頻率轉(zhuǎn)換成模擬頻率ws=2*pi*15/fs; %通帶截至數(shù)字頻率轉(zhuǎn)換成模擬頻率rp=2; %通帶最大衰減rs=20; %阻帶最小衰減[N,wc]=buttord(wp,ws,rp,rs,'s');%確定最小階數(shù) N和頻率參數(shù) Wc[Bz,Az]=butter(N,wc);%得巴特沃斯 歸一化低通原型[H,w]=freqz(Bz,Az);%生成頻率響應(yīng)參數(shù)f1=w/pi*fs/2; %采樣頻率轉(zhuǎn)換成模擬采樣頻率y1=filter(Bz,Az,c);%使用filter函數(shù)對信號進(jìn)行濾波figure(2);subplot(211);plot(f1,angle(H));xlabel('頻率/Hz');ylabel('幅度');title('低通濾波器相頻特性 ');subplot(212);plot(f1,abs(H));xlabel('頻率/Hz');ylabel('幅度');title('低通濾波器幅頻特性 ');%低通IIR濾波后圖形figure(3);subplot(311);plot(t,y1);title('濾波后時域波形 ');xlabel('時