wave_signif.pro (Documentation for /home/davin/idl/socware/)

;************************************************************** 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
	
	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
	IF (N_ELEMENTS(y) EQ 1) THEN variance=y ELSE variance=(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
			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