;******************************************************************* 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. ; Its units are sigma^2 (the time series variance). ; ; ; 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 ; ; DOF = degrees-of-freedom for signif test. Default is 2 (1 for DOG), which ; indicates no smoothing. If one has averaged over N times, ; then /GLOBAL should be set, and DOF should be set to N (*not* 2N). ; The actual DOFs are then computed using Eqn(23). ; ; GLOBAL = if set, then assume that wavelet has been globally smoothed, or ; averaged in time (see DOF above), and use Eqn(23). ; ; VERBOSE = if set, then print out info for each analyzed scale. ; ; ; 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 red-noise spectrum as fn 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-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/''. ; ; Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to ; Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, xx-xx. ; ; Please send a copy of such publications to either C. Torrence or G. Compo: ; Dr. Christopher Torrence Dr. Gilbert P. Compo ; Advanced Study Program Campus Box 311 ; National Center for Atmos. Research Program in Atmos. and Oceanic Sciences ; P.O. Box 3000 University of Colorado at Boulder ; Boulder CO 80307--3000, USA. Boulder CO 80309-0311, USA. ; E-mail: torrence@ucar.edu E-mail: compo@monsoon.colorado.edu ;---------------------------------------------------------------------------- ;- FUNCTION morlet,k0,a,k,period,coi,dofmin,gamma ;*********************** MORLET IF (k0 EQ -1) THEN k0 = 6d n = N_ELEMENTS(k) expnt = -(a*k - k0)^2/2d*(k GT 0.) norm = SQRT(a*k(1))*(!PI^(-0.25))*SQRT(n) ; 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 = a*fourier_factor coi = fourier_factor/SQRT(2) ; Cone-of-influence [Sec.3g] dofmin = 2 ; Degrees of freedom with no smoothing IF (k0 EQ 6) THEN gamma = 2.32 ELSE gamma=-1 ; decorrelation factor ; PRINT,a,SQRT(TOTAL(morlet^2,/DOUBLE))*SQRT(n) RETURN,morlet END FUNCTION paul,m,a,k,period,coi,dofmin,gamma ;************************** PAUL IF (m EQ -1) THEN m = 4d n = N_ELEMENTS(k) expnt = -(a*k)*(k GT 0.) norm = SQRT(a*k(1))*(2^m/SQRT(m*FACTORIAL(2*m-1)))*SQRT(n) paul = norm*((a*k)^m)*EXP(expnt > (-100d))*(expnt GT -100) paul = paul*(k GT 0.) fourier_factor = 4*!PI/(2*m+1) period = a*fourier_factor coi = fourier_factor dofmin = 2 ; Degrees of freedom with no smoothing IF (m EQ 4) THEN gamma=1.17 ELSE gamma=-1 ; decorrelation factor ; PRINT,a,N_ELEMENTS(k),SQRT(TOTAL(paul^2,/DOUBLE))*SQRT(n) RETURN,paul END FUNCTION dog,m,a,k,period,coi,dofmin,gamma ;*************************** DOG IF (m EQ -1) THEN m = 2 n = N_ELEMENTS(k) expnt = -(a*k)^2/2d norm = SQRT(a*k(1)/GAMMA(m+0.5))*SQRT(n) I = DCOMPLEX(0,1) gauss = -norm*(I^m)*(a*k)^m*EXP(expnt > (-100d))*(expnt GT -100) fourier_factor = 2*!PI*SQRT(2./(2*m+1)) period = a*fourier_factor coi = fourier_factor dofmin = 1 ; Degrees of freedom with no smoothing gamma = -1 IF (m EQ 2) THEN gamma=1.43 ; decorrelation factor IF (m EQ 6) THEN gamma=1.37 ; decorrelation factor ; PRINT,a,N_ELEMENTS(k),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, $ ;*** optional inputs SCALE=scale,PERIOD=period,YPAD=ypad, $ ;*** optional outputs DAUGHTER=daughter,COI=coi, $ ;*** optional outputs LAG1=lag1,SIGLVL=siglvl,DOF=dof,GLOBAL=global, $ ;*** optional inputs SIGNIF=signif,FFT_THEOR=fft_theor, $ ;*** optional outputs YFFT = yfft, $ ;*** optional output NO_WAVE=no_wave, $ ;*** optional inputs OCT=oct,VOICE=voice ;*** defunct inputs r = CHECK_MATH(0,1) n = N_ELEMENTS(y1) n1 = n ;....check keywords & optional inputs IF (N_ELEMENTS(s0) LT 1) THEN s0 = 2*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 = oct/dj IF (N_ELEMENTS(J) LT 1) THEN J=FIX((ALOG(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 verbose = KEYWORD_SET(verbose) do_daughter = KEYWORD_SET(daughter) do_wave = NOT KEYWORD_SET(no_wave) ;....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 if pad eq 1 then base2 = FIX(ALOG(n)/ALOG(2) + 0.4999)+1 ; power of 2 nearest to N if pad eq 2 then base2= FIX(ALOG(n)/ALOG(2))+1 ; next highest power of 2base2-1 ypad = [ypad,FLTARR(2L^(base2) - 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)/n/dt*2*!PI 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=',J IF (n1 NE n) THEN PRINT,'(padded with ',n-n1,' zeroes)' PRINT,['a1','a','period','totvar','mathflag'], $ FORMAT='(/,A2,A6,A10,A9,A9)' ENDIF ;verbose ;....loop thru each SCALE FOR a1=0,na-1 DO BEGIN ;scale psi_fft=CALL_FUNCTION(mother,param,scale(a1),k,period1,coi,dofmin,gamma) IF (do_wave) THEN $ wave(*,a1) = FFT(yfft*psi_fft,1,/DOUBLE) ;wavelet transform[Eqn(4)] period(a1) = period1 ; save period 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='(I2,f,f,f,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,*)] ;....confidence 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 = STDEV(y1)^2*fft_theor ; include time-series variance signif = fft_theor IF (N_ELEMENTS(dof) LT 1) THEN dof = dofmin IF (NOT(KEYWORD_SET(global))) THEN BEGIN ; no smoothing, DOF=2 signif = fft_theor*CHISQR_CVF(1. - siglvl,dof)/dof ; [Eqn(18)] ENDIF ELSE BEGIN ; smoothing, DOFs depend upon scale [Sec.5] 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(na) + dof dof = dof > 1 dof = dofmin*SQRT( 1 + (dof*dt/gamma/scale)^2 ) ; [Eqn(23)] dof = dof > dofmin ; minimum DOF is dofmin FOR a1=0,na-1 DO BEGIN chisqr = CHISQR_CVF(1. - siglvl,dof(a1))/dof(a1) signif(a1) = fft_theor(a1)*chisqr ENDFOR ENDELSE RETURN,wave(0:n1-1,*) ; get rid of padding before returning END