心電信號處理算法設計-實驗要求

data4 是一段實際采樣得到的心電數據, 采樣頻率為 100Hz, 波形如下圖所示。設計相應的算法, 計算心率, 單位為: 次/分鐘。可能會用到的知識為數字濾波器的設計, 離散傅里葉變換等。要求該算法復雜度盡量簡單, 可以移植到 Cortex-M3性能的單片機上實現(本次不需要做移植, 只需要考慮算法的復雜度即可)。可選擇自己熟悉的語言進行算法設計, 建議選擇 Matlab 進行算法設計。

同樣是 data4 實際采樣得到的心電數據, 對該數據進行相應的處理, 得到標準的心電波形。要求能夠顯示完整的 PQRST 的心臟跳動過程, 能夠清楚的看到處理后的心電信號顯示 P 波、R 波和 T 波。幅值(y 軸)可以不做處理, 或者為了計算方便也可以做歸一化處理, 不做具體要求。

認真撰寫設計報告, 要求有詳細的算法設計說明, 可以選擇文字、流程圖等圖文配合對算法進行說明。若參考了學術論文, 請注明參考文獻。

準備工作:

數據備注

采樣頻率	100Hz
數據數量	1868條
Matlab版本	2015b
操作系統	XUbuntu20.04
流程圖繪制工具	LibreOffice Draw

原始信號組成:

信號頻帶范圍濾波器

心電信號	5~20HZ[1]	-
肌電噪聲	30~300Hz[1]	低通
工頻噪聲	50Hz[1]	低通
基線漂移噪聲	0.03Hz[2]	零相移濾波器

實驗總流程圖:

#########第1部分.設計相應的算法, 計算心率####################

概述:

①肌電噪聲和工頻噪聲處理

主要代碼來自[1],[1]中濾除肌電噪聲和工頻噪聲各自用了一個濾波器，

屬于多余，保留最低截止頻率的低通濾波器即可.

過濾后效果如下:

fp=10;%通帶截止
fs=15;%阻帶截止頻率

可以看到:
①50Hz工頻被濾除
②20Hz以上信號被削弱
③之所以兩邊呈現對稱都有是因為傅里葉變換有雙邊效果.
----------------------------------------------------------------------

②基線漂移噪聲處理

[1]中的零相移(分子分母各自的相位偏移理論計算上一致稱為零相移)

濾波器代碼不完整，且沒給截止角頻率和階數

[2]中式(4)為:
$H(z)=\frac{09876\mathrm{z}?0.9876}{z?0.9752}$

可知零相移濾波器的階數為1

∵基線漂移噪聲的頻率遠遠低于心電信號，

∴基線漂移濾波器的截止頻率設為比"心電信號的頻率下限"小一些即可.

----------------------------------------------------------------------
③經過上面的兩步處理，最終得到如下波形:

在上面兩圖中對比波峰，可以看到波峰所在時刻點位置對應，所以確實是零相移。
----------------------------------------------------------------------

④對比[3]中的心患波形與③中波形(為了便于比較，經過放大處理):

實驗結果[3]房撲圖

依據上述比較可以判斷:
該心電信號來源于患有"房撲"的心臟病人.
----------------------------------------------------------------------
⑤波峰檢測:

采樣時刻(s)波峰

0.0100	229.6991
0.8600	62.4065
1.7100	72.1880
2.7900	82.7408
3.9400	81.1390
5.0200	72.2794
6.1200	65.1390
7.2800	55.6109
8.4200	62.6077
9.5300	75.2257
10.6400	71.2959
11.7800	69.1398
12.8900	66.7628
13.9900	58.2763
15.1600	86.7839
16.2800	64.3633
17.4000	77.9292
18.5800	76.4417

對照③中的圖，確實是18個波峰，說明代碼輸出正確。

數據清洗:
①第一條數據是異常數據，刪除.
②第二條數據和第三條數據時間間隔小于1s，所以也要刪除。

計算心率的處理步驟:
①計算第一列的其余數據的相對于上一條數據的時間間隔(詳細細節請見[6])，得到一個序列
②將①中的序列計算加權平均,得到心跳一次的耗時
③60s/心跳一次的耗時，可以得到平均心率(每分鐘心跳多少次)。

實驗第1部分結論:
該心電信號的心率為53.35次/分

算法復雜度分析:
這個要看底層被調用的代碼有哪些地方用for循環了，
沒有來得及去看。

第1部分代碼不靠譜的地方:
使用了零相移濾波器,
實際中不可能零相移(不知道硬件中效果如何),
以及階數512(硬件開銷較大),
可能需要后面再改。

#############-第2部分-PQRST模擬#############

根據[4],目前沒有能直接處理實驗第1部分中③中結果的波形.查看了文獻[5]，都是很理想很光滑的輸入波形，不適合本例。
[4][5]的共同特點都是對光滑密集的輸入波形進行PQRST模擬。③中波形布滿噪聲，故下面試圖處理該問題。

開始觀察波形并進行分析:

∵實驗報告第1部分的③最終判定該心電信號屬于"房撲"
∴所以上述除了最高尖峰以外的部分波形中，
判定波谷屬于有效信號,予以保留,
高頻波峰視為需要去除或平滑的噪聲.

這里有個矛盾:
峰值(心拍)所處的位置，至少有兩種諧波，
一種是窄諧波，一種是寬諧波。
窄諧波(高頻):構成峰值的形狀.
寬諧波(低頻):構成心率周期波形

信號分析:
①
∵因為窄諧波和噪聲諧波寬度很接近，也就意味著兩者頻率接近。
∴如果濾除圖中標記的噪聲，那么窄諧波也會被"誤傷"濾除，“誤傷”會導致心率難以計算.
所以只能是先計算心率，再進行PQRST模擬,不可能反著來

②如果不濾除圖中標記的噪聲，那么會導致無法使用[4]中的方式進行PQRST模擬,想要進行PQRST模擬必須濾掉房撲部分的噪聲，同時不能誤傷峰值(心拍).

下面考慮幾種方案來設法平滑該信號：

可能的方案存在的問題

峰值信號增強+低通濾波	肯定不行，濾波器無視幅值，只認頻率
多項式擬合后再恢復峰值(心拍)	階數會很高，肯定不靠譜
指數平滑+峰值(心拍)恢復	①會滯后一個采樣時刻，需要恢復 ②需要恢復峰值(心拍)
小波變換+峰值(心拍)恢復	需要恢復峰值(心拍) 好處是濾噪后沒有相移
零相移低通濾波+峰值(心拍)恢復	不靠譜，硬件上不可能完美零相移(實際效果不明)

最省事的目前想到的應該是小波變換了。

要如何恢復被小波變換誤傷的峰值信號(心拍)呢？
我們觀察到波峰到波谷大概是3~4個采樣間隔,如下:

①實驗第1部分峰值(心拍)檢測時，留下了峰值(心拍)對應的時刻位置
②小波變換過濾噪聲
③通過①中的峰值(心拍)位置,在該位置±3個采樣間隔(本實驗中±0.03s)，對信號進行增強處理。增強處理辦法如下:
峰值中心增強70mv,
峰值中心±0.01s處增強45mV
峰值中心±0.02s處增強25mV
峰值中心±0.03s處增強10mV

小波變換處理后效果如下:

峰值信號(心拍)增強后效果:

為什么不消除上圖中標記的噪聲諧波？是不是實驗處理有問題?
因為您提供的是"房撲"心電信號，而不是"健康"心電信號。
如果去除了，會把心臟病誤診為心臟健康。
所以不應該去除。

實驗第2部分結論-總體效果

實驗第2部分結論-局部放大效果

---

除了上面論述的通過小波變換進行PQRST模擬,
還可以通過指數平滑進行PQRST模擬，相關實驗報告請見參考文獻[7]
可以得到與上面類似的效果.

---

附錄

原始數據、完整代碼、運行步驟:
https://gitee.com/fastsource/heart_rate_pqrst

---

Reference:
[1]基于MATLAB的心電信號預處理
[2]Filters in the ECG Signal Processing
[3]常見的14種異常心電圖的波形特點
[4]心電信號的PQRST模擬matlab代碼(轉載+自己調研匯總)
[5]傅里葉級數在心電信號模擬中的應用
[6]根據心電信號計算心率的matlab代碼
[7]基于指數平滑對心電信號進行PQRST模擬

總結

以上是生活随笔為你收集整理的Matlab心电信号的PQRST模拟-实验报告的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。