该程序将使用小波多分辨分析对信号进行降噪,主要是降噪流程,为以后的小波更复杂的降噪算法打下良好的基础。降噪算法流程大致如下:
(1)去趋势项(如直流电流),并将数据归一化到区[0, 1];
(2)进行多级小波分解;
(3)使用步骤 (2)中的细节系数 cD 确定合适的阈值,给出5种不同的方法确定阈值;
(4)将简单的软阈值或硬阈值方法应用于细节系数;
(5)重建信号。
阈值确定方法,更多的细节请查看相关论文,很多
1. universal
在这种情况下,阈值由公式 MAD x sqrt{2 x log(m)} 给出,其中 MAD 是中值绝对偏差,m 是信号的长度。
2. sqtwolog
和universal一样,只是不使用MAD。
3. energy
在这种情况下,阈值算法估计细节系数的能量,并使用它们来估计最佳阈值。
4. stein
此方法实现了 Stein 的无偏风险估计。
5. heurstein
这是 Stein 的无偏风险估计的启发式实现。
首先导入相关模块
import numpy as npimport pandas as pdimport matplotlib.pylab as pltfrom scipy.signal import butter, filtfiltfrom scipy.signal import spectrogramfrom denoising import WaveletDenoising
写个函数用于绘制所有小波分解的系数
def plot_coeffs_distribution(coeffs): fig = plt.figure() size_ = int(len(coeffs) // 2) + 1 if size_ % 2 != 0: size_ = size_+1 for i in range(len(coeffs)): ax = fig.add_subplot(size_, 2, i+1) ax.hist(coeffs[i], bins=50)def pretty_plot(data, titles, palet, fs=1, length=100, nperseg=256): fig = plt.figure(figsize=(13, 13)) fig.subplots_adjust(hspace=0.5, wspace=0.5) index = 1 for i, d in enumerate(data): ax = fig.add_subplot(8, 2, index) ax.plot(d[:length], color=palet[i]) ax.set_title(titles[i]) ax = fig.add_subplot(8, 2, index+1) f, t, Sxx = spectrogram(d, fs=fs, nperseg=nperseg) ax.pcolormesh(t, f, Sxx, shading='auto') index += 2
5种阈值方法的计算函数
def run_experiment(data, level=2, fs=1, nperseg=256, length=100): titles = ['Original data', 'Universal Method', 'SURE Method', 'Energy Method', 'SQTWOLOG Method', 'Heursure Method'] experiment = ['universal', 'stein', 'energy', 'sqtwolog', 'heurstein'] wd = WaveletDenoising(normalize=False, wavelet='db3', level=level, thr_mode='soft', selected_level=level, method="universal", resolution=100, energy_perc=0.90) res = [data] for i, e in enumerate(experiment): wd.method = experiment[i] res.append(wd.fit(data)) palet = ['r', 'b', 'k', 'm', 'c', 'orange', 'g', 'y'] pretty_plot(res, titles, palet, fs=fs, length=length, nperseg=nperseg)
搞个最简单的数据,以便于初级分析
if __name__ == '__main__': data = np.zeros((128,)) data[:16] = 4 data += np.random.normal(0, 1, (128,)) run_experiment(data, level=3, length=100, nperseg=32) plt.show()
详细代码如下
https://mianbaoduo.com/o/bread/Yp2al5Zt