小波变换(测试中)

Function

该方法仍在测试中，欢迎大家使用，如有问题请及时与我联系，联系方式详见右侧边栏about页面。

算法主体：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
import numpy as np
from numpy import fft
from numpy.fft import ifft
from scipy.special import gamma
from scipy.stats.distributions import chi2
import matplotlib.pyplot as plt
def wave_bases(mother,k,scale,param):
    n = k.shape[0]
    k_tmp = k.copy()
    k_tmp[k_tmp>0]=1
    k_tmp[k_tmp<=0]=0
    if (mother =='MORLET'):
        k0 = param
        expnt = (-1)*(scale*k - k0)**2/2*k_tmp
        norm = np.sqrt(scale*k[1])*(np.pi**(-0.25))*np.sqrt(n)
        daughter = norm*np.exp(expnt)
        daughter = daughter*k_tmp     
        fourier_factor = (4*np.pi)/(k0 + np.sqrt(2 + k0**2))
        coi = fourier_factor/np.sqrt(2)
        dofmin = 2
    elif (mother =='PAUL'):
        m = param
        expnt = (-1)*(scale*k)*k_tmp
        norm = np.sqrt(scale*k[1])*(2**m/np.sqrt(m*np.prod(np.arange(2,2*m))))*np.sqrt(n)
        daughter = norm*((scale*k)**m)*np.exp(expnt)
        daughter = daughter*k_tmp  
        fourier_factor = 4*np.pi/(2*m+1)
        coi = fourier_factor*np.sqrt(2)
        dofmin = 2
    elif (mother =='DOG'):
        m = param
        expnt = (-1)*(scale*k)**2/2.0
        norm = np.sqrt(scale*k[1]/gamma(m+0.5))*np.sqrt(n)
        daughter = (-1)*norm*((1j)**m)*((scale*k)**m)*np.exp(expnt)
        fourier_factor = 2*np.pi*np.sqrt(2/(2*m+1))
        coi = fourier_factor/np.sqrt(2)
        dofmin = 1
    return daughter,fourier_factor,coi,dofmin
 
def wavelet(Y,dt,pad,dj,s0,J1,mother,param):
    n1 = Y.shape[0]
    x = (Y - Y.mean())
    if (pad == 1):
        base2 = np.trunc(np.log(n1)/np.log(2) + 0.4999).astype('int32')
        x = np.hstack((x,np.zeros((2**(base2+1)-n1))))
    n = x.shape[0]
    k = np.arange(1,n//2+1,1)
    k = k*((2*np.pi)/(n*dt))
    k = np.hstack((0,k,((-1)*k[np.arange((n-1)//2-1,-1,-1)])))
    f = np.fft.fft(x)
    scale = s0*2**(np.arange(0,J1+1)*dj)
    period = scale
    wave = np.zeros((np.array(J1+1).astype('int32'),n))
    wave = wave + (1j)*wave
    for i in range(1,np.array(J1+2).astype('int32')):
        daughter,fourier_factor,coi,dofmin = wave_bases(mother,k,scale[i-1],param)
        wave[i-1,:] = ifft(f*daughter)
    period = fourier_factor*scale
    coi = coi*dt*np.hstack((1e-5,np.arange(1,(n1+1)//2,1),np.arange(1,n1//2,1)[::-1],1e-5))
    wave = wave[:,:n1]
    return wave,period,scale,coi

def wave_signif(Y,dt,scale1,sigtest,lag1,siglvl,dof,mother,param):
    n1 = np.array([Y]).shape
    J1 = scale1.shape[0] - 1
    scale[:J1+1] = scale1
    s0 = np.min(scale)
    dj = np.log(scale[1]/scale[0])/np.log(2)
    variance = Y
    if (mother =='MORLET'):
        k0 = param
        fourier_factor = (4*np.pi)/(k0 + np.sqrt(2 + k0**2))
        empir = np.array([2.,-1,-1,-1])
        if (k0==6):
            empir[1:]= np.array([0.776,2.32,0.60])
    elif (mother =='PAUL'):
        m = param
        fourier_factor = 4*np.pi/(2*m+1)
        empir = np.array([2.,-1,-1,-1])
        if (m==4):
            empir[1:]= np.array([1.132,1.17,1.5])
    elif (mother =='DOG'):
        m = param
        fourier_factor = 2*np.pi*np.sqrt(2/(2*m+1))
        empir = np.array([1,-1,-1,-1])
        if (m==2):
            empir[1:]= np.array([3.541,1.43,1.4])
        elif (m==6):
            empir[1:]= np.array([1.966,1.37,0.97])
    period = scale*fourier_factor
    dofmin = empir[0]
    Cdelta = empir[1]
    gamma_fac = empir[2]
    dj0 = empir[3]
    freq = dt / period
    fft_theor = (1-lag1**2) / (1-2*lag1*np.cos(freq*2*np.pi)+lag1**2)
    fft_theor = variance*fft_theor
    signif = fft_theor
    if (np.array(dof).sum() == -1):
        dof = dofmin
    if (sigtest == 0):
        dof = dofmin
        chisquare = chi2.ppf(siglvl,dof)/dof
        signif = fft_theor*chisquare
    elif (sigtest == 1):
        if (dof.shape[0] == 1):
            dof=np.zeros((J1+1))+dof
        dof[dof<1]=1
        dof = dofmin*np.sqrt(1 + (dof*dt/gamma_fac/ scale)**2 )
        dof[dof<dofmin]=dofmin
        for i in np.arange(1,J1+2):
            chisquare = chi2.ppf(siglvl,dof[i-1])/dof[i-1]
            signif[i-1] = fft_theor[i-1]*chisquare
    elif (sigtest == 2):
        s1 = dof[0]
        s2 = dof[1]
        avg = np.array(np.where((scale >= s1)&(scale < s2))).reshape(-1)
        navg = avg.shape[0]
        Savg = 1/np.sum(1/ scale[avg],axis=(0))
        Smid = np.exp((np.log(s1)+np.log(s2))/2.)
        dof = (dofmin*navg*Savg/Smid)*np.sqrt(1 + (navg*dj/dj0)**2)
        fft_theor = Savg*np.sum(fft_theor[avg]/ scale[avg])
        chisquare = chi2.ppf(siglvl,dof)/dof
        signif = (dj*dt/Cdelta/Savg)*fft_theor*chisquare    
    return signif,fft_theor

  1. wave_bases(mother,k,scale,param)函数（生成小波函数，仅在wavelet函数中调用，无需处理）

  2. wavelet(Y,dt,pad,dj,s0,J1,mother,param)函数 （小波变换本体函数）

    各参数含义如下：

    Y：输入数据，一维序列

    pad：是否开启边缘扩展。若为1，则将边缘值填补为0

    dt : 数据分辨率

    dj: 最小尺度，越小越精细，但计算绘图越慢

    s0：小波最小尺度，一般为2*dt

    j1: 刻度数(不太好理解)，范围从s0到s0×2^(j1×dj)；一般取j1 = (log2(N×dt/s0))/dj

    mother: 小波族，可选’MORLET’,’PAUL’,’DOG’三种小波

    param：小波族对应的常量，morlet为6；paul 为4； dog为2

    函数返回值：

    wave,period,scale,coi

  1. wave_signif(Y,dt,scale1,sigtest,lag1,siglvl,dof,mother,param) 函数 (显著性检验函数)

    各参数含义如下：

    scale1:wavelet函数返回的scale

    sigtest:显著性检验方法： 0为常规卡方检验；1为时间平均检验，对于局部小波，dof=np.nan,对于全局小波dof=N(N为时间序列的点数)；2为尺度平均检验，dof为一个[s1,s2]型，例如一个尺度的尺度平均在2-8之间，则自由度dof=[2,8]

    lag1:噪声自相关，用于显著性水平

    siglvl:显著性水平，0.95为95%显著性

    dof： 自由度

参数设置十分复杂，就目前而言还没找到能简单设置的方法，但是大部分参数是不需要修改的。

给出一个具体示例，首先是计算部分：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
np.set_printoptions(threshold=np.inf)
#读数据，与数据标准化，数据为126年季节平均的nino3指数（1871-1997年）
import pandas as pd
sst = np.array(pd.read_table('/data/home/zenggang/yxy/sst_nino3.dat',header =None)).reshape((-1))
variance = sst.std()*sst.std()
sst = ((sst - sst.mean())/sst.std())
#参数设置
n = sst.shape[0]
dt = 0.25 
time = np.arange(n)*dt + 1871.0 ;  #年份数据
pad = 1;     # 边缘效应
dj = 0.25    # 
s0 = 2*dt    # 最小尺度从逐6个月开始
j1 = 7/dj    # 
lag1 = 0.72  # 
mother = 'MORLET'
#计算小波
wave,period,scale,coi = wavelet(sst,dt,pad,dj,s0,j1,mother,6)
power = (np.abs(wave))**2
#显著性检验
signif,fft_theor = wave_signif(1,dt,scale,0,lag1,0.95,0,mother,6);
sig95 = (signif.T).reshape((-1,1))*(np.zeros((n))+1).reshape((1,-1))
sig95 = power/ sig95

global_ws = variance*(np.sum(power.T,axis = 0)/n)
dof = n - scale
global_signif,_ = wave_signif(variance,dt,scale,1,lag1,0.95,dof,mother,6)
avg = np.where((scale >= 2)&(scale < 8))
Cdelta = 0.776
scale_avg = (scale.T).reshape((-1,1))*(np.zeros((n))+1).reshape((1,-1))
scale_avg = power/ scale_avg
scale_avg = variance*dj*dt/Cdelta*np.sum(scale_avg[avg,:],axis = (0,1))
scaleavg_signif,_ = wave_signif(variance,dt,scale,2,lag1,0.95,np.array([2,7.9]),mother,6)

接下来是作图部分：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#原始数据序列
xlim = np.array([1870,2000])
fig = plt.figure(figsize=(15,15))
f_ax1 = fig.add_axes([0.1,0.75,0.65,0.2])
f_ax1.plot(time,sst)
f_ax1.set_xlabel('Time (year)')
f_ax1.set_ylabel('NINO3 SST (degC)')
f_ax1.set_title('a) NINO3 Sea Surface Temperature (seasonal)')
#小波
f_ax2 = fig.add_axes([0.1, 0.37 ,0.65 ,0.28])
levels = [0.0625,0.125,0.25,0.5,1,2,4,8,16]
Yticks = np.arange(2**np.trunc(np.log2(np.min(period))),2**np.trunc(np.log2(np.max(period)))+1,1)
f_ax2.contourf(time,np.log2(period),np.log2(power),np.log2(levels),extend='both')
f_ax2.set_xlabel('Time (year)')
f_ax2.set_ylabel('Period (years)')
f_ax2.set_title('b) NINO3 SST Wavelet Power Spectrum')
f_ax2.set_ylim(np.log2(np.max(period)),np.log2(np.min(period)))
f_ax2.contour(time,np.log2(period),sig95,levels=[-99,1],colors ='k')
f_ax2.plot(time,np.log2(coi),'k')
f_ax2.set_yticks(np.log2(Yticks)[[0,1,2,3,4,9,14,19,29,49]])
f_ax2.set_yticklabels(Yticks[[0,1,2,3,4,9,14,19,29,49]])
#全局功率谱
f_ax3 = fig.add_axes([0.77, 0.37, 0.2, 0.28])
f_ax3.plot(global_ws,np.log2(period))
f_ax3.plot(global_signif,np.log2(period),'--')
f_ax3.set_xlabel('Power (degC^2)')
f_ax3.set_title('c) Global Wavelet Spectrum')
f_ax3.set_ylim(np.log2([np.max(period),np.min(period)]))
f_ax3.set_yticks(np.log2(Yticks))
f_ax3.set_yticklabels('')
f_ax3.set_xlim(0,1.25*np.max(global_ws))
#局部功率谱
f_ax4 = fig.add_axes([0.1,0.07 ,0.65, 0.2])
f_ax4.plot(xlim,scaleavg_signif+np.array([0,0]),'--')
f_ax4.plot(time,scale_avg)
f_ax4.set_xlabel('Time (year)')
f_ax4.set_ylabel('Avg variance (degC^2)')
f_ax4.set_title('d) 2-8 yr Scale-average Time Series')

plt.show()

输出图形如下：

image-20200527155848450

测试数据下载：sst_nino3.dat

作为一个数学基础很差的人，理解这些算法实在是太难了，我在编译的过程中很多东西完全是从最后的图形结果倒推算法每一步的含义。程序几乎是参考C. Torrence(1998)的matlab程序完成的。如果大家在程序运行过程中发现Bug请及时向我反馈，关于参数设置问题，还是不要来为难我了，很多东西我只是有个模糊的概念，很难具体解释出来。不过在翻译过程中还是发现python相对于matlab在有些地方还是很方便的，比如傅里叶转换和逆转换，以及卡方检验的实现等等。

最后感谢大家一直以来的的支持和关注。