;************************************************************** WAVE_SIGNIF ;+ ; NAME: WAVE_SIGNIF ; ; PURPOSE: Compute the significance levels for a wavelet transform. ; ; ; CALLING SEQUENCE: ; ; result = WAVE_SIGNIF(y,dt,scale,sigtest) ; ; ; INPUTS: ; ; Y = the time series, or, the VARIANCE of the time series. ; (If this is a single number, it is assumed to be the variance...) ; ; DT = amount of time between each Y value, i.e. the sampling time. ; ; SCALE = the vector of scale indices, from previous call to WAVELET. ; ; SIGTEST = 0, 1, or 2. If omitted, then assume 0. ; ; If 0 (the default), then just do a regular chi-square test, ; i.e. Eqn (18) from Torrence & Compo. ; If 1, then do a "time-average" test, i.e. Eqn (23). ; In this case, DOF should be set to NA, the number ; of local wavelet spectra that were averaged together. ; For the Global Wavelet Spectrum, this would be NA=N, ; where N is the number of points in your time series. ; If 2, then do a "scale-average" test, i.e. Eqns (25)-(28). ; In this case, DOF should be set to a ; two-element vector [S1,S2], which gives the scale ; range that was averaged together. ; e.g. if one scale-averaged scales between 2 and 8, ; then DOF=[2,8]. ; ; ; OUTPUTS: ; ; result = significance levels as a function of SCALE, ; or if /CONFIDENCE, then confidence intervals ; ; ; OPTIONAL KEYWORD INPUTS: ; ; MOTHER = A string giving the mother wavelet to use. ; Currently, 'Morlet','Paul','DOG' (derivative of Gaussian) ; are available. Default is 'Morlet'. ; ; PARAM = optional mother wavelet parameter. ; For 'Morlet' this is k0 (wavenumber), default is 6. ; For 'Paul' this is m (order), default is 4. ; For 'DOG' this is m (m-th derivative), default is 2. ; ; LAG1 = LAG 1 Autocorrelation, used for SIGNIF levels. Default is 0.0 ; ; SIGLVL = significance level to use. Default is 0.95 ; ; DOF = degrees-of-freedom for signif test. ; IF SIGTEST=0, then (automatically) DOF = 2 (or 1 for MOTHER='DOG') ; IF SIGTEST=1, then DOF = NA, the number of times averaged together. ; IF SIGTEST=2, then DOF = [S1,S2], the range of scales averaged. ; ; Note: IF SIGTEST=1, then DOF can be a vector (same length as SCALEs), ; in which case NA is assumed to vary with SCALE. ; This allows one to average different numbers of times ; together at different scales, or to take into account ; things like the Cone of Influence. ; See discussion following Eqn (23) in Torrence & Compo. ; ; GWS = global wavelet spectrum. If input then this is used ; as the theoretical background spectrum, ; rather than white or red noise. ; ; CONFIDENCE = if set, then return a Confidence INTERVAL. ; For SIGTEST=0,2 this will be two numbers, the lower & upper. ; For SIGTEST=1, this will return an array (J+1)x2, ; where J+1 is the number of scales. ; ; ; OPTIONAL KEYWORD OUTPUTS: ; ; PERIOD = the vector of "Fourier" periods (in time units) that corresponds ; to the SCALEs. ; ; FFT_THEOR = output theoretical red-noise spectrum as fn of PERIOD. ; ; ;---------------------------------------------------------------------------- ; ; EXAMPLE: ; ; IDL> wave = WAVELET(y,dt,PERIOD=period,SCALE=scale) ; IDL> signif = WAVE_SIGNIF(y,dt,scale) ; IDL> signif = REBIN(TRANSPOSE(signif),ntime,nscale) ; IDL> CONTOUR,ABS(wave)^2/signif,time,period, \$ ; LEVEL=1.0,C_ANNOT='95%' ; ; ;---------------------------------------------------------------------------- ; Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo, ; University of Colorado, Program in Atmospheric and Oceanic Sciences. ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties whatsoever. ; ; Notice: Please acknowledge the use of the above software in any publications: ; ``Wavelet software was provided by C. Torrence and G. Compo, ; and is available at URL: http://paos.colorado.edu/research/wavelets/''. ; ;---------------------------------------------------------------------------- ;- ;************************************************************* WAVE_SIGNIF FUNCTION wave_signif,y,dt,scale,sigtest, \$ ;*** required inputs MOTHER=mother,PARAM=param, \$ ;*** optional inputs LAG1=lag1,SIGLVL=siglvl,DOF=dof, \$ ;*** optional inputs GWS=gws,CONFIDENCE=confidence, \$ ;*** optional inputs FFT_THEOR=fft_theor,PERIOD=period, \$ ;*** optional outputs SAVG=Savg,SMID=Smid,CDELTA=CDelta,PSI0=psi0 ;*** optional outputs compile_opt idl2 ON_ERROR,2 IF (N_ELEMENTS(y) LT 1) THEN MESSAGE,'Time series Y must be input' IF (N_ELEMENTS(dt) LT 1) THEN MESSAGE,'DT must be input' IF (N_ELEMENTS(scale) LT 1) THEN MESSAGE,'Scales must be input' IF (N_PARAMS() LT 4) THEN sigtest = 0 ; the default variance = N_ELEMENTS(y) eq 1 ? y : (MOMENT(y))[1] ;....check keywords & optional inputs IF (N_ELEMENTS(mother) LT 1) THEN mother = 'MORLET' IF (N_ELEMENTS(param) LT 1) THEN param = -1 IF (N_ELEMENTS(siglvl) LT 1) THEN siglvl = 0.95 IF (N_ELEMENTS(lag1) LT 1) THEN lag1 = 0.0 confidence = KEYWORD_SET(confidence) lag1 = lag1[0] J = N_ELEMENTS(scale) - 1 s0 = MIN(scale) dj = ALOG(scale[1]/scale[0])/ALOG(2) CASE (STRUPCASE(mother)) OF 'MORLET': BEGIN IF (param EQ -1) THEN k0=6d ELSE k0=param fourier_factor = (4*!PI)/(k0 + SQRT(2+k0^2)) ; [Sec.3h] empir = [2.,-1,-1,-1] IF (k0 EQ 6) THEN empir[1:3]=[0.776,2.32,0.60] END 'PAUL': BEGIN ;****************** PAUL IF (param EQ -1) THEN m=4d ELSE m=param fourier_factor = 4*!PI/(2*m+1) empir = [2.,-1,-1,-1] IF (m EQ 4) THEN empir[1:3]=[1.132,1.17,1.5] END 'DOG': BEGIN ;******************* DOG IF (param EQ -1) THEN m=2 ELSE m=param fourier_factor = 2*!PI*SQRT(2./(2*m+1)) empir = [1.,-1,-1,-1] IF (m EQ 2) THEN empir[1:3] = [3.541,1.43,1.4] IF (m EQ 6) THEN empir[1:3] = [1.966,1.37,0.97] END ENDCASE period = scale*fourier_factor dofmin = empir[0] ; Degrees of freedom with no smoothing Cdelta = empir[1] ; reconstruction factor gamma = empir[2] ; time-decorrelation factor dj0 = empir[3] ; scale-decorrelation factor ;....significance levels [Sec.4] freq = dt/period ; normalized frequency fft_theor = (1-lag1^2)/(1-2*lag1*COS(freq*2*!PI)+lag1^2) ; [Eqn(16)] fft_theor = variance*fft_theor ; include time-series variance IF (N_ELEMENTS(gws) EQ (J+1)) THEN fft_theor = gws signif = fft_theor CASE (sigtest) OF 0: BEGIN ; no smoothing, DOF=dofmin dof = dofmin signif = fft_theor*CHISQR_CVF(1. - siglvl,dof)/dof ; [Eqn(18)] IF confidence THEN BEGIN sig = (1. - siglvl)/2. chisqr = dof/[CHISQR_CVF(sig,dof),CHISQR_CVF(1.-sig,dof)] signif = fft_theor # chisqr ENDIF END 1: BEGIN ; time-averaged, DOFs depend upon scale [Sec.5a] IF (N_ELEMENTS(dof) LT 1) THEN dof = dofmin IF (gamma EQ -1) THEN MESSAGE, \$ 'Gamma (decorrelation factor) not defined for '+mother+ \$ ' with param='+STRTRIM(param,2) IF (N_ELEMENTS(dof) EQ 1) THEN dof = FLTARR(J+1) + dof dof = dof > 1 dof = dofmin*SQRT( 1 + (dof*dt/gamma/scale)^2 ) ; [Eqn(23)] dof = dof > dofmin ; minimum DOF is dofmin IF (NOT confidence) THEN BEGIN FOR a1=0,J DO BEGIN chisqr = CHISQR_CVF(1. - siglvl,dof[a1])/dof[a1] signif[a1] = fft_theor[a1]*chisqr ; print,"a1=",a1,"dof=",dof[a1],"chisquare=", chisqr ENDFOR ENDIF ELSE BEGIN signif = FLTARR(J+1,2) sig = (1. - siglvl)/2. FOR a1=0,J DO BEGIN chisqr = dof(a1)/ \$ [CHISQR_CVF(sig,dof(a1)),CHISQR_CVF(1.-sig,dof(a1))] signif[a1,*] = fft_theor[a1]*chisqr ENDFOR ENDELSE END 2: BEGIN ; scale-averaged, DOFs depend upon scale range [Sec.5b] IF (N_ELEMENTS(dof) NE 2) THEN \$ MESSAGE,'DOF must be set to [S1,S2], the range of scale-averages' IF (Cdelta EQ -1) THEN MESSAGE, \$ 'Cdelta & dj0 not defined for '+mother+ \$ ' with param='+STRTRIM(param,2) s1 = dof[0] s2 = dof[1] avg = WHERE((scale GE s1) AND (scale LE s2),navg) IF (navg LT 1) THEN MESSAGE,'No valid scales between ' + \$ STRTRIM(s1,2) + ' and ' + STRTRIM(s2,2) s1 = MIN(scale[avg]) s2 = MAX(scale[avg]) Savg = 1./TOTAL(1./scale[avg]) ; [Eqn(25)] Smid = EXP((ALOG(s1)+ALOG(s2))/2.) ; power-of-two midpoint dof = (dofmin*navg*Savg/Smid)*SQRT(1 + (navg*dj/dj0)^2) ; [Eqn(28)] fft_theor = Savg*TOTAL(fft_theor[avg]/scale[avg]) ; [Eqn(27)] chisqr = CHISQR_CVF(1. - siglvl,dof)/dof IF confidence THEN BEGIN sig = (1. - siglvl)/2. chisqr = dof/[CHISQR_CVF(sig,dof),CHISQR_CVF(1.-sig,dof)] ENDIF signif = (dj*dt/Cdelta/Savg)*fft_theor*chisqr ; [Eqn(26)] END ENDCASE RETURN,signif END