小波降噪基础-python版本

这篇小文将使用小波多分辨分析对一个简单信号进行降噪，主要是降噪流程，为以后的小波更复杂的降噪算法打下良好的基础。降噪算法流程大致如下：

（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 np
import pandas as pd
import matplotlib.pylab as plt
from scipy.signal import butter, filtfilt
from scipy.signal import spectrogram
from denoising import WaveletDenoising

写个函数用于绘制所有小波分解的系数

def plot_coeffs_distribution(coeffs):fig = plt.figure()size_ = int(len(coeffs) // 2) + 1if size_ % 2 != 0:size_ = size_+1for 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 = 1for 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] = 4data += np.random.normal(0, 1, (128,))run_experiment(data, level=3, length=100, nperseg=32)plt.show()