本文摘要(由AI生成):
本文介绍了小波降噪方法,包括小波分解和小波降噪的基本原理,阈值的确定方法,以及软阈值和硬阈值两种处理方式。文章还提供了Python代码实现。
接上期信号处理基础之噪声与降噪(四) | 进击的EMD族降噪及python代码实现,本期为大家介绍小波降噪,并给出python代码。后续将会介绍基于DTWT和EWT的降噪应用,敬请关注。由于作者水平有限,文章的原理介绍与python代码难免有误,欢迎交流讨论!
1 小波分解与小波降噪
1.1 小波分解的基本原理
1.2 小波降噪的基本原理
1.3 小波降噪的基本过程
2.1 固定阈值原则与python实现
2.2 无偏似然估计原则与python实现
2.3 极值阈值原则与python实现
2.4 启发式阈值原则
3.1 软阈值方法
3.2 硬阈值方法
4.1 测试用例
4.2 降噪结果
4.3 降噪指标
1.1 小波分解的基本原理
小波降噪方法通过小波分解将噪声信号从原始信号中分离出来,信号经小波分解后可以得到一系列高频信号和低频信号,小波分解的过程如图1所示,记原始信号为S,第i次分解后的低频信号为
1.2 小波降噪的基本原理
注:小波变换具有线性特性,含噪信号经过小波分解后的小波系数等价于噪声的小波变换和信号的小波变换的代数和[1]。
小波降噪的关键是确定降噪准则,目前的方法较多,本文仅对阈值法降噪进行简单介绍,常用的阈值选取原则有4种[2]:固定阈值原则、无偏似然估计原则、极值阈值原则和启发式阈值原则。在阈值确定后,需采用阈值量化函数对小波系数进行处理,常用的处理方式有硬阈值函数和软阈值函数两种。
1.3 小波降噪的基本过程
import pywt
import numpy as np
def denoised_with_pywt_universal(data, wavelet='db1', mode='soft', level=1):
# 计算多级小波变换的系数
coeffs = pywt.wavedec(data, wavelet, mode='symmetric', level=level)
# 计算阈值
sigma = np.median(np.abs(coeffs[-1])) / 0.6745
threshold = sigma * np.sqrt(2 * np.log(len(data)))
# 固定阈值去噪
coeffs[1:] = (pywt.threshold(i, value=threshold, mode=mode) for i in coeffs[1:])
# 重构信号
return pywt.waverec(coeffs, wavelet, mode='symmetric')
import numpy as np
import pywt
def sure_threshold(coeffs):
n = len(coeffs)
risks = np.zeros(n)
for i in range(n):
threshold = coeffs[i]
# 计算硬阈值处理后的风险
risks[i] = n - 2 * np.sum(coeffs > threshold) + np.sum((coeffs - threshold)**2)
# 选择最小风险对应的阈值
min_risk_index = np.argmin(risks)
return coeffs[min_risk_index]
def denoised_with_pywt_sure(data, wavelet='db1', mode='soft', level=1):
# 计算多级小波变换的系数
coeffs = pywt.wavedec(data, wavelet, mode='symmetric', level=level)
# 对每一层系数应用SURE方法计算阈值并去噪
for i in range(1, len(coeffs)):
threshold = sure_threshold(np.abs(coeffs[i]))
coeffs[i] = pywt.threshold(coeffs[i], value=threshold, mode=mode)
# 重构信号
return pywt.waverec(coeffs, wavelet, mode='symmetric')
import pywt
import numpy as np
def denoised_with_pywt_minimax(data, wavelet='db1', mode='soft', level=1):
# 计算多级小波变换的系数
coeffs = pywt.wavedec(data, wavelet, mode='symmetric', level=level)
# 计算阈值
sigma = np.median(np.abs(coeffs[-1])) / 0.6745
if len(data) >= 32:
threshold = sigma * (0.3936+(0.1829*np.log(len(data))/np.log(2)))
else:
threshold = 0
# 固定阈值去噪
coeffs[1:] = (pywt.threshold(i, value=threshold, mode=mode) for i in coeffs[1:])
# 重构信号
return pywt.waverec(coeffs, wavelet, mode='symmetric')
import matplotlib.pyplot as plt
import matplotlib
# 创建一个信号
n = 500 # 生成500个点的信号
t = np.linspace(0, 30*np.pi, n, endpoint=False)
s = np.cos(0.1*np.pi*t) + np.sin(0.3*np.pi*t) + np.cos(0.5*np.pi*t) + np.sin(0.7*np.pi*t) # 原始信号
r = 0.5 * np.random.randn(n)
y = s + r # 加噪声
denoised_data_universal = denoised_with_pywt_universal(y, wavelet='db4', mode='soft', level=8)
denoised_data_sure = denoised_with_pywt_sure(y, wavelet='db4', mode='soft', level=8)
denoised_data_minimax = denoised_with_pywt_minimax(y, wavelet='db4', mode='soft', level=8)
value_universal = denoise_evaluate(s, r, denoised_data_universal)
value_sure = denoise_evaluate(s, r, denoised_data_sure)
value_minimax = denoise_evaluate(s, r, denoised_data_minimax)
fig = plt.figure(figsize=(16, 12))
plt.subplots_adjust(hspace=0.5)
plt.subplots_adjust(left=0.05, right=0.98, top=0.95, bottom=0.05)
font = {'family': 'Times New Roman', 'size': 16, 'weight': 'normal',}
matplotlib.rc('font', **font)
plt.subplot(5, 1, 1)
plt.plot(t, s, linewidth=1.5, color='b')
plt.title('pure signal', fontname='Times New Roman', fontsize=16)
plt.subplot(5, 1, 2)
plt.plot(t, y, linewidth=1.5, color='b')
plt.title('noised signal', fontname='Times New Roman', fontsize=16)
plt.subplot(5, 1, 3)
plt.plot(t, denoised_data_universal, linewidth=1.5, color='b')
plt.title('固定阈值法', fontname='Simsun', fontsize=16)
plt.subplot(5, 1, 4)
plt.plot(t, denoised_data_sure, linewidth=1.5, color='b')
plt.title('无偏似然估计阈值', fontname='Simsun', fontsize=16)
plt.subplot(5, 1, 5)
plt.plot(denoised_data_minimax, linewidth=1.5, color='b')
plt.title('极值阈值', fontname='Simsun', fontsize=16
plt.show()