心電信號的噪聲
EGG信號具備微弱、低幅值、低頻、隨杋性的特色,很容易被噪聲干擾,而噪聲可能來自生物體內,如呼吸、肌肉顫抖,也可能由於接觸不良而引發體外干擾。是ECG信號主要的三種噪聲爲工頻干擾、肌電干擾和基線漂移3,也是在濾波過程當中急需被抑制去除的噪聲干擾。python
- 工頻干擾:是由採集心電信號的設備周身的供電環境引發的電磁干擾,幅值低,噪聲頻率爲50Hz左右,其波形很像一個正弦信號,該噪聲經常會淹沒有用的心電信號,也會影響P波和T波的檢測。
- 肌電干擾:在心電圖採集過程當中,由於人體運動肌肉不自主顫抖形成,這種干擾無規律可言,波形形態會急速變化,頻率很高,而且分佈很廣,範圍在0-2000Hz內,能量集中在30-300Hz內,持續時間通常爲50ms,肌電干擾與心電信號會重合在一塊兒,這會致使有用的心電信號細微的變化極可能被忽視。
- 基線漂移:屬於低頻干擾,頻率分佈在0.15-0.3Hz內,因爲電極位置的滑動變化或者人體的呼吸運動形成心電信號隨時間緩慢變化而偏離正常基線位置產生基線漂移,幅度和頻率都會時刻變化着。心電信號中的PR波段和ST波段很是容易受到影響產生失真。
心電信號的預處理
小波變換(Wavelet Transform, WT)能夠進行時頻變換,是對信號進行時域以及頻域分析的最爲理想工具。本文對含噪心電信號採用基於小波變換的去噪處理方法,分爲如下3個步驟:數組
-
因爲噪聲和信號混雜在一塊兒,首先選擇一個小波基函數,因爲噪聲和信號混雜在一塊兒,因此要用小波變換對含噪心電信號進行某尺度分解獲得各尺度上的小波係數。函數
-
心電信號通過小波變換尺度分解後,幅值比較大的小波係數就是有用的信號,幅值比較小的小波係數就是噪聲,根據心電信號和夾雜噪聲的頻率分佈,對各尺度上的小波係數進行閾值處理,把小於閾值的小波係數置零或用閾值函數處理。工具
-
分別處理完小波尺度分解後的低頻係數和高頻係數,再重構信號。url
4尺度小波分解所得各尺度係數示意圖以下,9尺度小波分解能夠類比之:spa
小波係數處理的閾值函數有硬閾值和軟閾值之分。.net
- 硬閾值函數:若分解後的係數絕對值大於閾值,保證其值不變;當其小於給定的閾值時,令其爲零。
- 軟閾值函數:若分解後的係數絕對值大於閾值,令其值減去λ;當其小於給定的閾值時,令其爲零。
其中w爲原始小波係數,W爲處理後的小波係數,λ爲給定的閾值,N爲信號長度。λ的計算公式爲
$$
\lambda=\frac{median|w|\sqrt{2lnN}}{0.6745}
$$
3d
median爲中位數code
由上文分析可知,軟閾值去噪法的小波係數在連續性上優於硬閾值法,故本文采起了軟閾值法結合小波基對信號進行仿真實驗。orm
關於小波變換的介紹能夠參考這篇文章:https://blog.csdn.net/tbl1234567/article/details/52662985
代碼實戰
去噪理論大體如此,下面來實戰編寫代碼。以編號爲100的心電數據記錄爲例,上一篇文章已經經過wfdb將心電數據讀取到了record變量中,record爲一個大小爲(65000,1)的numpy數組。要對其進行小波分解,首先要經過np.flatten()將其轉換爲shape=(65000,)的數組。
record = wfdb.rdrecord('ecg_data/' + 100, channel_names=['MLII']) data = record.p_signal.flatten()
Python中使用pywt包進行信號的時頻分析,其中包括傅里葉變換、小波變換等工具。使用pywt.wavedec()將時域信號進行9尺度變換到頻域,小波基選用db5,返回值即爲各尺度係數。
# 用db5做爲小波基,對心電數據進行9尺度小波變換 coeffs = pywt.wavedec(data=data, wavelet='db5', level=9) cA9, cD9, cD8, cD7, cD6, cD5, cD4, cD3, cD2, cD1 = coeffs
在對原始信號進行分解以後,咱們能夠知道,1-2層的細節份量的能量與原始信號的高頻干擾保持一致。代表1-2層是高頻噪聲(工頻干擾以及肌電噪聲)集中的主要地方。所以咱們須要濾除D1層和D2層的細節份量,經過將其置0來達到去除的目的。而後將信號分解獲得的3~9層小波係數經過軟閾值公式對信號的閾值進行處理。pywt.threshold()函數提供了閾值濾波功能,而且默認是mode='soft'的軟閾值濾波。
threshold = (np.median(np.abs(cD1)) / 0.6745) * (np.sqrt(2 * np.log(len(cD1)))) # 將高頻信號cD一、cD2置零 cD1.fill(0) cD2.fill(0) # 將其餘中低頻信號按軟閾值公式濾波 for i in range(1, len(coeffs) - 2): coeffs[i] = pywt.threshold(coeffs[i], threshold)
最後對小波係數進行反變換,得到去噪後的信號。
rdata = pywt.waverec(coeffs=coeffs, wavelet='db5')
對預處理先後的100號信號截取前1500個信號點進行對好比下,能夠看出降噪效果仍是不錯的:
plt.figure(figsize=(20, 4)) plt.plot(data) plt.show() plt.plot(rdata) plt.show()