Python-ECG小波变换去噪学习笔记

NeoX2003

已于 2022-12-27 20:08:44 修改

阅读量2.6k

点赞数 5

文章标签： python 学习人工智能

于 2022-12-27 20:03:40 首次发布

原文链接：http://t.csdn.cn/9WSJB

版权

一，小波去噪原理：
信号产生的小波系数含有信号的重要信息，
将信号经小波分解后小波系数较大，噪声的小波系数较小，并且噪声的小波系数要小于信号的小波系数 ，通过选取一个合适的阀值，大于阀值的小波系数被认为是有信号产生的，应予以保留，小于阀值的则认为是噪声产生的，置为零从而达到去噪的目的。

小波阀值去噪的基本问题包括三个方面：小波基的选择，阀值的选择，阀值函数的选择。
(1) 小波基的选择：通常我们希望所选取的小波满足以下条件：
--正交性、高消失矩、紧支性、对称性或反对称性。
--但事实上具有上述性质的小波是不可能存在的，因为【小波是对称或反对称的只有Haar小波】，并且高消失矩与紧支性是一对矛盾，
--所以在应用的时候一般选取具有紧支的小波以及根据信号的特征来选取较为合适的小波。（具有较短的支撑区间，并且具有较高的精度的小波）

(2) 阀值的选择：直接影响去噪效果的一个重要因素就是阀值的选取，不同的阀值选取将有不同的去噪效果。
--目前主要有通用阀值（VisuShrink）、SureShrink阀值、Minimax阀值、BayesShrink阀值等。

(3) 阀值函数的选择：阀值函数是修正小波系数的规则，不同的反之函数体现了不同的处理小波系数的策略。 最常用的阀值函数有两种：
--一种是硬阀值函数，另一种是软阀值函数。还有一种介于软、硬阀值函数之间的Garrote函数。

另外，对于去噪效果好坏的评价，常用信号的信噪比（SNR）与估计信号同原始信号的均方根误差（RMSE）来判断。

补充：
--正交基通常用于向量空间的基的选择，从而使得向量空间中的向量可以被唯一地表示为基向量的线性组合
--小波基是一组函数，可以用来描述信号的时域和频域特征。通常使用小波基进行小波变换时，会对信号进行多级分解，从而得到多个不同频率成分。
--小波变换是一种信号处理方法，它能够将信号分解为不同尺度、不同频率的部分。
  在进行小波变换时，我们可以选择进行多少层分解，从而得到不同细节程度的分解结果。分解层数是对频率范围进行一定的划分,例如，在进行两层小波分解时，我们会得到一个由两层小波系数构成的分解结果。这些小波系数可以用于描述信号不同尺度、不同频率的部分。
  在进行小波变换时，能够对data数组进行小波变换的最大层数指的是，我们能够对data数组进行的最多层小波分解的层数。
  计算最大层数的方法可以使用pywt库中的dwt_max_level()函数，它会根据data数组的长度和所使用的小波的决策长度计算出最大层数。




二，在python中使用小波分析进行阈值去噪声,使用pywt.threshold函数
#coding=gbk
#使用小波分析进行阈值去噪声,使用pywt.threshold

import pywt
import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt 
import math 
 
data = np.linspace(1, 10, 10)
print(data)
# [ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
# pywt.threshold(data, value, mode, substitute) mode 模式有4种，soft, hard, greater, less; substitute是替换值
 
data_soft = pywt.threshold(data=data, value=6, mode='soft', substitute=12)
print(data_soft)
# [12. 12. 12. 12. 12.  0.  1.  2.  3.  4.] 将小于6 的值设置为12， 大于等于6 的值全部减去6
 
data_hard = pywt.threshold(data=data, value=6, mode='hard', substitute=12)
print(data_hard)
# [12. 12. 12. 12. 12.  6.  7.  8.  9. 10.] 将小于6 的值设置为12， 其余的值不变
 
data_greater = pywt.threshold(data, 6, 'greater', 12)
print(data_greater)
# [12. 12. 12. 12. 12.  6.  7.  8.  9. 10.] 将小于6 的值设置为12，大于等于阈值的值不变化
 
data_less = pywt.threshold(data, 6, 'less', 12)
print(data_less)
# [ 1.  2.  3.  4.  5.  6. 12. 12. 12. 12.] 将大于6 的值设置为12， 小于等于阈值的值不变

三，在python中使用ecg心电信号进行小波去噪实验

import matplotlib.pyplot as plt
import pywt

# Get data:
ecg = pywt.data.ecg()  # 生成心电信号
index = []
data = []
for i in range(len(ecg)-1):
    X = float(i)
    Y = float(ecg[i])
    index.append(X)
    data.append(Y)

# Create wavelet object and define parameters
w = pywt.Wavelet('db8')  # 选用Daubechies8小波
maxlev = pywt.dwt_max_level(len(data), w.dec_len)
print("maximum level is " + str(maxlev))
threshold = 0.04  # Threshold for filtering

# Decompose into wavelet components, to the level selected:
coeffs = pywt.wavedec(data, 'db8', level=maxlev)  # 将信号进行小波分解

plt.figure()
for i in range(1, len(coeffs)):
    coeffs[i] = pywt.threshold(coeffs[i], threshold*max(coeffs[i]))  # 将噪声滤波

datarec = pywt.waverec(coeffs, 'db8')  # 将信号进行小波重构

mintime = 0
maxtime = mintime + len(data) + 1

plt.figure()
plt.subplot(2, 1, 1)
plt.plot(index[mintime:maxtime], data[mintime:maxtime])
plt.xlabel('time (s)')
plt.ylabel('microvolts (uV)')
plt.title("Raw signal")
plt.subplot(2, 1, 2)
plt.plot(index[mintime:maxtime], datarec[mintime:maxtime-1])
plt.xlabel('time (s)')
plt.ylabel('microvolts (uV)')
plt.title("De-noised signal using wavelet techniques")

plt.tight_layout()
plt.show()

# coding=gbk
# 使用小波分析进行阈值去噪声,使用pywt.threshold
# 大于阀值的小波系数被认为是有信号产生的，应予以保留；小于阀值的则认为是噪声产生的，置为零从而达到去噪的目的。

import pywt
# Python 的 pywt 包是一个用于小波变换的 Python 库。它提供了用于进行小波变换的各种函数和算法，包括单小波、多小波、小波包变换、小波阈值处理等。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math

data = np.linspace(1, 10, 10)
# linspace 是 NumPy 库中的一个函数，用于创建等差数列。
# 这段代码中，np.linspace(1, 10, 10) 会创建一个包含 10 个数字的等差数列，该数列的第一项为 1，最后一项为 10，其余数字均以相等的间隔排列。
print(data)
# [ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
# pywt.threshold(data, value, mode, substitute) mode 模式有4种，soft, hard, greater, less; substitute是替换值

# ================================================================================================================
# 阀值函数的选择：阀值函数是修正小波系数的规则，不同的反之函数体现了不同的处理小波系数的策略。 最常用的阀值函数有两种：
# --一种是硬阀值函数，另一种是软阀值函数。还有一种介于软、硬阀值函数之间的Garrote函数。

data_soft = pywt.threshold(data=data, value=6, mode='soft', substitute=12)
# pywt.threshold 是 pywt 库中的一个函数，用于对数据进行小波阈值处理。函数的第一个参数是待处理的数据，第二个参数是阈值，第三个参数是处理方式（阈值函数），第四个参数是替代值。
print(data_soft)
# [12. 12. 12. 12. 12.  0.  1.  2.  3.  4.] 将小于6 的值设置为12， 大于等于6 的值全部减去6

data_hard = pywt.threshold(data=data, value=6, mode='hard', substitute=12)
print(data_hard)
# [12. 12. 12. 12. 12.  6.  7.  8.  9. 10.] 将小于6 的值设置为12， 其余的值不变

data_greater = pywt.threshold(data, 6, 'greater', 12)
print(data_greater)
# [12. 12. 12. 12. 12.  6.  7.  8.  9. 10.] 将小于6 的值设置为12，大于等于阈值的值不变化
# threshold 函数会将输入数据中大于阈值 6 的数设为 12，小于阈值的数保留不变。参数 'greater' 表示只对大于阈值的数进行处理，不对小于阈值的数进行处理。

data_less = pywt.threshold(data, 6, 'less', 12)
print(data_less)
# [ 1.  2.  3.  4.  5.  6. 12. 12. 12. 12.] 将大于6 的值设置为12， 小于等于阈值的值不变
# 它会将小于阈值的数字置零，并将大于阈值的数字替换为指定的替代值。这段代码的作用是对数据 data 进行小波阈值处理，将数据中小于 6 的数字置零，并将大于 6 的数字替换为 12。

import matplotlib.pyplot as plt
import pywt

# Get data:
ecg = pywt.data.ecg()
# 使用了pywt库中的data.ecg()函数生成的心电信号数据是一个由数值构成的数组，其中的每个数值代表对应时刻心电图的电压值。# 生成心电信号
index = []
data = []
for i in range(len(ecg)-1):
    X = float(i)
    # 数组下标从0 开始 将元素的索引（即数组下标）转换为浮点数，并赋值给变量X
    Y = float(ecg[i])
    # 将元素的值转换为浮点数，并赋值给变量Y。
    index.append(X)
    # 将X的值添加到index列表中
    data.append(Y)
    # 将Y的值添加到data列表中。

# Create wavelet object and define parameters
w = pywt.Wavelet('db8')  # 选用Daubechies8小波
# 创建小波对象w,选用Daubechies8小波
maxlev = pywt.dwt_max_level(len(data), w.dec_len)
# 使用pywt库中的dwt_max_level()函数计算出 能够对data数组进行小波变换的最大层数
print("maximum level is " + str(maxlev))
threshold = 0.04  # Threshold for filtering
# 定义了一个阈值threshold，用于后续的过滤操作。


# Decompose into wavelet components, to the level selected:
coeffs = pywt.wavedec(data, 'db8', level=maxlev)  # 将信号进行小波分解,返回值是小波系数列表
# wavedec()函数接收三个参数 返回一个数组：
# data：要进行分解的信号数组。
# 'db8'：使用的小波类型，这里使用的是Daubechies8小波。
# level=maxlev：分解的层数，这里使用的是之前计算出的最大层数。
# wavedec()函数会返回一个由小波系数构成的列表，这些小波系数可以用于描述信号不同尺度、不同频率的部分。将返回的结果赋值给变量coeffs

plt.figure()
for i in range(1, len(coeffs)):
    coeffs[i] = pywt.threshold(coeffs[i], threshold*max(coeffs[i]))  # 将噪声滤波
# 使用pywt库中的threshold()函数对元素进行阈值滤波。将滤波后的结果覆盖回coeffs列表中。
# 阈值是之前定义的threshold变量，乘上coeffs[i]中的最大值

datarec = pywt.waverec(coeffs, 'db8')  # 将信号进行小波重构
# 小波重构是将小波分解得到的小波系数合并为原信号的过程。它是小波变换的重要组成部分，与小波分解相对应。


mintime = 0
maxtime = mintime + len(data) + 1
# 用来后续确定绘制的数据范围

plt.figure()
# 使用的是原始信号数据
plt.subplot(2, 1, 1)
plt.plot(index[mintime:maxtime], data[mintime:maxtime])
plt.xlabel('time (s)')
plt.ylabel('microvolts (uV)')
plt.title("Raw signal")
# 使用的是经过小波变换和阈值滤波后得到的信号数据
plt.subplot(2, 1, 2)
plt.plot(index[mintime:maxtime], datarec[mintime:maxtime-1])
plt.xlabel('time (s)')
plt.ylabel('microvolts (uV)')
plt.title("De-noised signal using wavelet techniques")

plt.tight_layout()
plt.show()