小波降噪与重构例子 python

原理讲解

傅里叶变换
关于傅里叶变换的基本概念在此我们就不再赘述了，
下面我们主要将傅里叶变换的不足。即我们知道傅里叶变化可以分析信号的频谱，那么为什么还要提出小波变换？答案“对非平稳过程，傅里叶变换有局限性”。看如下一个简单的信号：

做完FFT（傅里叶变换）后，可以在频谱上看到清晰的四条线，信号包含四个频率部分。一切没有问题，但是，如果是非平稳信号呢？

如上图，最上边的是频率始终不变的平稳信号。而下边两个则是频率随着时间改变的非平稳信号，它们同样包含和最上信号相同频率的四个成分。做FFT后，我们发现这三个时域上有巨大差异的信号，频谱（幅值谱）却非常一致。尤其是下边两个非平稳信号，我们从频谱上无法区分它们，因为它们包含的四个频率的信号的成分确实是一样的，只是出现的先后顺序不同。
可见，傅里叶变换处理非平稳信号有先天缺陷。它只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此，时域相差很大的两个信号，可能频谱一样。

然而平稳信号大多是制造出来的，自然界的大量信号几乎都是非平稳的。
短时傅里叶变换

一个简单可行的方法就是——加窗。 “把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现了什么频率了。”这就是短时傅里叶变换。
看图：

时域上分成一段一段做FFT，不就知道频率成分随着时间的变化情况了吗？
使用STFT存在一个问题，我们应该用多宽的窗函数？
窗太宽太窄都有问题：

小波变换
那么你可能会想到，让窗口大小变起来，多做几次STFT不就可以了吗？！没错，小波变换就有着这样的思路。
但事实上小波并不是这么做的（有人认为“小波变换就是根据算法，加不等长的窗，对每一小部分进行傅里叶变换”，这是不准确的。小波变换并没有采用窗的思想，更没有做傅里叶变换。）
至于为什么不采用可变窗的STFT呢，我认为是因为这样做冗余会太严重，STFT做不到正交化，这也是它的一大缺陷。
于是小波变换的出发点和STFT还是不同的。STFT是给信号加窗，分段做FFT；而小波直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间了~
解释：
来我们回顾一下傅里叶变换吧，没弄清傅里叶变换为什么能得到信号各个频率成分的同学也可以再借我的图理解下。
傅里叶变换把无限长的三角函数作为基函数：

这个基函数会伸缩，平移（其实是两个正交基的分解）。缩得窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。某一个尺度（宽窄）下乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。那么我们就知道信号包含该频率的成分的多少。

看，这两种尺度能乘出一个大的值（相关度高），所以信号包含较多的这两个频率成分，在频谱上这两个频率会出现两个峰。

以上，就是粗浅意义上傅里叶变换的原理。
如前边所说，小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。

从公式可以看出，不同于傅里叶变换，变量只有频率ω，小波变换有两个变量：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函数的伸缩，平移量 τ控制小波函数的平移。尺度就对应于频率（反比），平移量 τ就对应于时间。

当伸缩、平移到这么一种重合情况时，也会相乘得到一个大的值。这时候和傅里叶变换不同的是，这不仅可以知道信号有这样频率的成分，而且知道它在时域上存在的具体位置。

而当我们在每个尺度下都平移着和信号乘过一遍后，我们就知道信号在每个位置都包含哪些频率成分。

原理部分参考文献
原理部分参考来源

小波变换原理

感谢大佬，感谢大佬。

代码部分

简单例子，python版
读取数据

import matplotlib.pyplot as plt
import matplotlib; matplotlib.use('TkAgg')
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
import pywt
import pywt.data
txt = pd.read_csv('电价.csv',header=None,names=['时间','当前时刻电价'])
print(txt.head())

#小波降噪与重构画图

txt=np.array(txt)
def plot_signal_decomp(data, w,title):data = data[:, -1]  # 只对电价数据进行去噪重构mode = pywt.Modes.smooth"""Decompose and plot a signal S.S = An + Dn + Dn-1 + ... + D1"""w = pywt.Wavelet(w)#选取小波函数a = dataca = []#近似分量cd = []#细节分量for i in range(1):(a, d) = pywt.dwt(a, w, mode)#进行1阶离散小波变换ca.append(a)cd.append(d)rec_a = []rec_d = []for i, coeff in enumerate(ca):coeff_list = [coeff, None] + [None] * irec_a.append(pywt.waverec(coeff_list, w))#重构for i, coeff in enumerate(cd):coeff_list = [None, coeff] + [None] * irec_d.append(pywt.waverec(coeff_list, w))fig = plt.figure()plt.subplot(3, 1,1)plt.plot(data,label='原始电价数据曲线')plt.legend()plt.subplot(3, 1,2)plt.plot(rec_a[-1],'r', label='电价数据曲线趋势')plt.legend()plt.subplot(3, 1,3)plt.plot(rec_d[0],'g', label='电价数据曲线噪声')plt.legend()plt.show()
plot_signal_decomp(txt, 'sym5','sym5')#这里选择sym5小波，小波还有

#获取重构后的趋势部分和噪声

def choose(data, w):data=data[:,-1]#只对电价数据进行去噪重构mode = pywt.Modes.smoothw = pywt.Wavelet(w)  # 选取小波函数a = dataca = []  # 近似分量cd = []  # 细节分量for i in range(1):#1阶变换(a, d) = pywt.dwt(a, w, mode)  # 进行1阶离散小波变换 dwt为分解ca.append(a)cd.append(d)rec_a = []rec_d = []for i, coeff in enumerate(ca):coeff_list = [coeff, None] +[None] * i#填充长度rec_a.append(pywt.waverec(coeff_list, w))  # waverec重构rec_a=np.array(rec_a[-1])#-1取最后一个基波rec_a=rec_a.reshape(-1,1)for i, coeff in enumerate(cd):coeff_list = [None, coeff] + [None] * irec_d.append(pywt.waverec(coeff_list, w))#取出所有噪声波，1阶变换，有3个噪声rec_d=np.array(rec_d)return  rec_a,rec_d#取出趋势，噪声
trainingrec_a, trainingrec_d=choose(txt, 'sym5')#训练集电价趋势 和噪声

trainingrec_a为趋势
trainingrec_d为噪声

算例2 多阶变换代码


import numpy as np
import matplotlib.pyplot as plt
import matplotlib; matplotlib.use('TkAgg')
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
import pywt #导入小波包#==============生成数据==============
X=np.random.uniform(90,100,100) #生成范围90-100的数，数据100个#============小波绘图==============
def plot_signal_decomp(data, w,title):""":param data: 要进行小波分解重构的数据:param w: 小波基名称:param title: 绘图名称:return: 图"""mode = pywt.Modes.smoothw = pywt.Wavelet(w)  # 选取小波函数a = dataca = []  # 近似分量 ，趋势cd = []  # 细节分量，噪声n=3   #定义阶数，这里定义3阶for i in range(n):#3阶(a, d) = pywt.dwt(a, w, mode)  # 进行n阶离散小波变换ca.append(a)#存放趋势cd.append(d)#存储噪声rec_a = []#存放重构后的趋势rec_d = []#存放重构后的噪声for i, coeff in enumerate(ca):coeff_list = [coeff, None] + [None] * irec_a.append(pywt.waverec(coeff_list, w))  # 重构for i, coeff in enumerate(cd):coeff_list = [None, coeff] + [None] * irec_d.append(pywt.waverec(coeff_list, w))print('趋势个数len(rec_a):',len(rec_a))print('噪声个数len(rec_d):',len(rec_d))#趋势只要最后一个，噪声全要。 加原始数据个数。所以子图个数=n+1+2fig = plt.figure()for  i  in range(n+2):if i==n+1:plt.subplot(n+2, 1, i+1)plt.plot(data, label='原数据曲线')plt.legend()elif i==n:plt.subplot(n+2, 1, i+1)plt.plot(rec_a[-1], 'g', label='数据曲线趋势')plt.legend()else:plt.subplot(n+2, 1, i+1)plt.plot(rec_d[i], 'r', label='数据曲线噪声%d'%(i+1))plt.legend()plt.show()plot_signal_decomp(X, 'sym5','sym5')#这里选择sym5小波，小波还有#=======获取重构后的数据========
def choose(data, w):""":param data: 要进行小波分解重构的数据:param w: 小波基:return: 趋势 和噪声"""mode = pywt.Modes.smoothw = pywt.Wavelet(w)  # 选取小波函数a = dataca = []  # 近似分量cd = []  # 细节分量n = 3  # 定义阶数，这里定义3阶for i in range(n):  # 3阶(a, d) = pywt.dwt(a, w, mode)  # 进行n阶离散小波变换ca.append(a)  # 存放趋势cd.append(d)  # 存储噪声rec_a = []  # 存放重构后的趋势rec_d = []  # 存放重构后的噪声for i, coeff in enumerate(ca):coeff_list = [coeff, None] + [None] * irec_a.append(pywt.waverec(coeff_list, w))  # 重构for i, coeff in enumerate(cd):coeff_list = [None, coeff] + [None] * irec_d.append(pywt.waverec(coeff_list, w))#返回最后一个趋势、所有噪声 return rec_a[-1],rec_dtrainingrec_a, trainingrec_d=choose(X, 'sym5')#所有趋势 和噪声
print('趋势\n',trainingrec_a)
print('噪声1\n',trainingrec_d[0])
print('噪声2\n',trainingrec_d[1])
print('噪声3\n',trainingrec_d[2])

最后得到的是趋势和3部分噪声。。重构一般是：舍弃一部分噪声，如（保留2部分噪声）。重构结果=噪声1+噪声2+趋势。。或者只要趋势。。原式数据=噪声1+噪声2+噪声3+趋势。若不舍弃噪声，小波分解没有意义。。