3篇小波分析的基础入门材料（绝对入门级！含PyWavelet实例）

引用 Rainyboy 2013-2-16 00:07

本帖最后由 Rainyboy 于 2013-2-16 00:15 编辑

自己按照那个网址的代码和说明玩了一下。顺便把代码写了个注释。


======原始信号的噪声=====


======原始信号的小波系数=====



======噪声信号的小波系数=====


======去噪之后的小波系数=====


======最终结果=====



======所有代码=====

1. # -*- coding: cp936 -*-
   
2. import pywt
   
3. import matplotlib.pyplot as plt
   
4. import numpy as np
   
5. from scipy import stats
   
6. from statsmodels.robust import stand_mad
   
7. 
   
8. wavtag = 'db8'
   
9. 
   
10. #===============================================================================
    
11. # 图1：绘出Haar小波母函数
    
12. #===============================================================================
    
13. 
    
14. #这里不是“函数调用”，二是“对象声明和创建”
    
15. #创建了一个pywt.Wavelet类，用以描述小波母函数的各种性质
    
16. w = pywt.Wavelet('Haar')
    
17. 
    
18. #调用Wavefun()成员函数，返回：
    
19. #phi - scaling function 尺度函数
    
20. #psi - wavelet function 母函数
    
21. phi,psi,x = w.wavefun(level=10)
    
22. 
    
23. #注意，此处采用“面对对象”的方式使用matplotlib
    
24. #而不是“状态机”的方式
    
25. fig = plt.figure()
    
26. ax = fig.add_subplot(111)
    
27. ax.set_xlim(-0.02,1.02)
    
28. ax.plot(x,psi)
    
29. ax.grid(True)
    
30. plt.show()
    
31. 
    
32. 
    
33. #===============================================================================
    
34. # 图2：Debauchies小波的尺度函数和母函数
    
35. #===============================================================================
    
36. 
    
37. db8 = pywt.Wavelet(wavtag)
    
38. scaling, wavelet, x = db8.wavefun()
    
39. 
    
40. fig = plt.figure(2)
    
41. ax1 = fig.add_subplot(121)
    
42. ax1.plot(x,scaling)
    
43. ax1.set_title('Scaling function,' + wavtag)
    
44. ax1.set_ylim(-1.2,1.2)
    
45. ax1.grid(True)
    
46. 
    
47. ax2 = fig.add_subplot(122,sharey = ax1)
    
48. ax2.set_title('Wavelet,'+wavtag)
    
49. ax2.plot(x,wavelet)
    
50. ax2.tick_params(labelleft=False)
    
51. ax2.grid(True)
    
52. 
    
53. plt.tight_layout()
    
54. plt.show()
    
55. 
    
56. #===============================================================================
    
57. # 图3：小波去噪模拟，原始信号和混合噪声的信号
    
58. #===============================================================================
    
59. def Blocks(x):
    
60.     K = lambda x : (1+np.sign(x))/2
    
61.     t = np.array([[0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81]]).T
    
62.     h = np.array([[4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,2.1,-4.2]]).T
    
63.     return 3.655606*np.sum(h*K(x-t), axis=0)
    
64. 
    
65. def bumps(x):
    
66.     K = lambda x : (1. + np.abs(x)) ** -4.
    
67.     t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
    
68.     h = np.array([[4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 2.1, 4.2]]).T
    
69.     w = np.array([[.005, .005, .006, .01, .01, .03, .01, .01, .005, .008, .005]]).T
    
70.     return np.sum(h*K((x-t)/w), axis=0)
    
71. 
    
72. #构造原始数据
    
73. x = np.linspace(0,1,2**15)
    
74. blk = bumps(x)
    
75. 
    
76. #构造含噪声的数据
    
77. np.random.seed(12345)
    
78. nblk = blk + stats.norm().rvs(2**15)*0.3
    
79. 
    
80. fig = plt.figure(3)
    
81. ax31 = fig.add_subplot(211)
    
82. ax31.plot(x,blk)
    
83. ax31.grid(True)
    
84. ax31.set_title('Original Data')
    
85. ax31.tick_params(labelbottom=False)
    
86. 
    
87. ax32 = fig.add_subplot(212)
    
88. ax32.plot(x,nblk)
    
89. ax32.grid(True)
    
90. ax32.set_title('Nosy Data')
    
91. 
    
92. plt.show()
    
93. 
    
94. #===============================================================================
    
95. # 图4,5：小波分析，及数据展示
    
96. #===============================================================================
    
97. def coef_pyramid_plot(coefs, first=0, scale='uniform', ax=None):
    
98.     """
    
99.     Parameters
    
100.     ----------
     
101.     coefs : array-like
     
102.         Wavelet Coefficients. Expects an iterable in order Cdn, Cdn-1, ...,
     
103.         Cd1, Cd0.
     
104.     first : int, optional
     
105.         The first level to plot.
     
106.     scale : str {'uniform', 'level'}, optional
     
107.         Scale the coefficients using the same scale or independently by
     
108.         level.
     
109.     ax : Axes, optional
     
110.         Matplotlib Axes instance
     
111. 
     
112.     Returns
     
113.     -------
     
114.     Figure : Matplotlib figure instance
     
115.         Either the parent figure of `ax` or a new pyplot.Figure instance if
     
116.         `ax` is None.
     
117.     """
     
118. 
     
119.     if ax is None:
     
120.         import matplotlib.pyplot as plt
     
121.         fig = plt.figure()
     
122.         ax = fig.add_subplot(111, axisbg='lightgrey')
     
123.     else:
     
124.         fig = ax.figure
     
125. 
     
126.     n_levels = len(coefs)
     
127.     n = 2**(n_levels - 1) # assumes periodic
     
128. 
     
129.     if scale == 'uniform':
     
130.         biggest = [np.max(np.abs(np.hstack(coefs)))] * n_levels
     
131.     else:
     
132.         # multiply by 2 so the highest bars only take up .5
     
133.         biggest = [np.max(np.abs(i))*2 for i in coefs]
     
134. 
     
135.     for i in range(first,n_levels):
     
136.         x = np.linspace(2**(n_levels - 2 - i), n - 2**(n_levels - 2 - i), 2**i)
     
137.         ymin = n_levels - i - 1 + first
     
138.         yheight = coefs[i]/biggest[i]
     
139.         ymax = yheight + ymin
     
140.         ax.vlines(x, ymin, ymax, linewidth=1.1)
     
141. 
     
142.     ax.set_xlim(0,n)
     
143.     ax.set_ylim(first - 1, n_levels)
     
144.     ax.yaxis.set_ticks(np.arange(n_levels-1,first-1,-1))
     
145.     ax.yaxis.set_ticklabels(np.arange(first,n_levels))
     
146.     ax.tick_params(top=False, right=False, direction='out', pad=6)
     
147.     ax.set_ylabel("Levels", fontsize=14)
     
148.     ax.grid(True, alpha=.85, color='white', axis='y', linestyle='-')
     
149.     ax.set_title('Wavelet Detail Coefficients', fontsize=16,
     
150.             position=(.5,1.05))
     
151.     fig.subplots_adjust(top=.89)
     
152. 
     
153.     return fig
     
154. 
     
155. 
     
156. fig = plt.figure(4)
     
157. ax4 = fig.add_subplot(111, axisbg='lightgrey')
     
158. fig = plt.figure(5)
     
159. ax5 = fig.add_subplot(111, axisbg='lightgrey')
     
160. 
     
161. #调用wavedec()函数对数据进行小波变换
     
162. #mode指定了数据补齐的方式
     
163. #‘per’指周期延拓数据
     
164. true_coefs = pywt.wavedec(blk,wavtag,level=15,mode='per')
     
165. noisy_coefs = pywt.wavedec(nblk,wavtag,level=15,mode='per')
     
166. 
     
167. #绘出‘coefficient pyramid’
     
168. #注意，这里只绘出了detail coefficients
     
169. #而没有展示approximation coefficient(s),该数据存在true_coefs[0]中
     
170. fig1 = coef_pyramid_plot(true_coefs[1:],scale = 'level',ax=ax4)
     
171. fig1.axes[0].set_title('Original Wavelet Detail Coefficients')
     
172. 
     
173. fig2 = coef_pyramid_plot(noisy_coefs[1:],scale = 'level',ax=ax5)
     
174. fig2.axes[0].set_title('Noisy Wavelet Detail Coefficients')
     
175. 
     
176. plt.show()
     
177. 
     
178. #===============================================================================
     
179. # 图6：降噪——全局阈值
     
180. # 图7：重构数据——对比效果
     
181. #===============================================================================
     
182. 
     
183. 
     
184. sigma = stand_mad(noisy_coefs[-1])
     
185. uthresh = sigma*np.sqrt(2*np.log(len(nblk)))
     
186. denoised_coefs = noisy_coefs[:]
     
187. denoised_coefs[1:] = (pywt.thresholding.soft(data,value=uthresh) for data in denoised_coefs[1:])
     
188. 
     
189. fig = plt.figure(6)
     
190. ax6 = fig.add_subplot(111, axisbg='lightgrey')
     
191. fig3 = coef_pyramid_plot(denoised_coefs[1:],scale = 'level',ax=ax6)
     
192. fig3.axes[0].set_title('Denoised Wavelet Detail Coefficients')
     
193. 
     
194. signal = pywt.waverec(denoised_coefs, wavtag, mode='per')
     
195. 
     
196. fig = plt.figure(7)
     
197. ax71 = fig.add_subplot(211)
     
198. ax71.plot(x,nblk)
     
199. ax71.grid(True)
     
200. ax71.set_title('Noisy Data')
     
201. ax71.tick_params(labelbottom=False)
     
202. 
     
203. ax72 = fig.add_subplot(212)
     
204. ax72.plot(x,signal,label = 'Denoised')
     
205. ax72.plot(x,blk,color='red',lw = 0.5, label = 'Original')
     
206. ax72.grid(True)
     
207. ax72.set_title('Denoised Data')
     
208. ax72.legend()
     
209. 
     
210. plt.show()

复制代码

引用 ChaChing 2013-2-16 22:32

Rainyboy 发表于 2013-2-16 00:07

                               
登录/注册后可看大图

自己按照那个网址的代码和说明玩了一下。顺便把代码写了个注释。

外行人仅能看结果!
从denoised与original的比较, 好像denoised的peak值都比原来的少? 是特例还是一般性?

引用 Rainyboy 2013-2-18 00:26

ChaChing 发表于 2013-2-16 22:32

                               
登录/注册后可看大图

外行人仅能看结果!
从denoised与original的比较, 好像denoised的peak值都比原来的少? 是特例还是 ...

我也是外行，哈哈！
峰值减小了这个问题恐怕是由于采用的是软阈值导致的吧（不过硬阈值也好不到哪里去……）
就是这句话：

1. denoised_coefs[1:] = (pywt.thresholding.soft(data,value=uthresh) for data in denoised_coefs[1:])

复制代码

这个问题我没怎么想明白，只能留给内行点拨一二了……

引用 hfuthj 2013-3-27 10:29

外行看的吃力呀

引用 沐雨柠檬 2013-3-27 19:35

本帖最后由 沐雨柠檬 于 2013-3-27 20:00 编辑

哎，下不了啊，体能怎么能多些？求附件123，请好心盟友发给我啊，正学呢，急需，10914439@qq.com，不胜感激！hht变换不是有虚假频率产生么？怎么解决？hht与小波变换相比，有何优缺点？
不用传了，自己解决了，需要干活赚钱才行。

引用 jayson2000 2013-4-10 13:07

感觉看不懂啊，好难啊

引用 寂寞的部落 2013-4-14 11:31

不知道小波尺度与频率到底什么关系？

引用 linyiming 2013-5-2 11:30

挺好的啊

引用 thomasten 2013-5-3 20:21

傻逼编的这程序全错的。。。无语

引用 Rainyboy 2013-5-3 22:28

thomasten 发表于 2013-5-3 20:21

                               
登录/注册后可看大图

傻逼编的这程序全错的。。。无语

兄弟指点一下，哪里有问题？

引用 worm 2013-5-30 16:55

体力值不够啊~

引用 雕金牡丹 2013-9-6 21:10

引用 哈哈一般 2013-9-14 01:09

真心不错的分享，很佩服楼主这样的人

引用 puma1 2013-9-27 14:48

找资料过来的,想下资料,却没体力,离开了

引用 和丹 2013-10-13 17:15

谢谢，学习下

引用 wangzi629 2013-11-18 00:34

引用 hanklx 2013-11-19 02:05

这个必须顶！！！

引用 珞珈山下 2013-12-6 09:26

非常感谢，好教材不容易啊

引用 忧郁的猪 2013-12-11 20:48

需要好多体力啊，新人木有体力～下次弄成一个吧～感谢楼主

查看全部评论(166)