ECG一维信号去噪方法 python实现总结

对于心电信号的预处理第一步一般都是去噪处理，但是很多论文对于这一步都只是简单带过，为了复现论文所述方法，我感觉走了很多弯路，这里总结一下现在有做出来的一些方法，包括有中值滤波，FIR滤波，butter滤波和小波滤波。

中值滤波实现

中值滤波去除基线漂移应该是最常见的一个方法，去除基线噪声的一个常用的方法就是，用200ms和600ms的中值滤波器处理。

import matplotlib.pyplot as plt
import pywt
import numpy as np
import os
import numpy as np
from scipy.signal import medfilt
import scipy.io as io
from scipy import signal

loadData=np.load('/mnt/data/acf/Bdata/Normal.npy',allow_pickle=True)#读取数据
ecg=loadData[0]
ecg1=ecg[:,1]#取通道II

def normalize(data):
    data = data.astype('float')
    mx = np.max(data, axis=0).astype(np.float64)
    mn = np.min(data, axis=0).astype(np.float64)
    # Workaround to solve the problem of ZeroDivisionError
    return np.true_divide(data - mn, mx - mn, out=np.zeros_like(data - mn), where=(mx - mn)!=0)
ecg1=normalize(ecg1)#归一化处理
t1=int(0.2*500)
t2=int(0.6*500)
ecg2=medfilt(ecg1,t1+1)
ecg3=medfilt(ecg2,t2+1)#分别用200ms和600ms的中值滤波得到基线
ecg4=ecg1-ecg3#得到基线滤波的结果

• 之前没有在滤波之前进行归一化处理，所以得到的结果很奇怪，所以个人认为normalize的步骤是很重要的
• 使用medfilt可以进行一维信号中值滤波

plt.figure(figsize=(20, 4))
plt.plot(ecg1)#输出原图像
plt.show()
plt.figure(figsize=(20, 4))
plt.plot(ecg3)#输出基线轮廓
plt.show()
plt.figure(figsize=(20, 4))
plt.plot(ecg4)#基线滤波结果
plt.show()

FIR滤波实现

参考官网用法firwin
包括了高通FIR，低通FIR，带通FIR，带阻FIR滤波
这里举个用法例子

import matplotlib.pyplot as plt
import pywt
import numpy as np
import os
import numpy as np
from scipy.signal import medfilt
import scipy.io as io
from scipy import signal
loadData=np.load('/mnt/data/acf/Bdata/Normal.npy',allow_pickle=True)#读取数据
ecg=loadData[0]
ecg1=ecg[:,1]#取通道II

#0.67Hz高通滤波
b = signal.firwin(51,0.67,pass_zero=False,fs=500)    # FIR高通滤波
ecg3 = signal.lfilter(b, 1, ecg1)
plt.figure(figsize=(20, 4))
plt.plot(ecg1)
plt.plot(ecg3)
plt.show()

巴特沃斯滤波

参考博客利用Python scipy.signal.filtfilt() 实现信号滤波
主要利用的是scipy.signal.butter函数

import matplotlib.pyplot as plt
import pywt
import numpy as np
import os
import numpy as np
from scipy.signal import medfilt
import scipy.io as io
loadData=np.load('/mnt/data/acf/Bdata/LVH.npy',allow_pickle=True)
ecg=loadData[0]
ecg1=ecg[:,1]
def normalize(data):
    data = data.astype('float')
    mx = np.max(data, axis=0).astype(np.float64)
    mn = np.min(data, axis=0).astype(np.float64)
    # Workaround to solve the problem of ZeroDivisionError
    return np.true_divide(data - mn, mx - mn, out=np.zeros_like(data - mn), where=(mx - mn)!=0)
ecg=normalize(ecg1)

#低通和高通滤波
from scipy import signal
b, a = signal.butter(8, 0.1, 'lowpass')   #配置滤波器 8 表示滤波器的阶数
#b, a = signal.butter(8, [0.01,0.4], 'bandpass')
ecg_1 = signal.filtfilt(b, a, ecg)  #data为要过滤的信号
b, a = signal.butter(8, 0.007, 'highpass')
ecg_outline_data = signal.filtfilt(b, a,ecg_1 ) 
plt.figure(figsize=(20, 4))
#plt.plot(rdata[:,1])
plt.plot(ecg_outline_data)
plt.show()
io.savemat('Normal.mat', {'Normal': ecg_outline_data})

小波滤波

参考博客小波去噪和小波阈值选择原理

去除基线漂移使用的是小波多分辨率分解，下面先进行小波多分辨率分解适宜。

#小波多分辨率分解
def normalize(data):
    data = data.astype('float')
    mx = np.max(data, axis=0).astype(np.float64)
    mn = np.min(data, axis=0).astype(np.float64)
    # Workaround to solve the problem of ZeroDivisionError
    return np.true_divide(data - mn, mx - mn, out=np.zeros_like(data - mn), where=(mx - mn)!=0)

if __name__ == "__main__":
    import pywt
    import numpy as np
    import pylab as plt
    import math as m

    import scipy.io.wavfile as wav

    #S = 6*T +np.cos(8*np.pi**T)+0.5*np.cos(40*np.pi*T)
    #print(S)
    loadData=np.load('/mnt/data/acf/Bdata/LVH.npy',allow_pickle=True)
    S1 = loadData[0][:,1]

    wavelet='db5'
    X = range(len(S1))
    w = pywt.Wavelet('db5') # 选用Daubechies6小波
    maxlev = pywt.dwt_max_level(len(S1), w.dec_len)#level:  分解阶次。使用dwt_max_level（）时，计算信号能达到的最高分解阶次。
    print("maximum level is " + str(maxlev))
    wave = pywt.wavedec(S1, 'db5', level=8) # 将信号进行小波分解
   #小波重构
    ya8 = pywt.waverec(np.multiply(wave,[1, 0, 0, 0, 0,0,0,0,0]).tolist(),wavelet)
    yd8 = pywt.waverec(np.multiply(wave, [0, 1, 0, 0, 0,0,0,0,0]).tolist(), wavelet)
    yd7 = pywt.waverec(np.multiply(wave, [0, 0, 1, 0, 0,0,0,0,0]).tolist(), wavelet)
    yd6 = pywt.waverec(np.multiply(wave, [0, 0, 0, 1, 0,0,0,0,0]).tolist(), wavelet)
    yd5 = pywt.waverec(np.multiply(wave, [0, 0, 0, 0, 1,0,0,0,0]).tolist(), wavelet)
    yd4 = pywt.waverec(np.multiply(wave, [0, 0, 0, 0, 0,1,0,0,0]).tolist(), wavelet)
    yd3 = pywt.waverec(np.multiply(wave, [0, 0, 0, 0, 0,0,1,0,0]).tolist(), wavelet)
    yd2 = pywt.waverec(np.multiply(wave, [0, 0, 0, 0, 0,0,0,1,0]).tolist(), wavelet)
    yd1 = pywt.waverec(np.multiply(wave, [0, 0, 0, 0, 0,0,0,0,1]).tolist(), wavelet)
    P = pywt.waverec(np.multiply(wave, [0, 1, 1, 1, 1,1,1,1,1]).tolist(), wavelet)
    plt.figure(figsize=(20, 4))
    plt.plot(X, S1)
    plt.title('original signal')
    plt.figure(figsize=(20, 4))
    plt.subplot(911)
    plt.plot(X, ya8)
    plt.title('approximated component in level 8')
    plt.subplot(912)
    plt.plot(X, yd1)
    plt.title('detailed component in level 1')
    plt.subplot(913)
    plt.plot(X, yd2)
    plt.title('detailed component in level 2')    
    plt.subplot(914)
    plt.plot(X, yd3)
    plt.title('detailed component in level 3')
    plt.subplot(915)
    plt.plot(X, yd4)
    plt.title('detailed component in level 4')    
    plt.subplot(916)
    plt.plot(X, yd5)
    plt.title('detailed component in level 5')   
    plt.subplot(917)
    plt.plot(X, yd6)
    plt.title('detailed component in level 6')   
    plt.subplot(918)
    plt.plot(X, yd7)
    plt.title('detailed component in level 7')   
    plt.subplot(919)
    plt.plot(X, yd8)
    plt.title('detailed component in level 8') 
    plt.figure(figsize=(20, 4))
    plt.plot(X, P)
    plt.title('P wave')    

去除approximated component就可以去除基线漂移。
接下来使用软硬阈值结合的方法去噪。

#模块调用
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
import pywt 

#封装成函数
def sgn(num):
    if(num > 0.0):
        return 1.0
    elif(num == 0.0):
        return 0.0
    else:
        return -1.0

def wavelet_noising(new_df):
    data = new_df
    data = data.T.tolist()  # 将np.ndarray()转为列表
    w = pywt.Wavelet('db8')
    # [ca3, cd3, cd2, cd1] = pywt.wavedec(data, w, level=3)  # 分解波
    [ca5, cd5, cd4, cd3, cd2, cd1] = pywt.wavedec(data, w, level=5)  # 分解波

    length1 = len(cd1)
    length0 = len(data)

    Cd1 = np.array(cd1)
    abs_cd1 = np.abs(Cd1)
    median_cd1 = np.median(abs_cd1)

    sigma = (1.0 / 0.6745) * median_cd1
    lamda = sigma * math.sqrt(2.0 * math.log(float(length0 ), math.e))
    usecoeffs = []
    usecoeffs.append(ca5)  # 向列表末尾添加对象

    #软硬阈值折中的方法
    a = 0.5

    for k in range(length1):
        if (abs(cd1[k]) >= lamda):
            cd1[k] = sgn(cd1[k]) * (abs(cd1[k]) - a * lamda)
        else:
            cd1[k] = 0.0

    length2 = len(cd2)
    for k in range(length2):
        if (abs(cd2[k]) >= lamda):
            cd2[k] = sgn(cd2[k]) * (abs(cd2[k]) - a * lamda)
        else:
            cd2[k] = 0.0

    length3 = len(cd3)
    for k in range(length3):
        if (abs(cd3[k]) >= lamda):
            cd3[k] = sgn(cd3[k]) * (abs(cd3[k]) - a * lamda)
        else:
            cd3[k] = 0.0

    length4 = len(cd4)
    for k in range(length4):
        if (abs(cd4[k]) >= lamda):
            cd4[k] = sgn(cd4[k]) * (abs(cd4[k]) - a * lamda)
        else:
            cd4[k] = 0.0

    length5 = len(cd5)
    for k in range(length5):
        if (abs(cd5[k]) >= lamda):
            cd5[k] = sgn(cd5[k]) * (abs(cd5[k]) - a * lamda)
        else:
            cd5[k] = 0.0

    usecoeffs.append(cd5)
    usecoeffs.append(cd4)
    usecoeffs.append(cd3)
    usecoeffs.append(cd2)
    usecoeffs.append(cd1)
    recoeffs = pywt.waverec(usecoeffs, w)
    return recoeffs

#主函数
loadData=np.load('/mnt/data/acf/Bdata/Normal.npy',allow_pickle=True)
data = loadData[0][:,1]
plt.figure(figsize=(20, 4))
plt.plot(data)
plt.show()
print(data)

data_denoising = wavelet_noising(data)#调用小波去噪函数
plt.figure(figsize=(20, 4))
plt.plot(data_denoising)#显示去噪结果    l