一维信号小波去噪原理解析及python实现方式
作者:衷科知眠
一维信号小波去噪原理及python实现
信号去噪是经常用到的信号预处理过程,以达到在保留原有信号真是信息的基础上尽可能低降低或者消除噪声,获得更高质量的信号,从而为下一步的处理奠定基础。
去噪方法可分为时域方法与频域方法。时域方法是指直接在原始信号上进行处理,比如均值滤波器、中值滤波器、EMD分解等方法。频域方法是指在信号的变换域进行去噪然后再恢复到时域得到去噪后的信号,比如小波变换、傅里叶变换等方法。
无论是一维信号还是二维信号其原理都是一样的,只不过实现上所有不同。这里主要记录一下我目前需要处理的一维信号的去噪。
噪声可分为加性噪声和乘性噪声两种形式,加性噪声一般假设其与原始信号无关,假设含噪信号为,其可表示为:
其中表示原始的干净信号,表示噪声信号。
而乘性噪声则假设噪声与信号有关,假如噪声和信号成正比,则含噪信号可表示为:
其中第二个噪声项是受原始信号影响的,原始信号的值越大,则噪声越大。
加性噪声的处理与计算都比较简单,所以一般均假设噪声为加性噪声。
基于小波变换的去噪方法其具体原理这里就不阐述了,其主要思想就是噪声在经过小波变换后的小波系数会远远小于原始信号的小波系数,因此可通过简单的阈值处理就能够去除绝大部分的噪声,同时又不会对原始信号造成很大的损失。对具体原理感兴趣的同学可以参考其他相关资料。
python 实现
使用python去噪,不需要我们自己动手实现相应的原理,安装相应的包即可完成。
首先安装小波变换的python包PyWavelets
pip install PyWavelets
安装好之后就可以直接使用,我们可通过定义去噪函数将去噪过程打包,假设去噪函数名为wavelet_denoising, 其去噪过程及参数设置代码如下:
def wavelet_denoising(signal): # signal: list wavelet_basis = 'db8' # 此处选择db8小波基 w = pywt.Wavelet(wavelet_basis) maxlevel = pywt.dwt_max_level(len(signal), w.dec_len) # 根据数据长度计算分解层数 threshold = 0.1 # threshold for filtering coeffs = pywt.wavedec(signal, wavelet_basis, level=maxlevel) # 将信号进行小波分解 # 使用阈值滤除噪声 for i in range(1, len(coeffs)): coeffs[i] = pywt.threshold(coeffs[i], threshold*max(coeffs[i])) # 将噪声滤波 # pywt.threshold 方法有三种阈值模式:硬阈值、软阈值、软硬结合 # 默认为soft 软阈值模式,其调用格式为 # pywt.threshold(signal, threshold_value, mode='soft', substitute=0) # 处理后值为 signal/np.abs(signal) * np.maximum(np.abs(signal) - threshold_value, 0) # 具体可参考pywavelets官方文档的说明 # 重建信号 signalrec= pywt.waverec(coeffs, wavelet_basis) return signalrec
通过调用wavelet_denoising方法,即可实现对原始信号的去噪。
值得注意的是去噪前后信号的长度会有所变化,若原始信号样本长度为奇数,则去噪后长度加1,若原始长度为偶数,则去噪后长度不变。
通过实验发现,如果要计算去噪前后信号的差,对于原始数据长度为奇数时,应该计算去噪后信号的第一个元素到倒数第二个元素这个子序列与原始信号之间的差值,也即舍去因去噪而多出的最后一个值。
如果使用去噪后的第二个元素到最后一个元素与原始信号的差,则可能出现偏差较大的情况。
如下图所示,其中蓝色线为去噪信号的第一个值到倒数第二个值与原始信号的差,红色线为去噪信号的第二个值到最后一个值与原始信号的差值。
小波去噪代码
先放代码,代码如下:
a = np.loadtxt('C:\\Users\\Administrator\\Desktop\\数据文件\\20201008009.txt', dtype=float)# 20201008009为一个一维列表 index = [] data = [] coffs = [] for i in range(len(a) - 1): X = float(i) Y = float(a[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.1 # Threshold for filtering # Decompose into wavelet components,to the level selected: coffs = pywt.wavedec(data, 'db8', level=maxlev) # 将信号进行小波分解 for i in range(1, len(coffs)): coffs[i] = pywt.threshold(coffs[i], threshold * max(coffs[i])) datarec = pywt.waverec(coffs, 'db8') # 将信号进行小波重构 mintime = 0 maxtime = mintime + len(data) print(mintime, maxtime) plt.rcParams['font.sans-serif']=['SimHei'] #使得标题输出为中文 plt.rcParams['axes.unicode_minus'] = False #使得标题输出为中文 plt.figure(1) #plt.subplot(3, 1, 1) plt.plot(index[mintime:maxtime], data[mintime:maxtime]) plt.title("原始信号") plt.figure(2) #plt.subplot(3, 1, 2) plt.plot(index[mintime:maxtime], datarec[mintime:maxtime]) plt.title("小波降噪信号") plt.figure(3) #plt.subplot(3, 1, 3) plt.plot(index[mintime:maxtime], data[mintime:maxtime] - datarec[mintime:maxtime]) plt.xlabel('time (s)') plt.ylabel('error (uV)') plt.tight_layout() plt.show()
去噪之前的波形
去噪之后的波形
最后一个我估计应该是去除的噪声波形
总结
以上为个人经验,希望能给大家一个参考,也希望大家多多支持脚本之家。