wave_idl/ 0042775 0007030 0000404 00000000000 10256306452 0013016 5 ustar 00laurie web 0000166 0000006 wave_idl/wavelet.pro 0100664 0007030 0000404 00000030363 10144200076 0015176 0 ustar 00laurie web 0000166 0000006 ;******************************************************************* 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. Bull. Amer. Meteor. Soc., 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
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.wave_idl/wave_signif.pro 0100664 0007030 0000404 00000021430 06662644735 0016053 0 ustar 00laurie web 0000166 0000006 ;************************************************************** 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
ackground 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 SIGTEwave_idl/wavetest.pro 0100664 0007030 0000404 00000013527 07020021330 0015365 0 ustar 00laurie web 0000166 0000006 ;************************************************************** WAVETEST
;+
; NAME: WAVETEST
;
; PURPOSE: Example IDL program for WAVELET, using NINO3 SST dataset
;
; EXECUTION:
;
; IDL> .run wavetest
;
;
; See "http://paos.colorado.edu/research/wavelets/"
; Written January 1998 by C. Torrence
;
;-
;**************************************************************
n = 504
sst = FLTARR(n)
OPENR,1,'sst_nino3.dat' ; input SST time series
READF,1,sst
CLOSE,1
;------------------------------------------------------ Computation
; normalize by standard deviation (not necessary, but makes it easier
; to compare with plot on Interactive Wavelet page, at
; "http://paos.colorado.edu/research/wavelets/plot/"
sst = (sst - TOTAL(sst)/n)
dt = 0.25
time = FINDGEN(n)*dt + 1871.0 ; construct time array
xrange = [1870,2000] ; plotting range
pad = 1
s0 = dt ; this says start at a scale of 3 months
dj = 0.25 ; this will do 4 sub-octaves per octave
j1 = 9./dj ; this says do 9 powers-of-two with dj sub-octaves each
mother = 'Morlet'
recon_sst = sst ; save an extra copy, so we don't erase original sst
; estimate lag-1 autocorrelation, for red-noise significance tests
; Note that we actually use the global wavelet spectrum (GWS)
; for the significance tests, but if you wanted to use red noise,
; here's how you could calculate it...
lag1 = (A_CORRELATE(sst,1) + SQRT(A_CORRELATE(sst,2)))/2.
; Wavelet transform:
wave = WAVELET(recon_sst,dt,PERIOD=period,SCALE=scale,S0=s0, $
PAD=pad,COI=coi,DJ=dj,J=j1,MOTHER=mother,/RECON)
power = (ABS(wave))^2 ; compute wavelet power spectrum
global_ws = TOTAL(power,1)/n ; global wavelet spectrum (GWS)
J = N_ELEMENTS(scale) - 1
; Significance levels, assuming the GWS as background spectrum:
signif = WAVE_SIGNIF(sst,dt,scale,0, $
GWS=global_ws,SIGLVL=0.90,MOTHER=mother)
signif = REBIN(TRANSPOSE(signif),n,J+1) ; expand signif --> (J+1)x(N) array
signif = power/signif ; where ratio > 1, power is significant
; GWS significance levels:
dof = n - scale ; the -scale corrects for padding at edges
global_signif = WAVE_SIGNIF(sst,dt,scale,1, $
LAG1=0.0,DOF=dof,MOTHER=mother,CDELTA=Cdelta,PSI0=psi0)
; check total variance (Parseval's theorem) [Eqn(14)]
scale_avg = REBIN(TRANSPOSE(scale),n,J+1) ; expand scale-->(J+1)x(N) array
power_norm = power/scale_avg
variance = (MOMENT(sst))(1)
recon_variance = dj*dt/(Cdelta*n)*TOTAL(power_norm) ; [Eqn(14)]
IF (N_ELEMENTS(recon_sst) GT 1) THEN BEGIN
recon_variance = (MOMENT(recon_sst))(1)
; RMS of Reconstruction [Eqn(11)]
rms_error = SQRT(TOTAL((sst - recon_sst)^2)/n)
PRINT
PRINT,' ******** RECONSTRUCTION ********'
PRINT,'original variance =',variance,' degC^2'
PRINT,'reconstructed var =',FLOAT(recon_variance),' degC^2'
PRINT,'Ratio = ',recon_variance/variance
PRINT,'root-mean-square error of reconstructed sst = ',rms_error,' degC'
PRINT
IF (mother EQ 'DOG') THEN BEGIN
PRINT,'Note: for better reconstruction with the DOG, you need'
PRINT,' to use a very small s0.'
ENDIF
PRINT
ENDIF
; Scale-average between El Nino periods of 2--8 years
avg = WHERE((scale GE 2) AND (scale LT 8))
scale_avg = dj*dt/Cdelta*TOTAL(power_norm(*,avg),2) ; [Eqn(24)]
scaleavg_signif = WAVE_SIGNIF(sst,dt,scale,2, $
GWS=global_ws,SIGLVL=0.90,DOF=[2,7.9],MOTHER=mother)
;------------------------------------------------------ Plotting
printfile = 0
!P.FONT = -1
!P.CHARSIZE = 1
IF (printfile) THEN BEGIN
SET_PLOT,'ps'
DEVICE,/PORT,/INCH,XSIZE=6.5,XOFF=1,YSIZE=6,YOFF=3,/COLOR,BITS=8
!P.FONT = 0
!P.CHARSIZE = 0.75
ENDIF ELSE WINDOW,0,XSIZE=600,YSIZE=600
!P.MULTI = 0
!X.STYLE = 1
!Y.STYLE = 1
LOADCT,39
;--- Plot time series
pos1 = [0.1,0.75,0.7,0.95]
PLOT,time,sst,XRANGE=xrange, $
XTITLE='Time (year)',YTITLE='NINO3 SST (!Uo!NC)', $
TITLE='a) NINO3 Sea Surface Temperature (seasonal)', $
POSITION=pos1
IF (N_ELEMENTS(recon_sst) GT 1) THEN OPLOT,time,recon_sst,COLOR=144
XYOUTS,0.85,0.9,/NORMAL,ALIGN=0.5, $
'!5WAVELET ANALYSIS!X'+$
'!C!CC. Torrence & G.P. Compo'+$
'!C!Chttp://paos.colorado.edu/!Cresearch/wavelets/'
;--- Contour plot wavelet power spectrum
yrange = [64,0.5] ; years
levels = [0.5,1,2,4]
colors = [64,128,208,254]
period2 = FIX(ALOG(period)/ALOG(2)) ; integer powers of 2 in period
ytickv = 2.^(period2(UNIQ(period2))) ; unique powers of 2
pos2 = [pos1(0),0.35,pos1(2),0.65]
CONTOUR,power,time,period,/NOERASE,POSITION=pos2, $
XRANGE=xrange,YRANGE=yrange,/YTYPE, $
YTICKS=N_ELEMENTS(ytickv)-1,YTICKV=ytickv, $
LEVELS=levels,C_COLORS=colors,/FILL, $
XTITLE='Time (year)',YTITLE='Period (years)', $
TITLE='b) Wavelet Power Spectrum (contours at 0.5,1,2,4!Uo!NC!U2!N)'
; significance contour, levels at -99 (fake) and 1 (significant)
CONTOUR,signif,time,period,/OVERPLOT,LEVEL=1,THICK=2, $
C_LABEL=1,C_ANNOT='90%',C_CHARSIZE=1
; cone-of-influence, anything "below" is dubious
x = [time(0),time,MAX(time)]
y = [MAX(period),coi,MAX(period)]
color = 4
POLYFILL,x,y,ORIEN=+45,SPACING=0.5,COLOR=color,NOCLIP=0,THICK=1
POLYFILL,x,y,ORIEN=-45,SPACING=0.5,COLOR=color,NOCLIP=0,THICK=1
PLOTS,time,coi,COLOR=color,NOCLIP=0,THICK=1
;--- Plot global wavelet spectrum
pos3 = [0.74,pos2(1),0.95,pos2(3)]
blank = REPLICATE(' ',29)
PLOT,global_ws,period,/NOERASE,POSITION=pos3, $
THICK=2,XSTYLE=10,YSTYLE=9, $
YRANGE=yrange,/YTYPE,YTICKLEN=-0.02, $
XTICKS=2,XMINOR=2, $
YTICKS=N_ELEMENTS(ytickv)-1,YTICKV=ytickv,YTICKNAME=blank, $
XTITLE='Power (!Uo!NC!U2!N)',TITLE='c) Global'
OPLOT,global_signif,period,LINES=1
XYOUTS,1.7,60,'95%'
;--- Plot 2--8 yr scale-average time series
pos4 = [pos1(0),0.05,pos1(2),0.25]
PLOT,time,scale_avg,/NOERASE,POSITION=pos4, $
XRANGE=xrange,YRANGE=[0,MAX(scale_avg)*1.25],THICK=2, $
XTITLE='Time (year)',YTITLE='Avg variance (!Uo!NC!U2!N)', $
TITLE='d) 2-8 yr Scale-average Time Series'
OPLOT,xrange,scaleavg_signif+[0,0],LINES=1
IF (printfile) THEN DEVICE,/CLOSE
END
wave_idl/sst_nino3.dat 0100664 0007030 0000404 00000005720 06457461064 0015437 0 ustar 00laurie web 0000166 0000006 -0.15
-0.30
-0.14
-0.41
-0.46
-0.66
-0.50
-0.80
-0.95
-0.72
-0.31
-0.71
-1.04
-0.77
-0.86
-0.84
-0.41
-0.49
-0.48
-0.72
-1.21
-0.80
0.16
0.46
0.40
1.00
2.17
2.50
2.34
0.80
0.14
-0.06
-0.34
-0.71
-0.34
-0.73
-0.48
-0.11
0.22
0.51
0.51
0.25
-0.10
-0.33
-0.42
-0.23
-0.53
-0.44
-0.30
0.15
0.09
0.19
-0.06
0.25
0.30
0.81
0.26
0.10
0.34
1.01
-0.31
-0.90
-0.73
-0.92
-0.73
-0.31
-0.03
0.12
0.37
0.82
1.22
1.83
1.60
0.34
-0.72
-0.87
-0.85
-0.40
-0.39
-0.65
0.07
0.67
0.39
0.03
-0.17
-0.76
-0.87
-1.36
-1.10
-0.99
-0.78
-0.93
-0.87
-0.44
-0.34
-0.50
-0.39
-0.04
0.42
0.62
0.17
0.23
1.03
1.54
1.09
0.01
0.12
-0.27
-0.47
-0.41
-0.37
-0.36
-0.39
0.43
1.05
1.58
1.25
0.86
0.60
0.21
0.19
-0.23
-0.29
0.18
0.12
0.71
1.42
1.59
0.93
-0.25
-0.66
-0.95
-0.47
0.06
0.70
0.81
0.78
1.43
1.22
1.05
0.44
-0.35
-0.67
-0.84
-0.66
-0.45
-0.12
-0.20
-0.16
-0.47
-0.52
-0.79
-0.80
-0.62
-0.86
-1.29
-1.04
-1.05
-0.75
-0.81
-0.90
-0.25
0.62
1.22
0.96
0.21
-0.11
-0.25
-0.24
-0.43
0.23
0.67
0.78
0.41
0.98
1.28
1.45
1.02
0.03
-0.59
-1.34
-0.99
-1.49
-1.74
-1.33
-0.55
-0.51
-0.36
-0.99
0.32
1.04
1.41
0.99
0.66
0.50
0.22
0.71
-0.16
0.38
0.00
-1.11
-1.04
0.05
-0.64
-0.34
-0.50
-1.85
-0.94
-0.78
0.29
0.27
0.69
-0.06
-0.83
-0.80
-1.02
-0.96
-0.09
0.62
0.87
1.03
0.70
-0.10
-0.31
0.04
-0.46
0.04
0.24
-0.08
-0.28
0.06
0.05
-0.31
0.11
0.27
0.26
0.04
0.12
1.11
1.53
1.23
0.17
-0.18
-0.56
0.05
0.41
0.22
0.04
-0.19
-0.46
-0.65
-1.06
-0.54
0.14
0.25
-0.21
-0.73
-0.43
0.48
0.26
0.05
0.11
-0.27
-0.08
-0.10
0.29
-0.15
-0.28
-0.55
-0.44
-1.40
-0.55
-0.69
0.58
0.37
0.42
1.83
1.23
0.65
0.41
1.03
0.64
-0.07
0.98
0.36
-0.30
-1.33
-1.39
-0.94
0.34
-0.00
-0.15
0.06
0.39
0.36
-0.49
-0.53
0.35
0.07
-0.24
0.20
-0.22
-0.68
-0.44
0.02
-0.22
-0.30
-0.59
0.10
-0.02
-0.27
-0.60
-0.48
-0.37
-0.53
-1.35
-1.22
-0.99
-0.34
-0.79
-0.24
0.02
0.69
0.78
0.17
-0.17
-0.29
-0.27
0.31
0.44
0.38
0.24
-0.13
-0.89
-0.76
-0.71
-0.37
-0.59
-0.63
-1.47
-0.40
-0.18
-0.37
-0.43
-0.06
0.61
1.33
1.19
1.13
0.31
0.14
0.03
0.21
0.15
-0.22
-0.02
0.03
-0.17
0.12
-0.35
-0.06
0.38
-0.45
-0.32
-0.33
-0.49
-0.14
-0.56
-0.18
0.46
1.09
1.04
0.23
-0.99
-0.59
-0.92
-0.28
0.52
1.31
1.45
0.61
-0.11
-0.18
-0.39
-0.39
-0.36
-0.50
-0.81
-1.10
-0.29
0.57
0.68
0.78
0.78
0.63
0.98
0.49
-0.42
-1.34
-1.20
-1.18
-0.65
-0.42
-0.97
-0.28
0.77
1.77
2.22
1.05
-0.67
-0.99
-1.52
-1.17
-0.22
-0.04
-0.45
-0.46
-0.75
-0.70
-1.38
-1.15
-0.01
0.97
1.10
0.68
-0.02
-0.04
0.47
0.30
-0.55
-0.51
-0.09
-0.01
0.34
0.61
0.58
0.33
0.38
0.10
0.18
-0.30
-0.06
-0.28
0.12
0.58
0.89
0.93
2.39
2.44
1.92
0.64
-0.24
0.27
-0.13
-0.16
-0.54
-0.13
-0.37
-0.78
-0.22
0.03
0.25
0.31
1.03
1.10
1.05
1.11
1.28
0.57
-0.55
-1.16
-0.99
-0.38
0.01
-0.29
0.09
0.46
0.57
0.24
0.39
0.49
0.86
0.51
0.95
1.25
1.33
-0.00
0.34
0.66
1.11
0.34
0.48
0.56
0.39
-0.17
1.04
0.77
0.12
-0.35
-0.22
0.08
-0.08
-0.18
-0.06
wave_idl/monsoon.txt 0100664 0007030 0000404 00000007140 07433251306 0015244 0 ustar 00laurie web 0000166 0000006 $submit
Submit
$dataset
nino3sst
$title
All-India Monsoon Rainfall
$source
Parthasarathy et al. 1991: J. Climate, 4, 927-938
$dataset1
Rainfall
$units
mm
t0
1871.0
dt
0.25
$t_units
years
data
1.97
21.17
-15.97
-14.79
-5.46
5.00
11.89
1.65
-2.23
-20.67
-14.81
-13.02
-1.23
23.70
19.33
-1.45
-1.56
8.10
16.13
-16.75
-6.66
-19.90
-11.44
-20.62
12.57
1.03
-75.71
14.28
-2.90
-3.00
52.33
0.31
-6.00
14.67
6.36
-1.99
-1.13
5.40
-18.67
10.48
-0.46
-1.57
3.99
-10.02
-3.50
16.47
0.43
6.48
0.54
13.63
-14.57
2.18
-5.43
-8.03
28.26
7.91
-1.93
8.23
-11.74
11.91
-3.50
13.07
-4.01
16.75
-1.96
7.93
6.96
7.28
2.14
-10.23
-3.17
-7.22
-2.73
10.20
12.83
1.98
-5.96
15.93
-3.71
-2.25
6.94
-26.13
6.76
-11.89
-5.76
0.97
48.33
-0.42
19.57
39.80
8.19
14.38
1.74
11.83
21.99
23.61
-2.80
11.80
-21.11
-5.19
-7.63
6.53
-21.14
-14.75
1.20
-11.43
19.93
-6.65
-1.10
-3.90
8.06
-1.02
-6.36
19.33
-84.94
-20.99
-2.16
-11.40
23.96
-15.55
11.20
-18.83
-27.21
-4.55
-6.53
-20.33
-0.21
1.41
-4.90
-15.57
15.06
22.41
-1.33
12.40
-39.31
-10.25
5.64
-25.50
-21.04
-15.65
12.27
-4.70
4.86
-6.52
8.47
-3.07
-21.77
-21.05
-1.30
-19.87
26.86
-23.35
-4.66
26.73
-3.07
-16.29
-7.10
8.83
11.53
8.85
0.90
4.70
-48.31
1.75
-0.06
-23.00
3.79
1.41
3.20
21.20
-39.81
-2.79
-3.63
5.83
15.96
-14.35
14.77
1.23
-21.61
9.71
-8.36
9.50
20.93
26.95
4.54
22.33
33.73
25.15
-3.90
15.30
-73.04
-15.72
7.44
3.00
5.46
7.75
5.34
-10.43
-37.44
-13.69
1.30
-4.47
3.33
-10.09
-2.40
-2.57
-0.01
1.81
4.67
-29.30
13.59
-8.25
-4.16
-22.67
21.03
4.95
-7.63
18.73
-24.07
4.28
11.54
-27.23
44.13
-16.85
2.87
-4.50
1.09
3.85
2.64
-8.03
-25.04
12.45
1.27
9.07
-16.81
6.05
-2.36
3.23
-19.14
15.38
-2.56
-20.97
24.56
26.35
-3.40
-14.03
0.86
6.21
-0.03
32.43
28.26
2.78
-2.66
6.80
7.06
-3.92
-2.43
-17.63
4.86
-15.32
4.84
33.27
-6.64
7.58
5.00
3.10
-0.84
1.91
0.90
26.87
-6.44
-7.32
0.84
-12.37
-15.57
2.61
6.64
4.63
-1.61
4.18
-0.60
2.07
-41.61
-4.09
3.70
2.80
31.73
-13.09
6.80
7.67
8.73
0.08
19.57
-13.67
31.46
4.68
-1.66
-2.93
21.09
-4.85
-3.86
20.43
2.49
28.35
1.74
-16.37
44.69
-8.42
4.17
0.37
11.59
10.25
-6.16
4.47
25.69
-2.89
2.37
-14.40
13.79
-9.29
-1.00
-1.23
-35.01
-1.95
-3.26
3.90
-19.17
-4.95
-1.73
-7.27
25.96
-3.82
0.70
-13.30
19.46
-8.39
-2.63
10.67
20.53
22.35
-2.33
25.17
28.86
29.88
6.17
-9.40
-15.17
-9.35
-3.86
-12.90
28.46
7.25
0.34
-1.23
32.66
12.95
-1.86
-9.20
1.59
0.01
6.14
10.00
48.13
10.38
-2.40
-14.50
4.99
2.68
-4.13
-3.50
5.49
3.81
-7.80
-9.93
29.99
-4.92
-3.10
-24.23
-29.97
-15.39
-4.63
-2.70
-38.21
5.58
8.54
-14.97
9.93
-5.19
2.74
-15.70
-23.57
-6.39
-7.53
-13.73
7.09
-0.95
5.74
16.93
13.36
-7.92
-2.66
31.17
-7.14
2.91
-6.23
-18.07
-52.84
-3.39
-4.86
-10.90
27.79
14.38
-6.46
-17.80
-15.57
0.78
-2.96
-1.67
31.23
9.75
-4.50
-11.27
8.69
-7.59
-5.90
22.37
1.99
16.31
3.20
10.37
6.96
0.65
4.27
-14.13
-41.51
8.68
-4.10
13.83
-7.07
-8.79
3.17
-5.87
8.19
-9.09
4.57
-8.50
-27.74
-5.19
-2.80
-4.70
43.03
2.71
5.30
-0.43
-8.71
-13.05
-1.80
-11.40
-23.17
8.68
3.80
3.30
-41.14
1.95
-0.50
-19.97
-35.84
14.65
-2.16
-1.17
38.69
-12.45
-3.26
0.20
-1.91
-14.89
5.20
22.57
15.86
7.55
-2.66
6.33
-28.37
-3.62
-6.23
-18.67
-8.37
0.78
-0.50
12.00
7.89
6.58
2.77
16.57
11.36
0.91
wave_idl/soi.txt 0100664 0007030 0000404 00000005050 07433253567 0014357 0 ustar 00laurie web 0000166 0000006 $submit
Submit
$dataset
nino3sst
$title
Southern Oscillation Index
$source
Climate Prediction Center (http://nic.fb4.noaa.gov/data/cddb/)
$dataset1
SOI
$units
mb
t0
1896
dt
0.25
$t_units
years
data
0.75
-2.48
-2.41
-1.84
-2.55
-1.61
-0.04
0.16
1.85
0.19
1.19
0.46
2.05
-0.61
-1.01
0.73
-2.75
-0.11
-0.08
-1.84
0.55
0.89
0.29
-2.00
1.51
0.69
-1.48
-1.00
-0.59
1.05
0.66
1.20
2.41
1.22
-0.51
-0.90
-3.89
-4.55
-1.94
-2.37
-1.55
-0.81
1.99
1.54
0.11
-0.41
-0.66
-0.11
-0.76
0.09
1.34
0.06
-0.52
0.35
0.99
0.70
2.05
1.09
2.29
2.36
0.35
-0.88
-1.84
-1.44
-2.69
-1.75
-0.84
-1.04
-0.72
-0.88
-1.21
-1.77
0.11
-1.50
-2.29
-1.50
-3.61
-0.69
1.28
-0.01
-0.49
0.55
2.22
1.63
2.01
2.62
4.69
3.20
1.78
0.92
-1.54
-0.97
-2.89
-0.95
-1.28
-1.97
-0.55
0.05
0.92
0.16
1.54
0.75
-0.04
1.23
1.28
-0.31
0.16
1.26
0.98
0.35
-2.48
-1.00
-0.29
0.12
1.22
1.16
2.08
0.25
-1.68
-2.00
-4.05
-2.31
-2.08
-1.17
1.28
0.79
-0.42
-0.37
0.68
0.02
0.92
1.23
2.51
-0.35
-0.01
1.20
1.18
-0.35
-0.88
-0.00
-0.39
1.59
0.34
-0.26
-0.50
-0.89
-1.57
-0.47
-0.79
0.15
0.12
0.86
0.31
0.39
-1.54
0.46
0.71
-0.21
0.36
0.10
-0.25
0.29
-0.24
-0.84
0.45
0.05
-0.21
-0.17
0.38
1.42
1.92
1.50
2.11
0.22
-0.31
-2.10
-1.15
-1.88
-2.94
-2.87
-2.65
-1.41
-2.64
-2.37
-1.75
0.25
0.49
0.93
1.41
0.25
0.72
-0.10
-0.19
-0.51
-0.28
-0.77
1.38
-0.05
1.12
0.16
-0.22
-1.35
-1.78
-1.34
-0.15
-0.75
1.36
0.46
-0.89
-0.01
-0.74
-0.04
-0.25
-0.78
-0.34
0.26
2.41
1.99
1.96
2.83
0.85
-1.02
-1.75
-1.81
-1.18
0.25
-0.12
-0.81
-0.81
-1.38
-1.82
-0.71
-0.18
0.22
0.68
0.82
0.65
0.82
2.35
2.02
2.05
1.55
1.12
1.55
-0.11
-0.68
-1.05
-1.05
-1.88
-0.45
0.28
-0.91
-1.18
-0.05
-0.55
1.12
0.19
0.35
0.88
0.59
-1.28
0.22
-0.05
0.75
0.55
0.68
0.38
0.62
1.12
-0.12
-0.75
-2.28
-0.08
0.75
1.72
0.42
-0.11
-1.08
-2.62
-1.61
-2.11
-0.68
-0.02
-0.58
2.25
-0.15
0.55
-0.88
0.59
0.92
0.08
-0.41
-1.58
-0.75
-1.28
-0.65
-1.55
0.25
0.52
2.49
2.29
1.32
1.62
1.22
0.75
-1.92
-2.35
-1.81
-1.18
0.45
1.58
2.89
3.65
0.92
1.55
0.32
0.55
1.38
3.35
2.75
2.35
0.05
-2.02
0.19
-0.41
-1.65
-1.98
-2.35
-2.51
0.52
0.32
-0.58
-0.18
0.08
0.45
-1.08
-0.41
-0.92
-0.42
-0.55
-1.25
0.65
0.82
-0.11
0.69
-1.32
-3.38
-4.31
-6.45
-0.68
-0.05
-0.08
-0.08
-0.38
0.12
-0.41
0.45
0.22
0.22
-0.48
-0.35
0.12
-0.68
-1.51
-2.55
-2.82
-2.35
-0.81
-0.48
0.18
2.35
2.35
1.79
1.68
0.35
-0.28
-2.01
0.52
-0.52
-0.58
-0.45
-1.62
-1.48
-2.41
-4.11
-1.35
-0.42
-1.91
-1.85
-1.95
-1.82
-0.88
-0.91
-1.98
-2.85
-2.08
-0.45
-1.08
0.32
-0.58
9
2.36
0.35
-0.88
-1.84
-1.44
-2.69
-1.75
-0.84
-1.04
-0.72
-0.88
-1.21
-1.77
0.11
-1.50
-2.29
-1.50
-3.61
-0.69
1.28
-0.01
-0.49
0.55
2.22
1.63
2.01
2.62
4.69
3.20
1.78
0.92
-1.54
-0.97
-2.89
-0.95
-1.28
-1.97
-0.55
0.05
0.92
0.16
1.54
0.75
-0.04
1.23
1.28
-0.31
0.16
1.26
0.98
0.35
-2.48
-1.00
-0.29
0.12
1.22
1.16
2.08
0.25
-1.68
-2.00
-4.05
-2.31
-2.08
-1.17
1.28
0.79
-0.42
-0.37
0.68
0.02
0.92
1.23
2.51
-0.35
-0.01
1.20
1.18
-0wave_idl/sunspot.txt 0100664 0007030 0000404 00000015771 07433253673 0015311 0 ustar 00laurie web 0000166 0000006 $submit
Submit
$dataset
nino3sst
$title
Wolf's Sunspot Number
$source
National Geophysical Data Center (http://web.ngdc.noaa.gov/stp/SOLAR/solar.html)
$dataset1
Sunspots
$units
#
t0
1748
dt
0.25
$t_units
years
data
0.00
0.00
0.00
0.00
63.53
74.73
79.00
106.43
79.47
92.77
93.20
68.13
52.93
55.93
49.87
31.90
52.00
52.87
44.93
41.40
40.57
35.23
29.67
17.23
1.57
20.37
13.10
13.83
9.40
2.17
9.87
16.83
8.33
11.60
7.27
13.57
20.50
26.97
38.67
43.57
46.20
54.57
48.40
41.23
46.37
48.67
66.50
54.33
67.17
59.53
67.63
57.10
80.57
92.73
95.30
74.80
54.10
59.07
56.67
74.77
40.87
33.80
49.60
56.20
53.20
36.23
29.40
26.57
25.00
20.73
24.23
13.67
19.87
11.93
3.87
9.97
33.47
32.00
35.13
50.70
55.30
65.93
64.73
93.40
67.47
88.23
129.23
139.43
108.87
68.13
113.50
112.70
42.97
112.37
75.87
95.13
74.27
62.40
61.33
67.97
44.93
34.13
23.23
36.80
55.97
41.20
10.67
14.63
5.33
9.13
4.03
9.53
13.20
17.33
13.73
35.00
40.17
85.50
107.73
136.60
140.20
185.17
154.90
137.20
132.80
132.90
122.00
116.00
88.67
96.73
88.57
65.23
75.47
90.23
63.50
43.17
42.83
44.43
38.33
28.23
31.13
25.50
23.40
11.17
10.67
8.33
8.67
13.00
7.83
20.90
29.43
38.23
44.17
78.90
94.73
113.67
109.37
120.40
140.83
157.50
136.83
125.23
139.50
122.07
119.77
122.27
110.67
119.73
108.93
92.67
77.87
80.10
69.57
71.70
60.17
64.90
61.67
66.23
54.93
57.33
55.50
52.10
34.43
45.57
42.33
41.87
30.70
49.10
24.63
22.23
17.37
20.83
20.50
19.80
15.60
8.17
7.53
7.47
5.33
5.23
6.13
0.37
1.80
7.97
11.97
9.07
0.70
5.30
10.03
9.57
17.33
20.97
28.67
31.40
35.73
40.20
45.20
44.00
49.93
41.00
43.43
35.10
42.57
51.10
47.20
39.60
41.97
61.10
52.17
39.00
41.57
36.23
33.77
26.57
26.77
25.33
11.27
15.27
10.13
3.53
1.50
13.10
8.80
9.17
5.77
4.07
0.30
0.00
0.00
0.00
0.00
0.00
0.00
0.00
3.00
2.67
4.63
0.77
7.10
7.30
4.07
11.10
14.00
19.60
13.30
14.83
9.63
17.97
25.87
32.43
38.07
45.23
56.27
48.90
36.60
41.47
63.50
29.20
43.90
27.63
27.57
41.33
28.53
22.77
19.07
24.93
24.13
27.73
16.77
19.83
17.23
8.63
10.47
4.23
3.90
7.80
5.67
6.87
3.33
0.13
0.20
0.00
0.17
6.80
10.83
7.40
7.30
8.67
14.30
11.57
24.10
16.43
24.20
31.17
37.00
52.73
46.60
53.00
48.53
50.13
60.73
82.73
60.33
52.87
54.90
78.77
73.73
60.43
69.63
79.23
52.23
82.57
63.67
42.00
46.03
39.57
47.20
31.63
10.37
20.97
12.67
5.57
8.10
7.77
8.97
6.00
8.07
29.93
17.23
46.10
73.20
90.90
98.13
126.33
106.53
154.83
166.07
135.97
131.03
120.17
123.50
119.57
86.87
82.67
95.27
56.73
116.23
74.47
78.90
61.20
64.17
54.33
27.90
54.47
35.53
29.03
21.40
24.10
19.20
32.07
8.37
13.70
8.50
12.37
12.57
12.03
17.37
17.93
37.53
45.30
30.83
46.60
51.20
64.77
69.47
60.60
64.40
68.47
117.90
142.97
126.50
112.77
124.03
135.50
128.30
88.07
79.80
89.17
83.33
58.57
62.30
62.27
81.83
60.77
53.80
61.63
65.33
55.73
39.77
55.67
40.57
40.77
43.27
31.50
18.70
23.87
18.97
20.80
13.70
6.27
1.17
5.63
1.93
3.90
4.97
6.47
8.77
18.57
27.17
36.40
43.80
41.40
64.03
70.00
87.20
87.93
102.60
97.60
89.87
95.70
103.07
94.53
80.33
81.13
80.13
67.13
57.07
67.37
67.50
44.47
57.13
45.07
34.27
39.60
57.03
44.73
46.00
40.03
42.50
32.50
28.73
18.17
31.53
15.67
9.77
8.20
3.30
3.17
6.53
16.10
19.27
31.47
36.87
62.77
57.83
84.43
73.13
80.57
116.60
157.20
140.73
141.30
118.93
133.20
97.80
94.93
96.00
106.53
104.33
99.50
97.33
56.30
60.73
50.57
57.13
38.27
52.37
30.83
23.30
21.50
9.83
13.43
19.97
3.00
11.30
10.80
15.00
17.20
9.73
7.70
5.90
4.13
1.80
1.90
0.53
4.47
8.10
10.90
23.50
25.70
45.33
34.43
47.03
51.87
62.83
55.50
60.43
68.37
47.83
61.80
50.10
63.30
59.73
81.40
88.63
64.60
56.93
43.87
54.80
70.57
52.03
30.43
37.70
33.83
22.87
7.30
9.23
14.20
17.37
11.40
9.20
6.40
4.90
6.50
5.33
4.37
12.17
3.00
3.67
2.57
12.43
9.53
15.37
36.63
48.53
41.97
64.87
75.17
80.23
71.50
71.23
87.57
98.57
82.97
73.37
93.90
80.73
64.03
63.83
71.97
58.13
61.93
46.13
40.17
44.50
36.43
33.03
20.77
32.50
18.67
34.97
20.87
25.07
25.97
15.60
14.13
8.27
10.43
10.53
14.43
6.97
5.90
2.37
5.33
0.77
2.50
5.97
1.40
3.60
9.23
12.93
19.00
22.60
43.00
31.10
41.47
46.30
48.93
65.70
45.43
62.27
80.47
47.10
58.73
69.13
40.47
81.77
45.30
63.00
58.07
33.93
48.83
72.30
39.10
56.53
30.30
32.57
56.13
26.43
14.30
17.27
16.33
6.73
9.23
3.83
3.00
1.73
4.33
4.27
4.03
1.90
0.30
1.03
2.53
2.83
11.30
8.60
15.63
34.70
47.70
63.57
43.50
55.90
71.33
44.60
56.43
80.47
101.23
134.57
99.30
77.83
72.20
96.40
75.87
64.70
83.70
62.80
43.23
58.40
28.93
27.67
35.57
28.83
29.43
27.50
18.77
30.97
8.27
7.37
10.37
3.10
6.13
5.73
8.13
2.47
18.70
24.17
21.53
15.57
40.67
45.53
75.47
68.07
58.77
58.23
70.47
81.40
77.23
59.03
58.50
80.80
83.00
90.50
56.90
60.63
60.97
56.80
81.03
50.07
34.60
26.30
31.93
29.23
23.70
16.47
15.50
11.30
17.10
6.80
9.37
14.87
3.77
2.70
1.30
5.17
12.57
7.20
9.93
20.73
28.40
35.37
59.63
71.40
66.50
71.77
109.27
114.97
118.77
127.83
96.03
101.37
108.63
123.53
104.67
74.10
109.47
105.33
66.10
64.40
66.33
79.83
60.57
45.50
40.70
64.27
39.47
47.53
32.37
18.37
24.13
22.90
15.93
14.20
12.27
5.07
2.60
12.00
18.70
17.57
32.93
34.47
47.40
70.13
78.03
105.93
115.93
126.30
171.67
172.03
136.03
96.47
177.17
147.80
123.37
152.97
124.97
131.63
130.90
102.03
101.07
75.83
56.77
58.57
100.67
68.53
49.93
28.47
29.63
40.80
26.73
13.47
20.70
17.13
4.10
3.87
0.93
4.90
7.93
16.27
23.97
36.70
74.87
105.33
121.30
157.30
182.90
150.87
180.17
193.67
234.70
186.03
180.93
197.60
173.80
182.07
168.00
164.80
120.13
118.17
117.27
127.67
86.00
52.33
63.27
63.20
36.73
44.87
44.03
31.63
29.87
20.43
36.07
30.53
24.53
16.50
9.07
5.70
9.53
14.47
15.60
12.53
17.63
25.97
47.23
52.70
61.60
105.43
74.43
91.83
102.97
108.63
106.23
107.53
101.17
120.23
110.93
95.37
95.70
114.07
114.60
101.67
88.43
77.00
59.70
64.20
65.70
76.67
77.23
72.43
49.40
44.10
46.53
36.00
25.97
24.97
38.60
43.20
30.87
13.97
8.50
27.27
12.10
11.43
14.47
10.60
13.70
16.07
23.33
31.83
38.70
74.00
92.50
88.90
115.23
147.37
128.47
163.33
181.93
146.93
167.10
142.23
162.33
130.27
124.93
156.60
150.00
142.87
104.87
110.83
106.60
67.27
90.33
68.10
40.83
75.30
64.07
26.20
17.83
16.53
22.63
15.23
17.37
13.60
11.10
9.77
19.13
9.17
30.00
35.20
42.53
58.40
83.30
115.17
143.13
152.60
155.10
157.50
165.97
149.37
125.97
158.30
135.53
148.77
143.67
158.43
132.23
139.27
79.60
71.37
87.70
73.37
57.77
40.83
46.97
41.67
20.63
27.77
29.40
28.40
14.70
13.53
13.37
wave_idl/mauna.txt 0100664 0007030 0000404 00000006533 07433253640 0014665 0 ustar 00laurie web 0000166 0000006 $submit
Submit
$dataset
soi
$title
Mauna Loa Carbon Dioxide
$source
Carbon Dioxide Info. Anal. Center (CDIAC) (http://cdiac.esd.ornl.gov/trends/trends.htm)
$dataset1
CO2
$units
ppm
t0
1958.0
dt
0.08333333
$t_units
years
data
315.56
315.56
315.56
317.29
317.34
316.52
315.69
314.78
313.05
313.11
313.18
314.50
315.42
316.31
316.50
317.56
318.13
318.00
316.39
314.65
313.68
313.18
314.66
315.43
316.27
316.81
317.42
318.87
319.87
319.43
318.01
315.74
314.00
313.68
314.84
316.03
316.73
317.54
318.38
319.31
320.42
319.61
318.42
316.63
314.84
315.16
315.94
316.85
317.78
318.40
319.53
320.42
320.85
320.45
319.45
317.25
316.11
315.27
316.53
317.53
318.58
318.93
319.70
321.22
322.08
321.31
319.58
317.61
316.05
315.83
316.91
318.20
319.41
320.07
320.73
321.39
322.05
321.73
320.27
318.53
316.54
316.72
317.53
318.55
319.27
320.28
320.73
321.97
322.00
321.71
321.05
318.71
317.66
317.14
318.70
319.25
320.46
321.43
322.23
323.54
323.91
323.59
322.24
320.20
318.48
317.94
319.63
320.87
322.17
322.34
322.88
324.25
324.83
323.93
322.38
320.76
319.10
319.24
320.56
321.80
322.40
322.99
323.73
324.86
325.40
325.20
323.98
321.95
320.18
320.09
321.16
322.74
323.83
324.27
325.47
326.50
327.21
326.54
325.72
323.50
322.22
321.62
322.69
323.95
324.89
325.82
326.77
327.97
327.90
327.50
326.18
324.53
322.93
322.90
323.85
324.96
326.01
326.51
327.02
327.62
328.76
328.40
327.20
325.27
323.20
323.40
324.63
325.85
326.60
327.47
327.58
329.56
329.90
328.92
327.88
326.16
324.68
325.04
326.34
327.39
328.37
329.40
330.14
331.33
332.31
331.90
330.70
329.15
327.35
327.02
327.99
328.48
329.18
330.55
331.32
332.48
332.92
332.08
331.01
329.23
327.27
327.21
328.29
329.41
330.23
331.25
331.87
333.14
333.80
333.43
331.73
329.90
328.40
328.17
329.32
330.59
331.58
332.39
333.33
334.41
334.71
334.17
332.89
330.77
329.14
328.78
330.14
331.52
332.75
333.24
334.53
335.90
336.57
336.10
334.76
332.59
331.42
330.98
332.24
333.68
334.80
335.22
336.47
337.59
337.84
337.72
336.37
334.51
332.60
332.38
333.75
334.78
336.05
336.59
337.79
338.71
339.30
339.12
337.56
335.92
333.75
333.70
335.12
336.56
337.84
338.19
339.91
340.60
341.29
341.00
339.39
337.43
335.72
335.84
336.93
338.04
339.06
340.30
341.21
342.33
342.74
342.08
340.32
338.26
336.52
336.68
338.19
339.44
340.57
341.44
342.53
343.39
343.96
343.18
341.88
339.65
337.81
337.69
339.09
340.32
341.20
342.35
342.93
344.77
345.58
345.14
343.81
342.21
339.69
339.82
340.98
342.82
343.52
344.33
345.11
346.89
347.25
346.62
345.22
343.11
340.90
341.18
342.80
344.04
344.79
345.82
347.25
348.17
348.74
348.07
346.38
344.52
342.92
342.62
344.06
345.38
346.11
346.78
347.68
349.37
350.03
349.37
347.76
345.73
344.68
343.99
345.48
346.72
347.84
348.29
349.24
350.80
351.66
351.08
349.33
347.92
346.27
346.18
347.64
348.78
350.25
351.54
352.05
353.41
354.04
353.63
352.22
350.27
348.55
348.72
349.91
351.18
352.60
352.92
353.53
355.26
355.52
354.97
353.75
351.52
349.64
349.83
351.14
352.37
353.50
354.55
355.23
356.04
357.00
356.07
354.67
352.76
350.82
351.04
352.70
354.07
354.59
355.63
357.03
358.48
359.22
358.12
356.06
353.92
352.05
352.11
353.64
354.89
355.88
356.63
357.72
359.07
359.58
359.17
356.94
354.92
352.94
353.23
354.09
355.33
356.63
357.10
358.32
359.41
360.23
359.55
357.53
355.48
353.67
353.95
355.30
356.78
358.34
358.89
359.95
361.25
361.67
360.94
359.55
357.49
355.84
356.00
357.59
359.05
359.98
361.03
361.66
363.48
363.82
363.30
361.93
359.49
358.08
357.77
359.57
360.69
wave_idl/nino3sst.txt 0100664 0007030 0000404 00000006212 07433251306 0015333 0 ustar 00laurie web 0000166 0000006 $submit
Submit
$dataset
nino3sst
$title
NINO3 Sea Surface Temperature
$source
UK Meteorological Office GISST2.3
$dataset1
NINO3 SST
$units
!Uo!NC
t0
1871.0
dt
0.25
$t_units
years
data
-0.15
-0.30
-0.14
-0.41
-0.46
-0.66
-0.50
-0.80
-0.95
-0.72
-0.31
-0.71
-1.04
-0.77
-0.86
-0.84
-0.41
-0.49
-0.48
-0.72
-1.21
-0.80
0.16
0.46
0.40
1.00
2.17
2.50
2.34
0.80
0.14
-0.06
-0.34
-0.71
-0.34
-0.73
-0.48
-0.11
0.22
0.51
0.51
0.25
-0.10
-0.33
-0.42
-0.23
-0.53
-0.44
-0.30
0.15
0.09
0.19
-0.06
0.25
0.30
0.81
0.26
0.10
0.34
1.01
-0.31
-0.90
-0.73
-0.92
-0.73
-0.31
-0.03
0.12
0.37
0.82
1.22
1.83
1.60
0.34
-0.72
-0.87
-0.85
-0.40
-0.39
-0.65
0.07
0.67
0.39
0.03
-0.17
-0.76
-0.87
-1.36
-1.10
-0.99
-0.78
-0.93
-0.87
-0.44
-0.34
-0.50
-0.39
-0.04
0.42
0.62
0.17
0.23
1.03
1.54
1.09
0.01
0.12
-0.27
-0.47
-0.41
-0.37
-0.36
-0.39
0.43
1.05
1.58
1.25
0.86
0.60
0.21
0.19
-0.23
-0.29
0.18
0.12
0.71
1.42
1.59
0.93
-0.25
-0.66
-0.95
-0.47
0.06
0.70
0.81
0.78
1.43
1.22
1.05
0.44
-0.35
-0.67
-0.84
-0.66
-0.45
-0.12
-0.20
-0.16
-0.47
-0.52
-0.79
-0.80
-0.62
-0.86
-1.29
-1.04
-1.05
-0.75
-0.81
-0.90
-0.25
0.62
1.22
0.96
0.21
-0.11
-0.25
-0.24
-0.43
0.23
0.67
0.78
0.41
0.98
1.28
1.45
1.02
0.03
-0.59
-1.34
-0.99
-1.49
-1.74
-1.33
-0.55
-0.51
-0.36
-0.99
0.32
1.04
1.41
0.99
0.66
0.50
0.22
0.71
-0.16
0.38
0.00
-1.11
-1.04
0.05
-0.64
-0.34
-0.50
-1.85
-0.94
-0.78
0.29
0.27
0.69
-0.06
-0.83
-0.80
-1.02
-0.96
-0.09
0.62
0.87
1.03
0.70
-0.10
-0.31
0.04
-0.46
0.04
0.24
-0.08
-0.28
0.06
0.05
-0.31
0.11
0.27
0.26
0.04
0.12
1.11
1.53
1.23
0.17
-0.18
-0.56
0.05
0.41
0.22
0.04
-0.19
-0.46
-0.65
-1.06
-0.54
0.14
0.25
-0.21
-0.73
-0.43
0.48
0.26
0.05
0.11
-0.27
-0.08
-0.10
0.29
-0.15
-0.28
-0.55
-0.44
-1.40
-0.55
-0.69
0.58
0.37
0.42
1.83
1.23
0.65
0.41
1.03
0.64
-0.07
0.98
0.36
-0.30
-1.33
-1.39
-0.94
0.34
-0.00
-0.15
0.06
0.39
0.36
-0.49
-0.53
0.35
0.07
-0.24
0.20
-0.22
-0.68
-0.44
0.02
-0.22
-0.30
-0.59
0.10
-0.02
-0.27
-0.60
-0.48
-0.37
-0.53
-1.35
-1.22
-0.99
-0.34
-0.79
-0.24
0.02
0.69
0.78
0.17
-0.17
-0.29
-0.27
0.31
0.44
0.38
0.24
-0.13
-0.89
-0.76
-0.71
-0.37
-0.59
-0.63
-1.47
-0.40
-0.18
-0.37
-0.43
-0.06
0.61
1.33
1.19
1.13
0.31
0.14
0.03
0.21
0.15
-0.22
-0.02
0.03
-0.17
0.12
-0.35
-0.06
0.38
-0.45
-0.32
-0.33
-0.49
-0.14
-0.56
-0.18
0.46
1.09
1.04
0.23
-0.99
-0.59
-0.92
-0.28
0.52
1.31
1.45
0.61
-0.11
-0.18
-0.39
-0.39
-0.36
-0.50
-0.81
-1.10
-0.29
0.57
0.68
0.78
0.78
0.63
0.98
0.49
-0.42
-1.34
-1.20
-1.18
-0.65
-0.42
-0.97
-0.28
0.77
1.77
2.22
1.05
-0.67
-0.99
-1.52
-1.17
-0.22
-0.04
-0.45
-0.46
-0.75
-0.70
-1.38
-1.15
-0.01
0.97
1.10
0.68
-0.02
-0.04
0.47
0.30
-0.55
-0.51
-0.09
-0.01
0.34
0.61
0.58
0.33
0.38
0.10
0.18
-0.30
-0.06
-0.28
0.12
0.58
0.89
0.93
2.39
2.44
1.92
0.64
-0.24
0.27
-0.13
-0.16
-0.54
-0.13
-0.37
-0.78
-0.22
0.03
0.25
0.31
1.03
1.10
1.05
1.11
1.28
0.57
-0.55
-1.16
-0.99
-0.38
0.01
-0.29
0.09
0.46
0.57
0.24
0.39
0.49
0.86
0.51
0.95
1.25
1.33
-0.00
0.34
0.66
1.11
0.34
0.48
0.56
0.39
-0.17
1.04
0.77
0.12
-0.35
-0.22
0.08
-0.08
-0.18
-0.06
wave_idl/wave_coherency.pro 0100644 0006724 0000404 00000015416 10256306452 0017112 0 ustar 00torrence web 0000166 0000006 ;****************************************************** 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. J. Climate, 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(TEMPORARY(cross_wavelet))^2)/ $
(power1*power2 > 1E-9)
; wave_phase = wave_phase + 360.*(wave_phase LT 0.)
END