查看: 18587|回复: 165

[小波] 3篇小波分析的基础入门材料(绝对入门级!含PyWavelet实例)

  [复制链接]
发表于 2013-2-15 12:46 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
本帖最后由 Rainyboy 于 2014-2-10 10:54 编辑

见附件,这些文献都没有从艰深的数学公式来讲述, 而是从尽量采用能够直觉到的,可以类比和想象的方式来讲述小波分析。
我觉得只要有傅立叶变换等信号处理的基本常识,都能够明白这些论文的大部分内容。

除了这三篇文献,还有一个网址结合Python中的PyWavelet库做了一个小波降噪的实例:
Wavelet Regression in Python
block_signal.png
==补一张图,@2013-3-27==
2.png


====关于PyWavelet====
pyWavelet已经被最新的Python(xy)套件包括在其中了,当然也包含了它的文档。
Python(xy)的网址:
http://code.google.com/p/pythonxy/






2_IEEE_An_Introduction_to_Wavelets.rar

1.01 MB, 下载次数: 1042

1_An_Introduction_to _Wavelets_and_some_Applications.rar

325.12 KB, 下载次数: 721

0_HP_An_Introduction_to_Wavelets.rar

1002.63 KB, 下载次数: 744

点评

赞成: 5.0
赞成: 5
好贴!!  发表于 2014-3-30 22:19
赞成: 5
真好!  发表于 2013-2-15 15:37

评分

2

查看全部评分

 楼主| 发表于 2013-2-16 00:07 | 显示全部楼层
本帖最后由 Rainyboy 于 2013-2-16 00:15 编辑

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


======原始信号的噪声=====
3.png

======原始信号的小波系数=====
4.png


======噪声信号的小波系数=====
5.png

======去噪之后的小波系数=====
6.png

======最终结果=====
7.png


======所有代码=====
  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. wavtag = 'db8'

  8. #===============================================================================
  9. # 图1:绘出Haar小波母函数
  10. #===============================================================================

  11. #这里不是“函数调用”,二是“对象声明和创建”
  12. #创建了一个pywt.Wavelet类,用以描述小波母函数的各种性质
  13. w = pywt.Wavelet('Haar')

  14. #调用Wavefun()成员函数,返回:
  15. #phi - scaling function 尺度函数
  16. #psi - wavelet function 母函数
  17. phi,psi,x = w.wavefun(level=10)

  18. #注意,此处采用“面对对象”的方式使用matplotlib
  19. #而不是“状态机”的方式
  20. fig = plt.figure()
  21. ax = fig.add_subplot(111)
  22. ax.set_xlim(-0.02,1.02)
  23. ax.plot(x,psi)
  24. ax.grid(True)
  25. plt.show()


  26. #===============================================================================
  27. # 图2:Debauchies小波的尺度函数和母函数
  28. #===============================================================================

  29. db8 = pywt.Wavelet(wavtag)
  30. scaling, wavelet, x = db8.wavefun()

  31. fig = plt.figure(2)
  32. ax1 = fig.add_subplot(121)
  33. ax1.plot(x,scaling)
  34. ax1.set_title('Scaling function,' + wavtag)
  35. ax1.set_ylim(-1.2,1.2)
  36. ax1.grid(True)

  37. ax2 = fig.add_subplot(122,sharey = ax1)
  38. ax2.set_title('Wavelet,'+wavtag)
  39. ax2.plot(x,wavelet)
  40. ax2.tick_params(labelleft=False)
  41. ax2.grid(True)

  42. plt.tight_layout()
  43. plt.show()

  44. #===============================================================================
  45. # 图3:小波去噪模拟,原始信号和混合噪声的信号
  46. #===============================================================================
  47. def Blocks(x):
  48.     K = lambda x : (1+np.sign(x))/2
  49.     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
  50.     h = np.array([[4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,2.1,-4.2]]).T
  51.     return 3.655606*np.sum(h*K(x-t), axis=0)

  52. def bumps(x):
  53.     K = lambda x : (1. + np.abs(x)) ** -4.
  54.     t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
  55.     h = np.array([[4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 2.1, 4.2]]).T
  56.     w = np.array([[.005, .005, .006, .01, .01, .03, .01, .01, .005, .008, .005]]).T
  57.     return np.sum(h*K((x-t)/w), axis=0)

  58. #构造原始数据
  59. x = np.linspace(0,1,2**15)
  60. blk = bumps(x)

  61. #构造含噪声的数据
  62. np.random.seed(12345)
  63. nblk = blk + stats.norm().rvs(2**15)*0.3

  64. fig = plt.figure(3)
  65. ax31 = fig.add_subplot(211)
  66. ax31.plot(x,blk)
  67. ax31.grid(True)
  68. ax31.set_title('Original Data')
  69. ax31.tick_params(labelbottom=False)

  70. ax32 = fig.add_subplot(212)
  71. ax32.plot(x,nblk)
  72. ax32.grid(True)
  73. ax32.set_title('Nosy Data')

  74. plt.show()

  75. #===============================================================================
  76. # 图4,5:小波分析,及数据展示
  77. #===============================================================================
  78. def coef_pyramid_plot(coefs, first=0, scale='uniform', ax=None):
  79.     """
  80.     Parameters
  81.     ----------
  82.     coefs : array-like
  83.         Wavelet Coefficients. Expects an iterable in order Cdn, Cdn-1, ...,
  84.         Cd1, Cd0.
  85.     first : int, optional
  86.         The first level to plot.
  87.     scale : str {'uniform', 'level'}, optional
  88.         Scale the coefficients using the same scale or independently by
  89.         level.
  90.     ax : Axes, optional
  91.         Matplotlib Axes instance

  92.     Returns
  93.     -------
  94.     Figure : Matplotlib figure instance
  95.         Either the parent figure of `ax` or a new pyplot.Figure instance if
  96.         `ax` is None.
  97.     """

  98.     if ax is None:
  99.         import matplotlib.pyplot as plt
  100.         fig = plt.figure()
  101.         ax = fig.add_subplot(111, axisbg='lightgrey')
  102.     else:
  103.         fig = ax.figure

  104.     n_levels = len(coefs)
  105.     n = 2**(n_levels - 1) # assumes periodic

  106.     if scale == 'uniform':
  107.         biggest = [np.max(np.abs(np.hstack(coefs)))] * n_levels
  108.     else:
  109.         # multiply by 2 so the highest bars only take up .5
  110.         biggest = [np.max(np.abs(i))*2 for i in coefs]

  111.     for i in range(first,n_levels):
  112.         x = np.linspace(2**(n_levels - 2 - i), n - 2**(n_levels - 2 - i), 2**i)
  113.         ymin = n_levels - i - 1 + first
  114.         yheight = coefs[i]/biggest[i]
  115.         ymax = yheight + ymin
  116.         ax.vlines(x, ymin, ymax, linewidth=1.1)

  117.     ax.set_xlim(0,n)
  118.     ax.set_ylim(first - 1, n_levels)
  119.     ax.yaxis.set_ticks(np.arange(n_levels-1,first-1,-1))
  120.     ax.yaxis.set_ticklabels(np.arange(first,n_levels))
  121.     ax.tick_params(top=False, right=False, direction='out', pad=6)
  122.     ax.set_ylabel("Levels", fontsize=14)
  123.     ax.grid(True, alpha=.85, color='white', axis='y', linestyle='-')
  124.     ax.set_title('Wavelet Detail Coefficients', fontsize=16,
  125.             position=(.5,1.05))
  126.     fig.subplots_adjust(top=.89)

  127.     return fig


  128. fig = plt.figure(4)
  129. ax4 = fig.add_subplot(111, axisbg='lightgrey')
  130. fig = plt.figure(5)
  131. ax5 = fig.add_subplot(111, axisbg='lightgrey')

  132. #调用wavedec()函数对数据进行小波变换
  133. #mode指定了数据补齐的方式
  134. #‘per’指周期延拓数据
  135. true_coefs = pywt.wavedec(blk,wavtag,level=15,mode='per')
  136. noisy_coefs = pywt.wavedec(nblk,wavtag,level=15,mode='per')

  137. #绘出‘coefficient pyramid’
  138. #注意,这里只绘出了detail coefficients
  139. #而没有展示approximation coefficient(s),该数据存在true_coefs[0]中
  140. fig1 = coef_pyramid_plot(true_coefs[1:],scale = 'level',ax=ax4)
  141. fig1.axes[0].set_title('Original Wavelet Detail Coefficients')

  142. fig2 = coef_pyramid_plot(noisy_coefs[1:],scale = 'level',ax=ax5)
  143. fig2.axes[0].set_title('Noisy Wavelet Detail Coefficients')

  144. plt.show()

  145. #===============================================================================
  146. # 图6:降噪——全局阈值
  147. # 图7:重构数据——对比效果
  148. #===============================================================================


  149. sigma = stand_mad(noisy_coefs[-1])
  150. uthresh = sigma*np.sqrt(2*np.log(len(nblk)))
  151. denoised_coefs = noisy_coefs[:]
  152. denoised_coefs[1:] = (pywt.thresholding.soft(data,value=uthresh) for data in denoised_coefs[1:])

  153. fig = plt.figure(6)
  154. ax6 = fig.add_subplot(111, axisbg='lightgrey')
  155. fig3 = coef_pyramid_plot(denoised_coefs[1:],scale = 'level',ax=ax6)
  156. fig3.axes[0].set_title('Denoised Wavelet Detail Coefficients')

  157. signal = pywt.waverec(denoised_coefs, wavtag, mode='per')

  158. fig = plt.figure(7)
  159. ax71 = fig.add_subplot(211)
  160. ax71.plot(x,nblk)
  161. ax71.grid(True)
  162. ax71.set_title('Noisy Data')
  163. ax71.tick_params(labelbottom=False)

  164. ax72 = fig.add_subplot(212)
  165. ax72.plot(x,signal,label = 'Denoised')
  166. ax72.plot(x,blk,color='red',lw = 0.5, label = 'Original')
  167. ax72.grid(True)
  168. ax72.set_title('Denoised Data')
  169. ax72.legend()

  170. plt.show()
复制代码




点评

赞成: 4.0
赞成: 4
  发表于 2014-6-25 10:16
赞成: 4
  发表于 2014-5-13 09:28
赞成: 4
  发表于 2014-3-30 22:19

评分

1

查看全部评分

回复 支持 2 反对 0

使用道具 举报

发表于 2013-2-16 22:32 | 显示全部楼层

外行人仅能看结果!
从denoised与original的比较, 好像denoised的peak值都比原来的少? 是特例还是一般性?
 楼主| 发表于 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:])
复制代码

这个问题我没怎么想明白,只能留给内行点拨一二了……
发表于 2013-3-27 10:29 | 显示全部楼层
外行看的吃力呀
发表于 2013-3-27 19:35 | 显示全部楼层
本帖最后由 沐雨柠檬 于 2013-3-27 20:00 编辑

哎,下不了啊,体能怎么能多些?求附件123,请好心盟友发给我啊,正学呢,急需,10914439@qq.com,不胜感激!hht变换不是有虚假频率产生么?怎么解决?hht与小波变换相比,有何优缺点?
不用传了,自己解决了,需要干活赚钱才行。
发表于 2013-4-10 13:07 | 显示全部楼层
感觉看不懂啊,好难啊
发表于 2013-4-14 11:31 | 显示全部楼层
不知道小波尺度与频率到底什么关系?
发表于 2013-5-2 11:30 | 显示全部楼层
挺好的啊
发表于 2013-5-3 20:21 | 显示全部楼层
傻逼编的这程序全错的。。。无语
 楼主| 发表于 2013-5-3 22:28 | 显示全部楼层
thomasten 发表于 2013-5-3 20:21
傻逼编的这程序全错的。。。无语

兄弟指点一下,哪里有问题?
发表于 2013-5-30 16:55 | 显示全部楼层
体力值不够啊~
发表于 2013-9-6 21:10 | 显示全部楼层
发表于 2013-9-14 01:09 | 显示全部楼层
真心不错的分享,很佩服楼主这样的人
发表于 2013-9-27 14:48 | 显示全部楼层
找资料过来的,想下资料,却没体力,离开了
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表