1 | ;************************************************************** WAVE_SIGNIF |
---|
2 | ;+ |
---|
3 | ; NAME: WAVE_SIGNIF |
---|
4 | ; |
---|
5 | ; PURPOSE: Compute the significance levels for a wavelet transform. |
---|
6 | ; |
---|
7 | ; |
---|
8 | ; CALLING SEQUENCE: |
---|
9 | ; |
---|
10 | ; result = WAVE_SIGNIF(y,dt,scale,sigtest) |
---|
11 | ; |
---|
12 | ; |
---|
13 | ; INPUTS: |
---|
14 | ; |
---|
15 | ; Y = the time series, or, the VARIANCE of the time series. |
---|
16 | ; (If this is a single number, it is assumed to be the variance...) |
---|
17 | ; |
---|
18 | ; DT = amount of time between each Y value, i.e. the sampling time. |
---|
19 | ; |
---|
20 | ; SCALE = the vector of scale indices, from previous call to WAVELET. |
---|
21 | ; |
---|
22 | ; SIGTEST = 0, 1, or 2. If omitted, then assume 0. |
---|
23 | ; |
---|
24 | ; If 0 (the default), then just do a regular chi-square test, |
---|
25 | ; i.e. Eqn (18) from Torrence & Compo. |
---|
26 | ; If 1, then do a "time-average" test, i.e. Eqn (23). |
---|
27 | ; In this case, DOF should be set to NA, the number |
---|
28 | ; of local wavelet spectra that were averaged together. |
---|
29 | ; For the Global Wavelet Spectrum, this would be NA=N, |
---|
30 | ; where N is the number of points in your time series. |
---|
31 | ; If 2, then do a "scale-average" test, i.e. Eqns (25)-(28). |
---|
32 | ; In this case, DOF should be set to a |
---|
33 | ; two-element vector [S1,S2], which gives the scale |
---|
34 | ; range that was averaged together. |
---|
35 | ; e.g. if one scale-averaged scales between 2 and 8, |
---|
36 | ; then DOF=[2,8]. |
---|
37 | ; |
---|
38 | ; |
---|
39 | ; OUTPUTS: |
---|
40 | ; |
---|
41 | ; result = significance levels as a function of SCALE, |
---|
42 | ; or if /CONFIDENCE, then confidence intervals |
---|
43 | ; |
---|
44 | ; |
---|
45 | ; OPTIONAL KEYWORD INPUTS: |
---|
46 | ; |
---|
47 | ; MOTHER = A string giving the mother wavelet to use. |
---|
48 | ; Currently, 'Morlet','Paul','DOG' (derivative of Gaussian) |
---|
49 | ; are available. Default is 'Morlet'. |
---|
50 | ; |
---|
51 | ; PARAM = optional mother wavelet parameter. |
---|
52 | ; For 'Morlet' this is k0 (wavenumber), default is 6. |
---|
53 | ; For 'Paul' this is m (order), default is 4. |
---|
54 | ; For 'DOG' this is m (m-th derivative), default is 2. |
---|
55 | ; |
---|
56 | ; LAG1 = LAG 1 Autocorrelation, used for SIGNIF levels. Default is 0.0 |
---|
57 | ; |
---|
58 | ; SIGLVL = significance level to use. Default is 0.95 |
---|
59 | ; |
---|
60 | ; DOF = degrees-of-freedom for signif test. |
---|
61 | ; IF SIGTEST=0, then (automatically) DOF = 2 (or 1 for MOTHER='DOG') |
---|
62 | ; IF SIGTEST=1, then DOF = NA, the number of times averaged together. |
---|
63 | ; IF SIGTEST=2, then DOF = [S1,S2], the range of scales averaged. |
---|
64 | ; |
---|
65 | ; Note: IF SIGTEST=1, then DOF can be a vector (same length as SCALEs), |
---|
66 | ; in which case NA is assumed to vary with SCALE. |
---|
67 | ; This allows one to average different numbers of times |
---|
68 | ; together at different scales, or to take into account |
---|
69 | ; things like the Cone of Influence. |
---|
70 | ; See discussion following Eqn (23) in Torrence & Compo. |
---|
71 | ; |
---|
72 | ; GWS = global wavelet spectrum. If input then this is used |
---|
73 | ; as the theoretical background spectrum, |
---|
74 | ; rather than white or red noise. |
---|
75 | ; |
---|
76 | ; CONFIDENCE = if set, then return a Confidence INTERVAL. |
---|
77 | ; For SIGTEST=0,2 this will be two numbers, the lower & upper. |
---|
78 | ; For SIGTEST=1, this will return an array (J+1)x2, |
---|
79 | ; where J+1 is the number of scales. |
---|
80 | ; |
---|
81 | ; |
---|
82 | ; OPTIONAL KEYWORD OUTPUTS: |
---|
83 | ; |
---|
84 | ; PERIOD = the vector of "Fourier" periods (in time units) that corresponds |
---|
85 | ; to the SCALEs. |
---|
86 | ; |
---|
87 | ; FFT_THEOR = output theoretical red-noise spectrum as fn of PERIOD. |
---|
88 | ; |
---|
89 | ; |
---|
90 | ;---------------------------------------------------------------------------- |
---|
91 | ; |
---|
92 | ; EXAMPLE: |
---|
93 | ; |
---|
94 | ; IDL> wave = WAVELET(y,dt,PERIOD=period,SCALE=scale) |
---|
95 | ; IDL> signif = WAVE_SIGNIF(y,dt,scale) |
---|
96 | ; IDL> signif = REBIN(TRANSPOSE(signif),ntime,nscale) |
---|
97 | ; IDL> CONTOUR,ABS(wave)^2/signif,time,period, $ |
---|
98 | ; LEVEL=1.0,C_ANNOT='95%' |
---|
99 | ; |
---|
100 | ; |
---|
101 | ;---------------------------------------------------------------------------- |
---|
102 | ; Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo, |
---|
103 | ; University of Colorado, Program in Atmospheric and Oceanic Sciences. |
---|
104 | ; This software may be used, copied, or redistributed as long as it is not |
---|
105 | ; sold and this copyright notice is reproduced on each copy made. This |
---|
106 | ; routine is provided as is without any express or implied warranties whatsoever. |
---|
107 | ; |
---|
108 | ; Notice: Please acknowledge the use of the above software in any publications: |
---|
109 | ; ``Wavelet software was provided by C. Torrence and G. Compo, |
---|
110 | ; and is available at URL: http://paos.colorado.edu/research/wavelets/''. |
---|
111 | ; |
---|
112 | ;---------------------------------------------------------------------------- |
---|
113 | ;- |
---|
114 | ;************************************************************* WAVE_SIGNIF |
---|
115 | FUNCTION wave_signif,y,dt,scale,sigtest, $ ;*** required inputs |
---|
116 | MOTHER=mother,PARAM=param, $ ;*** optional inputs |
---|
117 | LAG1=lag1,SIGLVL=siglvl,DOF=dof, $ ;*** optional inputs |
---|
118 | GWS=gws,CONFIDENCE=confidence, $ ;*** optional inputs |
---|
119 | FFT_THEOR=fft_theor,PERIOD=period, $ ;*** optional outputs |
---|
120 | SAVG=Savg,SMID=Smid,CDELTA=CDelta,PSI0=psi0 ;*** optional outputs |
---|
121 | |
---|
122 | ON_ERROR,2 |
---|
123 | IF (N_ELEMENTS(y) LT 1) THEN MESSAGE,'Time series Y must be input' |
---|
124 | IF (N_ELEMENTS(dt) LT 1) THEN MESSAGE,'DT must be input' |
---|
125 | IF (N_ELEMENTS(scale) LT 1) THEN MESSAGE,'Scales must be input' |
---|
126 | IF (N_PARAMS() LT 4) THEN sigtest = 0 ; the default |
---|
127 | IF (N_ELEMENTS(y) EQ 1) THEN variance=y ELSE variance=(MOMENT(y))(1) |
---|
128 | |
---|
129 | ;....check keywords & optional inputs |
---|
130 | IF (N_ELEMENTS(mother) LT 1) THEN mother = 'MORLET' |
---|
131 | IF (N_ELEMENTS(param) LT 1) THEN param = -1 |
---|
132 | IF (N_ELEMENTS(siglvl) LT 1) THEN siglvl = 0.95 |
---|
133 | IF (N_ELEMENTS(lag1) LT 1) THEN lag1 = 0.0 |
---|
134 | confidence = KEYWORD_SET(confidence) |
---|
135 | |
---|
136 | lag1 = lag1(0) |
---|
137 | |
---|
138 | J = N_ELEMENTS(scale) - 1 |
---|
139 | s0 = MIN(scale) |
---|
140 | dj = ALOG(scale(1)/scale(0))/ALOG(2) |
---|
141 | |
---|
142 | CASE (STRUPCASE(mother)) OF |
---|
143 | 'MORLET': BEGIN |
---|
144 | IF (param EQ -1) THEN k0=6d ELSE k0=param |
---|
145 | fourier_factor = (4*!PI)/(k0 + SQRT(2+k0^2)) ; [Sec.3h] |
---|
146 | empir = [2.,-1,-1,-1] |
---|
147 | IF (k0 EQ 6) THEN empir(1:3)=[0.776,2.32,0.60] |
---|
148 | END |
---|
149 | 'PAUL': BEGIN ;****************** PAUL |
---|
150 | IF (param EQ -1) THEN m=4d ELSE m=param |
---|
151 | fourier_factor = 4*!PI/(2*m+1) |
---|
152 | empir = [2.,-1,-1,-1] |
---|
153 | IF (m EQ 4) THEN empir(1:3)=[1.132,1.17,1.5] |
---|
154 | END |
---|
155 | 'DOG': BEGIN ;******************* DOG |
---|
156 | IF (param EQ -1) THEN m=2 ELSE m=param |
---|
157 | fourier_factor = 2*!PI*SQRT(2./(2*m+1)) |
---|
158 | empir = [1.,-1,-1,-1] |
---|
159 | IF (m EQ 2) THEN empir(1:3) = [3.541,1.43,1.4] |
---|
160 | IF (m EQ 6) THEN empir(1:3) = [1.966,1.37,0.97] |
---|
161 | END |
---|
162 | ENDCASE |
---|
163 | |
---|
164 | period = scale*fourier_factor |
---|
165 | dofmin = empir(0) ; Degrees of freedom with no smoothing |
---|
166 | Cdelta = empir(1) ; reconstruction factor |
---|
167 | gamma = empir(2) ; time-decorrelation factor |
---|
168 | dj0 = empir(3) ; scale-decorrelation factor |
---|
169 | |
---|
170 | ;....significance levels [Sec.4] |
---|
171 | freq = dt/period ; normalized frequency |
---|
172 | fft_theor = (1-lag1^2)/(1-2*lag1*COS(freq*2*!PI)+lag1^2) ; [Eqn(16)] |
---|
173 | fft_theor = variance*fft_theor ; include time-series variance |
---|
174 | IF (N_ELEMENTS(gws) EQ (J+1)) THEN fft_theor = gws |
---|
175 | signif = fft_theor |
---|
176 | |
---|
177 | CASE (sigtest) OF |
---|
178 | |
---|
179 | 0: BEGIN ; no smoothing, DOF=dofmin |
---|
180 | dof = dofmin |
---|
181 | signif = fft_theor*CHISQR_CVF(1. - siglvl,dof)/dof ; [Eqn(18)] |
---|
182 | IF confidence THEN BEGIN |
---|
183 | sig = (1. - siglvl)/2. |
---|
184 | chisqr = dof/[CHISQR_CVF(sig,dof),CHISQR_CVF(1.-sig,dof)] |
---|
185 | signif = fft_theor # chisqr |
---|
186 | ENDIF |
---|
187 | END |
---|
188 | |
---|
189 | 1: BEGIN ; time-averaged, DOFs depend upon scale [Sec.5a] |
---|
190 | IF (N_ELEMENTS(dof) LT 1) THEN dof = dofmin |
---|
191 | IF (gamma EQ -1) THEN MESSAGE, $ |
---|
192 | 'Gamma (decorrelation factor) not defined for '+mother+ $ |
---|
193 | ' with param='+STRTRIM(param,2) |
---|
194 | IF (N_ELEMENTS(dof) EQ 1) THEN dof = FLTARR(J+1) + dof |
---|
195 | dof = dof > 1 |
---|
196 | dof = dofmin*SQRT( 1 + (dof*dt/gamma/scale)^2 ) ; [Eqn(23)] |
---|
197 | dof = dof > dofmin ; minimum DOF is dofmin |
---|
198 | IF (NOT confidence) THEN BEGIN |
---|
199 | FOR a1=0,J DO BEGIN |
---|
200 | chisqr = CHISQR_CVF(1. - siglvl,dof(a1))/dof(a1) |
---|
201 | signif(a1) = fft_theor(a1)*chisqr |
---|
202 | ENDFOR |
---|
203 | ENDIF ELSE BEGIN |
---|
204 | signif = FLTARR(J+1,2) |
---|
205 | sig = (1. - siglvl)/2. |
---|
206 | FOR a1=0,J DO BEGIN |
---|
207 | chisqr = dof(a1)/ $ |
---|
208 | [CHISQR_CVF(sig,dof(a1)),CHISQR_CVF(1.-sig,dof(a1))] |
---|
209 | signif(a1,*) = fft_theor(a1)*chisqr |
---|
210 | ENDFOR |
---|
211 | ENDELSE |
---|
212 | END |
---|
213 | |
---|
214 | 2: BEGIN ; scale-averaged, DOFs depend upon scale range [Sec.5b] |
---|
215 | IF (N_ELEMENTS(dof) NE 2) THEN $ |
---|
216 | MESSAGE,'DOF must be set to [S1,S2], the range of scale-averages' |
---|
217 | IF (Cdelta EQ -1) THEN MESSAGE, $ |
---|
218 | 'Cdelta & dj0 not defined for '+mother+ $ |
---|
219 | ' with param='+STRTRIM(param,2) |
---|
220 | s1 = dof(0) |
---|
221 | s2 = dof(1) |
---|
222 | avg = WHERE((scale GE s1) AND (scale LE s2),navg) |
---|
223 | IF (navg LT 1) THEN MESSAGE,'No valid scales between ' + $ |
---|
224 | STRTRIM(s1,2) + ' and ' + STRTRIM(s2,2) |
---|
225 | s1 = MIN(scale(avg)) |
---|
226 | s2 = MAX(scale(avg)) |
---|
227 | Savg = 1./TOTAL(1./scale(avg)) ; [Eqn(25)] |
---|
228 | Smid = EXP((ALOG(s1)+ALOG(s2))/2.) ; power-of-two midpoint |
---|
229 | dof = (dofmin*navg*Savg/Smid)*SQRT(1 + (navg*dj/dj0)^2) ; [Eqn(28)] |
---|
230 | fft_theor = Savg*TOTAL(fft_theor(avg)/scale(avg)) ; [Eqn(27)] |
---|
231 | chisqr = CHISQR_CVF(1. - siglvl,dof)/dof |
---|
232 | IF confidence THEN BEGIN |
---|
233 | sig = (1. - siglvl)/2. |
---|
234 | chisqr = dof/[CHISQR_CVF(sig,dof),CHISQR_CVF(1.-sig,dof)] |
---|
235 | ENDIF |
---|
236 | signif = (dj*dt/Cdelta/Savg)*fft_theor*chisqr ; [Eqn(26)] |
---|
237 | END |
---|
238 | |
---|
239 | ENDCASE |
---|
240 | |
---|
241 | RETURN,signif |
---|
242 | |
---|
243 | END |
---|