;******************************************************************* WAVELET
;+
; NAME:   WAVELET
;
; PURPOSE:   Compute the WAVELET transform of a 1D time series.
;       
;
; CALLING SEQUENCE:
;
;      wave = WAVELET(Y,DT)
;
;
; INPUTS:
;
;    Y = the time series of length N.
;
;    DT = amount of time between each Y value, i.e. the sampling time.
;
;
; OUTPUTS:
;
;    WAVE is the WAVELET transform of Y. This is a complex array
;    of dimensions (N,J+1). FLOAT(WAVE) gives the WAVELET amplitude,
;    ATAN(IMAGINARY(WAVE),FLOAT(WAVE)) gives the WAVELET phase.
;    The WAVELET power spectrum is ABS(WAVE)^2.
;
;
; OPTIONAL KEYWORD INPUTS:
;
;    S0 = the smallest scale of the wavelet.  Default is 2*DT.
;
;    DJ = the spacing between discrete scales. Default is 0.125.
;         A smaller # will give better scale resolution, but be slower to plot.
;
;    J = the # of scales minus one. Scales range from S0 up to S0*2^(J*DJ),
;        to give a total of (J+1) scales. Default is J = (LOG2(N DT/S0))/DJ.
;
;    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.
;
;    PAD = if set, then pad the time series with enough zeroes to get
;         N up to the next higher power of 2. This prevents wraparound
;         from the end of the time series to the beginning, and also
;         speeds up the FFT's used to do the wavelet transform.
;         This will not eliminate all edge effects (see COI below).
;
;    LAG1 = LAG 1 Autocorrelation, used for SIGNIF levels. Default is 0.0
;
;    SIGLVL = significance level to use. Default is 0.95
;
;    VERBOSE = if set, then print out info for each analyzed scale.
;
;    RECON = if set, then reconstruct the time series, and store in Y.
;            Note that this will destroy the original time series,
;            so be sure to input a dummy copy of Y.
;
;    FFT_THEOR = theoretical background spectrum as a function of
;                Fourier frequency. This will be smoothed by the
;                wavelet function and returned as a function of PERIOD.
;
;
; OPTIONAL KEYWORD OUTPUTS:
;
;    PERIOD = the vector of "Fourier" periods (in time units) that corresponds
;           to the SCALEs.
;
;    SCALE = the vector of scale indices, given by S0*2^(j*DJ), j=0...J
;            where J+1 is the total # of scales.
;
;    COI = if specified, then return the Cone-of-Influence, which is a vector
;        of N points that contains the maximum period of useful information
;        at that particular time.
;        Periods greater than this are subject to edge effects.
;        This can be used to plot COI lines on a contour plot by doing:
;            IDL>  CONTOUR,wavelet,time,period
;            IDL>  PLOTS,time,coi,NOCLIP=0
;
;    YPAD = returns the padded time series that was actually used in the
;         wavelet transform.
;
;    DAUGHTER = if initially set to 1, then return the daughter wavelets.
;         This is a complex array of the same size as WAVELET. At each scale
;         the daughter wavelet is located in the center of the array.
;
;    SIGNIF = output significance levels as a function of PERIOD
;
;    FFT_THEOR = output theoretical background spectrum (smoothed by the
;                wavelet function), as a function of PERIOD.
;
;
; [ Defunct INPUTS:
; [   OCT = the # of octaves to analyze over.           ]
; [         Largest scale will be S0*2^OCT.             ]
; [         Default is (LOG2(N) - 1).                   ]
; [   VOICE = # of voices in each octave. Default is 8. ]
; [          Higher # gives better scale resolution,    ]
; [          but is slower to plot.                     ]
; ]
;
;----------------------------------------------------------------------------
;
; EXAMPLE:
;
;    IDL> ntime = 256
;    IDL> y = RANDOMN(s,ntime)       ;*** create a random time series
;    IDL> dt = 0.25
;    IDL> time = FINDGEN(ntime)*dt   ;*** create the time index
;    IDL> 
;    IDL> wave = WAVELET(y,dt,PERIOD=period,COI=coi,/PAD,SIGNIF=signif)
;    IDL> nscale = N_ELEMENTS(period)
;    IDL> LOADCT,39
;    IDL> CONTOUR,ABS(wave)^2,time,period, $
;       XSTYLE=1,XTITLE='Time',YTITLE='Period',TITLE='Noise Wavelet', $
;       YRANGE=[MAX(period),MIN(period)], $   ;*** Large-->Small period
;       /YTYPE, $                             ;*** make y-axis logarithmic
;       NLEVELS=25,/FILL
;    IDL>
;    IDL> signif = REBIN(TRANSPOSE(signif),ntime,nscale)
;    IDL> CONTOUR,ABS(wave)^2/signif,time,period, $
;           /OVERPLOT,LEVEL=1.0,C_ANNOT='95%'
;    IDL> PLOTS,time,coi,NOCLIP=0   ;*** anything "below" this line is dubious
;
;
;----------------------------------------------------------------------------
; Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
;
; 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/''.
;
; Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
;            Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
;
; Please send a copy of such publications to either C. Torrence or G. Compo:
;  Dr. Christopher Torrence               Dr. Gilbert P. Compo
;  Research Systems, Inc.                 Climate Diagnostics Center
;  4990 Pearl East Circle                 325 Broadway R/CDC1
;  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
;  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
;----------------------------------------------------------------------------
;-

FUNCTION morlet, $ ;*********************************************** MORLET
	k0,scale,k,period,coi,dofmin,Cdelta,psi0

	IF (k0 EQ -1) THEN k0 = 6d
	n = N_ELEMENTS(k)
	expnt = -(scale*k - k0)^2/2d*(k GT 0.)
	dt = 2*!PI/(n*k(1))
	norm = SQRT(2*!PI*scale/dt)*(!PI^(-0.25))   ; total energy=N   [Eqn(7)]
	morlet = norm*EXP(expnt > (-100d))
	morlet = morlet*(expnt GT -100)  ; avoid underflow errors
	morlet = morlet*(k GT 0.)  ; Heaviside step function (Morlet is complex)
	fourier_factor = (4*!PI)/(k0 + SQRT(2+k0^2)) ; Scale-->Fourier [Sec.3h]
	period = scale*fourier_factor
	coi = fourier_factor/SQRT(2)   ; Cone-of-influence [Sec.3g]
	dofmin = 2   ; Degrees of freedom with no smoothing
	Cdelta = -1
	IF (k0 EQ 6) THEN Cdelta = 0.776 ; reconstruction factor
	psi0 = !PI^(-0.25)
;	PRINT,scale,n,SQRT(TOTAL(ABS(morlet)^2,/DOUBLE))
	RETURN,morlet
END

FUNCTION paul, $ ;************************************************** PAUL
	m,scale,k,period,coi,dofmin,Cdelta,psi0

	IF (m EQ -1) THEN m = 4d
	n = N_ELEMENTS(k)
	expnt = -(scale*k)*(k GT 0.)
	dt = 2d*!PI/(n*k(1))
	norm = SQRT(2*!PI*scale/dt)*(2^m/SQRT(m*FACTORIAL(2*m-1)))
	paul = norm*((scale*k)^m)*EXP(expnt > (-100d))*(expnt GT -100)
	paul = paul*(k GT 0.)
	fourier_factor = 4*!PI/(2*m+1)
	period = scale*fourier_factor
	coi = fourier_factor*SQRT(2)
	dofmin = 2   ; Degrees of freedom with no smoothing
	Cdelta = -1
	IF (m EQ 4) THEN Cdelta = 1.132 ; reconstruction factor
	psi0 = 2.^m*FACTORIAL(m)/SQRT(!PI*FACTORIAL(2*m))
;	PRINT,scale,n,norm,SQRT(TOTAL(paul^2,/DOUBLE))*SQRT(n)
	RETURN,paul
END

FUNCTION dog, $ ;*************************************************** DOG
	m,scale,k,period,coi,dofmin,Cdelta,psi0

	IF (m EQ -1) THEN m = 2
	n = N_ELEMENTS(k)
	expnt = -(scale*k)^2/2d
	dt = 2d*!PI/(n*k(1))
	norm = SQRT(2*!PI*scale/dt)*SQRT(1d/GAMMA(m+0.5))
	I = DCOMPLEX(0,1)
	gauss = -norm*(I^m)*(scale*k)^m*EXP(expnt > (-100d))*(expnt GT -100)
	fourier_factor = 2*!PI*SQRT(2./(2*m+1))
	period = scale*fourier_factor
	coi = fourier_factor/SQRT(2)
	dofmin = 1   ; Degrees of freedom with no smoothing
	Cdelta = -1
	psi0 = -1
	IF (m EQ 2) THEN BEGIN
		Cdelta = 3.541 ; reconstruction factor
		psi0 = 0.867325
	ENDIF
	IF (m EQ 6) THEN BEGIN
		Cdelta = 1.966 ; reconstruction factor
		psi0 = 0.88406
	ENDIF
;	PRINT,scale,n,norm,SQRT(TOTAL(ABS(gauss)^2,/DOUBLE))*SQRT(n)
	RETURN,gauss
END


;****************************************************************** WAVELET
FUNCTION wavelet,y1,dt, $   ;*** required inputs
	S0=s0,DJ=dj,J=j, $   ;*** optional inputs
	PAD=pad,MOTHER=mother,PARAM=param, $
	VERBOSE=verbose,NO_WAVE=no_wave,RECON=recon, $
	LAG1=lag1,SIGLVL=siglvl,DOF=dof,GLOBAL=global, $   ;*** optional inputs
	SCALE=scale,PERIOD=period,YPAD=ypad, $  ;*** optional outputs
	DAUGHTER=daughter,COI=coi, $
	SIGNIF=signif,FFT_THEOR=fft_theor, $
	OCT=oct,VOICE=voice   ;*** defunct inputs
	
	ON_ERROR,2
	r = CHECK_MATH(0,1)
	n = N_ELEMENTS(y1)
	n1 = n
	base2 = FIX(ALOG(n)/ALOG(2) + 0.4999)   ; power of 2 nearest to N

;....check keywords & optional inputs
	IF (N_ELEMENTS(s0) LT 1) THEN s0 = 2.0*dt
	IF (N_ELEMENTS(voice) EQ 1) THEN dj = 1./voice
	IF (N_ELEMENTS(dj) LT 1) THEN dj = 1./8
	IF (N_ELEMENTS(oct) EQ 1) THEN J = FLOAT(oct)/dj
	IF (N_ELEMENTS(J) LT 1) THEN J=FIX((ALOG(FLOAT(n)*dt/s0)/ALOG(2))/dj)  ;[Eqn(10)]
	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
	lag1 = lag1(0)
	verbose = KEYWORD_SET(verbose)
	do_daughter = KEYWORD_SET(daughter)
	do_wave = NOT KEYWORD_SET(no_wave)
	recon = KEYWORD_SET(recon)
	IF KEYWORD_SET(global) THEN MESSAGE, $
		'Please use WAVE_SIGNIF for global significance tests'
		
;....construct time series to analyze, pad if necessary
	ypad = y1 - TOTAL(y1)/n    ; remove mean
	IF KEYWORD_SET(pad) THEN BEGIN   ; pad with extra zeroes, up to power of 2
		ypad = [ypad,FLTARR(2L^(base2 + 1) - n)]
		n = N_ELEMENTS(ypad)
	ENDIF

;....construct SCALE array & empty PERIOD & WAVE arrays
	na = J + 1                  ; # of scales
	scale = DINDGEN(na)*dj      ; array of j-values
	scale = 2d0^(scale)*s0      ; array of scales  2^j   [Eqn(9)]
	period = FLTARR(na,/NOZERO) ; empty period array (filled in below)
	wave = COMPLEXARR(n,na,/NOZERO)  ; empty wavelet array
	IF (do_daughter) THEN daughter = wave   ; empty daughter array

;....construct wavenumber array used in transform [Eqn(5)]
	k = (DINDGEN(n/2) + 1)*(2*!PI)/(DOUBLE(n)*dt)
	k = [0d,k,-REVERSE(k(0:(n-1)/2 - 1))]

;....compute FFT of the (padded) time series
	yfft = FFT(ypad,-1,/DOUBLE)  ; [Eqn(3)]

	IF (verbose) THEN BEGIN  ;verbose
		PRINT
		PRINT,mother
		PRINT,'#points=',n1,'   s0=',s0,'   dj=',dj,'   J=',FIX(J)
		IF (n1 NE n) THEN PRINT,'(padded with ',n-n1,' zeroes)'
		PRINT,['j','scale','period','variance','mathflag'], $
			FORMAT='(/,A3,3A11,A10)'
	ENDIF  ;verbose
	IF (N_ELEMENTS(fft_theor) EQ n) THEN fft_theor_k = fft_theor ELSE $
		fft_theor_k = (1-lag1^2)/(1-2*lag1*COS(k*dt)+lag1^2)  ; [Eqn(16)]
	fft_theor = FLTARR(na)
	
;....loop thru each SCALE
	FOR a1=0,na-1 DO BEGIN  ;scale
		psi_fft=CALL_FUNCTION(mother, $
			param,scale(a1),k,period1,coi,dofmin,Cdelta,psi0)
		IF (do_wave) THEN $
			wave(*,a1) = FFT(yfft*psi_fft,1,/DOUBLE)  ;wavelet transform[Eqn(4)]
		period(a1) = period1   ; save period
		fft_theor(a1) = TOTAL((ABS(psi_fft)^2)*fft_theor_k)/n
		IF (do_daughter) THEN $
			daughter(*,a1) = FFT(psi_fft,1,/DOUBLE)   ; save daughter
		IF (verbose) THEN PRINT,a1,scale(a1),period(a1), $
				TOTAL(ABS(wave(*,a1))^2),CHECK_MATH(0), $
				FORMAT='(I3,3F11.3,I6)'
	ENDFOR  ;scale

	coi = coi*[FINDGEN((n1+1)/2),REVERSE(FINDGEN(n1/2))]*dt   ; COI [Sec.3g]
	
	IF (do_daughter) THEN $   ; shift so DAUGHTERs are in middle of array
		daughter = [daughter(n-n1/2:*,*),daughter(0:n1/2-1,*)]

;....significance levels [Sec.4]
	sdev = (MOMENT(y1))(1)
	fft_theor = sdev*fft_theor  ; include time-series variance
	dof = dofmin
	signif = fft_theor*CHISQR_CVF(1. - siglvl,dof)/dof   ; [Eqn(18)]

	IF (recon) THEN BEGIN  ; Reconstruction [Eqn(11)]
		IF (Cdelta EQ -1) THEN BEGIN
			y1 = -1
			MESSAGE,/INFO, $
				'Cdelta undefined, cannot reconstruct with this wavelet'
		ENDIF ELSE BEGIN
			y1=dj*SQRT(dt)/(Cdelta*psi0)*(FLOAT(wave) # (1./SQRT(scale)))
			y1 = y1[0:n1-1]
		ENDELSE
	ENDIF
	
	RETURN,wave(0:n1-1,*)    ; get rid of padding before returning

END


;************************************************************** 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



;****************************************************** WAVE_COHERENCY
;+
;   WAVE_COHERENCY
;
; PURPOSE:   Compute the wavelet coherency between two time series.
;
;
; CALLING SEQUENCE:
;
;     WAVE_COHERENCY, $
;         wave1,time1,scale1, $
;         wave2,time2,scale2, $
;         WAVE_COHER=wave_coher,WAVE_PHASE=wave_phase, $
;         TIME_OUT=time_out,SCALE_OUT=scale_out
;
;
; INPUTS:
;
;    WAVE1 = wavelet power spectrum for time series #1
;    TIME1 = a vector of times for time series #1
;    SCALE1 = a vector of scales for time series #1
;    WAVE2 = wavelet power spectrum for time series #2
;    TIME2 = a vector of times for time series #2
;    SCALE2 = a vector of scales for time series #2
;
;
; OPTIONAL KEYWORD INPUTS:
;
;    DT = amount of time between each Y value, i.e. the sampling time.
;         If not input, then calculated from TIME1(1)-TIME1(0)
;
;    DJ = the spacing between discrete scales.
;         If not input, then calculated from SCALE1
;
;   VERBOSE = if set, then print out the scales and system time
;
;   NOSMOOTH = if set, then just compute the GLOBAL_COHER, GLOBAL_PHASE,
;              and the unsmoothed CROSS_WAVELET and return
;
;
; OPTIONAL KEYWORD OUTPUTS:
;
;   WAVE_COHER = the wavelet coherency, as a function of
;       TIME_OUT and SCALE_OUT
;
;   TIME_OUT = the time vector, given by the overlap of TIME1 and TIME2
;
;   SCALE_OUT = the scale vector of scale indices, given by the overlap
;               of SCALE1 and SCALE2
;
;   COI_OUT = the vector of the cone-of-influence
;
;	GLOBAL_COHER = the global (or mean) coherence averaged over all times.
;
;   GLOBAL_PHASE = the global (or mean) phase averaged over all times
;
;	CROSS_WAVELET = the cross wavelet between the time series
;
;   POWER1 = the wavelet power spectrum; should be the same as WAVE1
;            if TIME1 and TIME2 are identical, otherwise it is only the
;            overlapping portion. If NOSMOOTH is set,
;            then this is unsmoothed, otherwise it is smoothed.
;
;   POWER2 = same as POWER1 but for time series #2
;
;
;----------------------------------------------------------------------------
; Copyright (C) 1998-2005, Christopher Torrence
; 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.
;
; Reference: Torrence, C. and P. J. Webster, 1999: Interdecadal changes in the
;            ENSO-monsoon system. <I>J. Climate</I>, 12, 2679-2690.
;
; Please send a copy of any publications to C. Torrence:
;  Dr. Christopher Torrence
;  Research Systems, Inc.
;  4990 Pearl East Circle
;  Boulder, CO 80301, USA
;  E-mail: chris[AT]rsinc[DOT]com
;----------------------------------------------------------------------------
;-



;****************************************************************** WAVELET
PRO wave_coherency, $
	wave1,time1,scale1,wave2,time2,scale2, $   ;*** required inputs
	COI1=coi1, $
	DT=dt,DJ=dj, $
	WAVE_COHER=wave_coher,WAVE_PHASE=wave_phase, $
	TIME_OUT=time_out,SCALE_OUT=scale_out,COI_OUT=coi_out, $
	GLOBAL_COHER=global_coher,GLOBAL_PHASE=global_phase, $
	CROSS_WAVELET=cross_wavelet,POWER1=power1,POWER2=power2, $
	NOSMOOTH=nosmooth, $
	VERBOSE=verbose

;	ON_ERROR,2

	verbose = KEYWORD_SET(verbose)

;*** find overlapping times
	time_start = MIN(time1) > MIN(time2)
	time_end = MAX(time1) < MAX(time2)
	time1_start = MIN(WHERE((time1 GE time_start)))
	time1_end = MAX(WHERE((time1 LE time_end)))
	time2_start = MIN(WHERE((time2 GE time_start)))
	time2_end = MAX(WHERE((time2 LE time_end)))

;*** find overlapping scales
	scale_start = MIN(scale1) > MIN(scale2)
	scale_end = MAX(scale1) < MAX(scale2)
	scale1_start = MIN(WHERE((scale1 GE scale_start)))
	scale1_end = MAX(WHERE((scale1 LE scale_end)))
	scale2_start = MIN(WHERE((scale2 GE scale_start)))
	scale2_end = MAX(WHERE((scale2 LE scale_end)))

;*** cross wavelet & individual wavelet power
	cross_wavelet = wave1(time1_start:time1_end,scale1_start:scale1_end)*$
		CONJ(wave2(time2_start:time2_end,scale2_start:scale2_end))
	power1 = ABS(wave1(time1_start:time1_end,scale1_start:scale1_end))^2
	power2 = ABS(wave2(time2_start:time2_end,scale2_start:scale2_end))^2

	IF (N_ELEMENTS(dt) LE 0) THEN dt = time1(1) - time1(0)
	ntime = time1_end - time1_start + 1
	nj = scale1_end - scale1_start + 1
	IF (N_ELEMENTS(dj) LE 0) THEN dj = ALOG(scale1(1)/scale1(0))/ALOG(2)
	scale = scale1(scale1_start:scale1_end)
	IF (verbose) THEN PRINT,dt,ntime,dj,nj
	time_out = time1(time1_start:time1_end)
	scale_out = scale1(scale1_start:scale1_end)
	IF (N_ELEMENTS(coi1) EQ N_ELEMENTS(time1)) THEN $
		coi_out = coi1(time1_start:time1_end)

; calculate global coherency before doing local smoothing
	global1 = TOTAL(power1,1)
	global2 = TOTAL(power2,1)
	global_cross = TOTAL(cross_wavelet,1)
	global_coher = ABS(global_cross)^2/(global1*global2)
	global_phase = 180./!PI*ATAN(IMAGINARY(global_cross),FLOAT(global_cross))

	IF KEYWORD_SET(nosmooth) THEN RETURN

	FOR j=0,nj-1 DO BEGIN ;*** time-smoothing
		st1 = SYSTIME(1)
		nt = LONG(4L*scale(j)/dt)/2L*4 + 1L
		time_wavelet = (FINDGEN(nt) - nt/2)*dt/scale(j)
		wave_function = EXP(-time_wavelet^2/2.)   ;*** Morlet
		wave_function = FLOAT(wave_function/TOTAL(wave_function)) ; normalize
		nz = nt/2
		zeros = COMPLEX(FLTARR(nz),FLTARR(nz))
		cross_wave_slice = [zeros,cross_wavelet(*,j),zeros]
		cross_wave_slice = CONVOL(cross_wave_slice,wave_function)
		cross_wavelet(*,j) = cross_wave_slice(nz:ntime+nz-1)
		zeros = FLOAT(zeros)
		power_slice = [zeros,power1(*,j),zeros]
		power_slice = CONVOL(power_slice,wave_function)
		power1(*,j) = power_slice(nz:ntime + nz - 1)
		power_slice = [zeros,power2(*,j),zeros]
		power_slice = CONVOL(power_slice,wave_function)
		power2(*,j) = power_slice(nz:ntime + nz - 1)
		IF (verbose) THEN PRINT,j,scale(j),SYSTIME(1)-st1;,FORMAT='(I4,$)'
	ENDFOR  ;*** time-smoothing

;*** normalize by scale
	scales = REBIN(TRANSPOSE(scale),ntime,nj)
	cross_wavelet = TEMPORARY(cross_wavelet)/scales
	power1 = TEMPORARY(power1)/scales
	power2 = TEMPORARY(power2)/scales

	nweights = FIX(0.6/dj/2 + 0.5)*2 - 1   ; closest (smaller) odd integer
	weights = REPLICATE(1.,nweights)
	weights = weights/TOTAL(weights) ; normalize
	FOR i=0,ntime-1 DO BEGIN ;*** scale-smoothing
		cross_wavelet(i,*) = CONVOL((cross_wavelet(i,*))(*),weights)
		power1(i,*) = CONVOL((power1(i,*))(*),weights)
		power2(i,*) = CONVOL((power2(i,*))(*),weights)
	ENDFOR ;*** scale-smoothing

	wave_phase = 180./!PI*ATAN(IMAGINARY(cross_wavelet),FLOAT(cross_wavelet))
	wave_coher = (ABS(cross_wavelet)^2)/(power1*power2 > 1E-9)
;	wave_phase = wave_phase + 360.*(wave_phase LT 0.)
END