Morlet小波变换脚本示例参考

qxsw2016

根据https://wenku.baidu.com/view/8d32e32058fb770bf78a5585.html提供的数据计算某流域径流量的小波分析，采用Morlet小波。w0取为5，脚本如下：

1. def cpsai(t): #Morlet 小波函数
   
2.     w0 = 5
   
3.     a = cos(w0*t)*exp(-t*t/2)
   
4.     b = sin(w0*t)*exp(-t*t/2)
   
5.     a = a*(pi**(-0.25))
   
6.     b = b*(pi**(-0.25))
   
7.     return complex(a,b)
   
8.    
   
9. n = 39 #样本量
   
10. a = 32 #时间尺度因子
    
11. m = 1
    
12. datafn = 'F:\DataSet\sw.txt'
    
13. swdata = asciiread(datafn, shape=(n,m))
    
14. swdata = swdata - sum(swdata)/39 #距平化
    
15. an = 32
    
16. x =  arange(0,n,0.5) #时间平移因子离散
    
17. y =  arange(0.5,a,0.5) #时间尺度因子离散
    
18. an = len(y)
    
19. bn = len(x)
    
20. wf = zeros([an,bn])
    
21. deltat = 1
    
22. for b in arange(0,bn):
    
23.     for a in arange(0,an):
    
24.         s = 0
    
25.         for k in arange(0,n):
    
26.             s = s + swdata[k]*cpsai((k*deltat-x)/y[a])
    
27.         wf[a,b] = (deltat*s/sqrt(abs(y[a]))).real #小波系数
    
28.         #r = (deltat*s/sqrt(abs(y[a]))).real
    
29.         #i = (deltat*s/sqrt(abs(y[a]))).imag
    
30.         #wf[a,b] = sqrt(r*r+i*i) #小波变换的模
    
31.     print a,b
    
32. xtick = arange(1966,1966+bn)
    
33. ytick = arange(0,an)
    
34. layer = contour(xtick,ytick,wf,15,smooth = True)
    
35. clabel(layer)
    
36. xlabel('b(Year)')
    
37. ylabel('a')
    
38. title(u'小波系数实部等值线图',fontname='宋体', fontsize=16)

def cpsai(t): #Morlet 小波函数<br />     w0 = 5<br />     a = cos(w0*t)*exp(-t*t/2)<br />     b = sin(w0*t)*exp(-t*t/2)<br />     a = a*(pi**(-0.25))<br />     b = b*(pi**(-0.25))<br />     return complex(a,b)<br />     <br /> n = 39 #样本量<br /> a = 32 #时间尺度因子<br /> m = 1<br /> datafn = 'F:\DataSet\sw.txt'<br /> swdata = asciiread(datafn, shape=(n,m))<br /> swdata = swdata - sum(swdata)/39 #距平化<br /> an = 32<br /> x =  arange(0,n,0.5) #时间平移因子离散<br /> y =  arange(0.5,a,0.5) #时间尺度因子离散<br /> an = len(y)<br /> bn = len(x)<br /> wf = zeros([an,bn])<br /> deltat = 1<br /> for b in arange(0,bn):<br />     for a in arange(0,an):<br />         s = 0<br />         for k in arange(0,n):<br />             s = s + swdata[k]*cpsai((k*deltat-x<strong>)/y[a])<br />         wf[a,b] = (deltat*s/sqrt(abs(y[a]))).real #小波系数<br />         #r = (deltat*s/sqrt(abs(y[a]))).real<br />         #i = (deltat*s/sqrt(abs(y[a]))).imag<br />         #wf[a,b] = sqrt(r*r+i*i) #小波变换的模<br />     print a,b<br /> xtick = arange(1966,1966+bn)<br /> ytick = arange(0,an)<br /> layer = contour(xtick,ytick,wf,15,smooth = True)<br /> clabel(layer)<br /> xlabel('b(Year)')<br /> ylabel('a')<br /> title(u'小波系数实部等值线图',fontname='宋体', fontsize=16)

qxsw2016

计算的实部等值线图

︶ㄣ安定■丶

很好的例子，谢谢分享啊。

MeteoInfo

sw.txt文件能提供一下吗？多谢！

qxsw2016

matlab有个数据的处理，对称性延伸，这个步骤进行了省略。另外，原值分析、距平值分析会不同。

MeteoInfo

qxsw2016 发表于 2018-5-23 12:39
matlab有个数据的处理，对称性延伸，这个步骤进行了省略。另外，原值分析、距平值分析会不同。

数据的年份扩展了那么多如何解释？

qxsw2016

数据 序列是1966-2004年共39年的实测数据，语句x =  arange(0,n,0.5) #时间平移因子离散，取间隔0.5，横轴范围[0,n-1],语句xtick = arange(1966,1966+bn)改成xtick = arange(1966,1966+n,0.5)

qxsw2016

纵轴刻度也有问题，改成ytick = arange(0.5,a+0.5,0.5)