3篇小波分析的基础入门材料（绝对入门级！含PyWavelet实例）
本帖最后由 Rainyboy 于 2014-2-10 10:54 编辑见附件，这些文献都没有从艰深的数学公式来讲述， 而是从尽量采用能够直觉到的，可以类比和想象的方式来讲述小波分析。
我觉得只要有傅立叶变换等信号处理的基本常识，都能够明白这些论文的大部分内容。
除了这三篇文献，还有一个网址结合Python中的PyWavelet库做了一个小波降噪的实例：
Wavelet Regression in Python
==补一张图，@2013-3-27==
====关于PyWavelet====
pyWavelet已经被最新的Python(xy)套件包括在其中了，当然也包含了它的文档。
Python(xy)的网址：
http://code.google.com/p/pythonxy/
本帖最后由 Rainyboy 于 2013-2-16 00:15 编辑
自己按照那个网址的代码和说明玩了一下。顺便把代码写了个注释。
======原始信号的噪声=====
======原始信号的小波系数=====
======噪声信号的小波系数=====
======去噪之后的小波系数=====
======最终结果=====
======所有代码=====
# -*- coding: cp936 -*-
import pywt
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from statsmodels.robust import stand_mad
wavtag = 'db8'
#===============================================================================
# 图1：绘出Haar小波母函数
#===============================================================================
#这里不是“函数调用”，二是“对象声明和创建”
#创建了一个pywt.Wavelet类，用以描述小波母函数的各种性质
w = pywt.Wavelet('Haar')
#调用Wavefun()成员函数，返回：
#phi - scaling function 尺度函数
#psi - wavelet function 母函数
phi,psi,x = w.wavefun(level=10)
#注意，此处采用“面对对象”的方式使用matplotlib
#而不是“状态机”的方式
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlim(-0.02,1.02)
ax.plot(x,psi)
ax.grid(True)
plt.show()
#===============================================================================
# 图2：Debauchies小波的尺度函数和母函数
#===============================================================================
db8 = pywt.Wavelet(wavtag)
scaling, wavelet, x = db8.wavefun()
fig = plt.figure(2)
ax1 = fig.add_subplot(121)
ax1.plot(x,scaling)
ax1.set_title('Scaling function,' + wavtag)
ax1.set_ylim(-1.2,1.2)
ax1.grid(True)
ax2 = fig.add_subplot(122,sharey = ax1)
ax2.set_title('Wavelet,'+wavtag)
ax2.plot(x,wavelet)
ax2.tick_params(labelleft=False)
ax2.grid(True)
plt.tight_layout()
plt.show()
#===============================================================================
# 图3：小波去噪模拟，原始信号和混合噪声的信号
#===============================================================================
def Blocks(x):
K = lambda x : (1+np.sign(x))/2
t = np.array([]).T
h = np.array([]).T
return 3.655606*np.sum(h*K(x-t), axis=0)
def bumps(x):
K = lambda x : (1. + np.abs(x)) ** -4.
t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
h = np.array([]).T
w = np.array([[.005, .005, .006, .01, .01, .03, .01, .01, .005, .008, .005]]).T
return np.sum(h*K((x-t)/w), axis=0)
#构造原始数据
x = np.linspace(0,1,2**15)
blk = bumps(x)
#构造含噪声的数据
np.random.seed(12345)
nblk = blk + stats.norm().rvs(2**15)*0.3
fig = plt.figure(3)
ax31 = fig.add_subplot(211)
ax31.plot(x,blk)
ax31.grid(True)
ax31.set_title('Original Data')
ax31.tick_params(labelbottom=False)
ax32 = fig.add_subplot(212)
ax32.plot(x,nblk)
ax32.grid(True)
ax32.set_title('Nosy Data')
plt.show()
#===============================================================================
# 图4,5：小波分析，及数据展示
#===============================================================================
def coef_pyramid_plot(coefs, first=0, scale='uniform', ax=None):
"""
Parameters
----------
coefs : array-like
Wavelet Coefficients. Expects an iterable in order Cdn, Cdn-1, ...,
Cd1, Cd0.
first : int, optional
The first level to plot.
scale : str {'uniform', 'level'}, optional
Scale the coefficients using the same scale or independently by
level.
ax : Axes, optional
Matplotlib Axes instance
Returns
-------
Figure : Matplotlib figure instance
Either the parent figure of `ax` or a new pyplot.Figure instance if
`ax` is None.
"""
if ax is None:
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111, axisbg='lightgrey')
else:
fig = ax.figure
n_levels = len(coefs)
n = 2**(n_levels - 1) # assumes periodic
if scale == 'uniform':
biggest = * n_levels
else:
# multiply by 2 so the highest bars only take up .5
biggest =
for i in range(first,n_levels):
x = np.linspace(2**(n_levels - 2 - i), n - 2**(n_levels - 2 - i), 2**i)
ymin = n_levels - i - 1 + first
yheight = coefs/biggest
ymax = yheight + ymin
ax.vlines(x, ymin, ymax, linewidth=1.1)
ax.set_xlim(0,n)
ax.set_ylim(first - 1, n_levels)
ax.yaxis.set_ticks(np.arange(n_levels-1,first-1,-1))
ax.yaxis.set_ticklabels(np.arange(first,n_levels))
ax.tick_params(top=False, right=False, direction='out', pad=6)
ax.set_ylabel("Levels", fontsize=14)
ax.grid(True, alpha=.85, color='white', axis='y', linestyle='-')
ax.set_title('Wavelet Detail Coefficients', fontsize=16,
position=(.5,1.05))
fig.subplots_adjust(top=.89)
return fig
fig = plt.figure(4)
ax4 = fig.add_subplot(111, axisbg='lightgrey')
fig = plt.figure(5)
ax5 = fig.add_subplot(111, axisbg='lightgrey')
#调用wavedec()函数对数据进行小波变换
#mode指定了数据补齐的方式
#‘per’指周期延拓数据
true_coefs = pywt.wavedec(blk,wavtag,level=15,mode='per')
noisy_coefs = pywt.wavedec(nblk,wavtag,level=15,mode='per')
#绘出‘coefficient pyramid’
#注意，这里只绘出了detail coefficients
#而没有展示approximation coefficient(s),该数据存在true_coefs中
fig1 = coef_pyramid_plot(true_coefs,scale = 'level',ax=ax4)
fig1.axes.set_title('Original Wavelet Detail Coefficients')
fig2 = coef_pyramid_plot(noisy_coefs,scale = 'level',ax=ax5)
fig2.axes.set_title('Noisy Wavelet Detail Coefficients')
plt.show()
#===============================================================================
# 图6：降噪——全局阈值
# 图7：重构数据——对比效果
#===============================================================================
sigma = stand_mad(noisy_coefs[-1])
uthresh = sigma*np.sqrt(2*np.log(len(nblk)))
denoised_coefs = noisy_coefs[:]
denoised_coefs = (pywt.thresholding.soft(data,value=uthresh) for data in denoised_coefs)
fig = plt.figure(6)
ax6 = fig.add_subplot(111, axisbg='lightgrey')
fig3 = coef_pyramid_plot(denoised_coefs,scale = 'level',ax=ax6)
fig3.axes.set_title('Denoised Wavelet Detail Coefficients')
signal = pywt.waverec(denoised_coefs, wavtag, mode='per')
fig = plt.figure(7)
ax71 = fig.add_subplot(211)
ax71.plot(x,nblk)
ax71.grid(True)
ax71.set_title('Noisy Data')
ax71.tick_params(labelbottom=False)
ax72 = fig.add_subplot(212)
ax72.plot(x,signal,label = 'Denoised')
ax72.plot(x,blk,color='red',lw = 0.5, label = 'Original')
ax72.grid(True)
ax72.set_title('Denoised Data')
ax72.legend()
plt.show()
Rainyboy 发表于 2013-2-16 00:07 static/image/common/back.gif
自己按照那个网址的代码和说明玩了一下。顺便把代码写了个注释。
外行人仅能看结果!{:{33}:}
从denoised与original的比较, 好像denoised的peak值都比原来的少? 是特例还是一般性? ChaChing 发表于 2013-2-16 22:32 static/image/common/back.gif
外行人仅能看结果!
从denoised与original的比较, 好像denoised的peak值都比原来的少? 是特例还是 ...
我也是外行，哈哈！
峰值减小了这个问题恐怕是由于采用的是软阈值导致的吧（不过硬阈值也好不到哪里去……）
就是这句话：
denoised_coefs = (pywt.thresholding.soft(data,value=uthresh) for data in denoised_coefs)
这个问题我没怎么想明白，只能留给内行点拨一二了…… 外行看的吃力呀 本帖最后由 沐雨柠檬 于 2013-3-27 20:00 编辑
哎，下不了啊，体能怎么能多些？求附件123，请好心盟友发给我啊，正学呢，急需，10914439@qq.com，不胜感激！hht变换不是有虚假频率产生么？怎么解决？hht与小波变换相比，有何优缺点？
不用传了，自己解决了，需要干活赚钱才行。{:3_61:}
感觉看不懂啊，好难啊 不知道小波尺度与频率到底什么关系？ 挺好的啊 傻逼编的这程序全错的。。。无语 thomasten 发表于 2013-5-3 20:21 static/image/common/back.gif
傻逼编的这程序全错的。。。无语
兄弟指点一下，哪里有问题？ 体力值不够啊~ {:{10}:} 真心不错的分享，很佩服楼主这样的人
找资料过来的,想下资料,却没体力,离开了