C****************************************************************************
C
C WAVEPACK: routines to compute the wavelet transform of a time series,
C and significance levels.
C
C Written by: Christopher Torrence and Gilbert P. Compo
C
C Available from: http://paos.colorado.edu/research/wavelets/
C
C Requires the following packages: CFFTPACK, CHISQR
C
C Notes:
C
C (1) All routines are written in single precision (DOUBLE PRECISION),
C except for CHISQR, which requires double precision input.
C Single precision should be sufficient for most applications.
C
C (2) The CFFTPACK and CHISQR routines were not written by us,
C and no guarentees are made as to their reliability or efficiency.
C
C (3) No provision is made for output to files or to graphics packages.
C The user is expected to call these routines from within
C their own programs. See sample program "wavetest.f".
C
C (4) Little error checking is done. Check your input carefully.
C The programs are not completely ANSI compatible, so use caution.
C
C (5) Time series are currently limited to 65535 points.
C
C
C Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
C Wavelet Analysis. *Bull. Amer. Meteor. Soc.*, 79, 61-78.
C
C Notice: Please acknowledge the use of this software in any publications:
C ``Wavelet software was provided by C. Torrence and G. Compo,
C and is available at URL: http://paos.colorado.edu/research/wavelets/''.
C
C Copyright (C) 1998, Christopher Torrence
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made. This
C routine is provided as is without any express or implied warranties
C whatsoever.
C
C
C Modified: November 1999 by Arjan van Dijk to include IMPLICIT NONE and
C to convert all routines to DOUBLE precision.
C****************************************************************************
C****************************************************************************
C WAVELET: computes the wavelet transform of a time series,
C with appropriate parameters.
C
C
C INPUTS:
C
C n [INT] = the number of points in "y".
C
C y [DOUBLE PRECISION] = the time series of length "n".
C
C dt [DOUBLE PRECISION] = amount of time between each Y value, i.e. the sampling time.
C
C mother [INT] = An integer giving the mother wavelet to use.
C 0='Morlet'
C 1='Paul'
C 2='DOG' (derivative of Gaussian)
C If (mother<0 or >2) then default is 'Morlet'.
C
C param [DOUBLE PRECISION] = mother wavelet parameter. If <0 then default is used.
C For 'Morlet' this is k0 (wavenumber), default is 6.
C For 'Paul' this is m (order), default is 4.
C For 'DOG' this is m (m-th derivative), default is 2.
C
C
C s0 [DOUBLE PRECISION] = the smallest scale of the wavelet. Typically = 2*dt.
C Note: for accurate reconstruction and variance computation
C set s0=dt for Morlet; s0=dt/4 for Paul
C
C dj [DOUBLE PRECISION] = the spacing between discrete scales. Typically = 0.25.
C A smaller # will give better scale resolution, but be slower.
C
C jtot [INT] = the # of scales.
C Scales range from s0 up to s0*2^[(jtot-1)*dj],
C Typically jtot=1+(LOG2(n dt/s0))/dj
C
C npad [INT] = the total number of points (including padding) to
C use for the wavelet transform. Typically this is some
C power of 2. It must be greater or equal to "n".
C If npad>n, then zeroes are padded onto the end
C of the time series.
C
C
C OUTPUTS:
C
C wave [DCMPLX(n,jtot)] = 2D array of the real & imaginary parts
C of the wavelet transform, versus time & scale.
C CABS(wave) gives the WAVELET amplitude,
C ATAN2(AIMAG(wave),DBLE(wave)) gives WAVELET phase.
C The WAVELET power spectrum is CABS(wave)**2.
C
C scale [DOUBLE PRECISION(jtot)] = the wavelet scales that were used.
C
C period [DOUBLE PRECISION(jtot)] = the "Fourier" periods (in time units) corresponding
C to "scale".
C
C coi [DOUBLE PRECISION(n)] = the e-folding factor used for the cone of influence.
C
C
C REQUIRES: WAVE_FUNCTION, CFFTPACK
C
C Copyright (C) 1998, Christopher Torrence and Gilbert P. Compo
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made. This
C routine is provided as is without any express or implied warranties
C whatsoever.
SUBROUTINE WAVELET (n,y,dt,mother,param,s0,dj,jtot,npad,
& wave,scale,period,coi)
IMPLICIT none
INTEGER n,mother,jtot,npad
DOUBLE PRECISION y(n),dt,param,s0,dj,scale(jtot),period(jtot),
& coi(n)
DOUBLE COMPLEX wave(n,jtot)
INTEGER i,j,k,nk
DOUBLE PRECISION ymean,freq1,pi,period1,coi1
C** initialize work arrays
PARAMETER (nk=65535)
DOUBLE PRECISION wsave(4*nk+15),kwave(nk)
DOUBLE COMPLEX yfft(nk),daughter(nk)
pi = 4.D0*ATAN(1.D0)
IF (npad.LT.n) THEN
PRINT*,'**WAVELET: "npad" must be greater than or equal to "n"'
RETURN
END IF
IF ((mother.LT.0).OR.(mother.GT.2)) mother = 0
C** find the time-series mean & remove it
ymean = 0.D0
DO 10 i=1,n
ymean = ymean + y(i)
10 CONTINUE
ymean = ymean/n
DO 15 i=1,n
yfft(i) = y(i) - ymean
15 CONTINUE
C** if desired, pad with extra zeroes
DO 20 i=n+1,npad
yfft(i) = 0.D0
20 CONTINUE
C** find the FFT of the time series [Eqn(3)]
CALL CFFTI(npad,wsave)
CALL CFFTF(npad,yfft,wsave)
DO 30 k=1,npad
yfft(k) = yfft(k)/npad
30 CONTINUE
C** construct the wavenumber array [Eqn(5)]
freq1 = 2.D0*pi/(DBLE(npad)*dt)
kwave(1) = 0.D0
DO 40 i=2,npad/2+1
kwave(i) = (DBLE(i)-1.D0)*freq1
40 CONTINUE
DO 50 i=npad/2+2,npad
kwave(i) = -kwave(npad-i+2)
50 CONTINUE
c**----- Main wavelet loop
c PRINT '(A8,2A12,/,5X,27("-"))','j','scale','period'
DO 100 j=1,jtot
scale(j) = s0*(2.D0**(DBLE(j-1)*dj))
CALL WAVE_FUNCTION(npad,dt,mother,param,scale(j),
& kwave,period1,coi1,daughter)
period(j) = period1
C** multiply the daughter by the time-series FFT
DO 60 k=1,npad
daughter(k) = daughter(k)*yfft(k)
60 CONTINUE
C** inverse FFT [Eqn(4)]
CALL CFFTB(npad,daughter,wsave)
C** store the wavelet transform, discard zero-padding at end
DO 70 i=1,n
wave(i,j) = daughter(i)
70 CONTINUE
c PRINT '(I8,2F12.3)',j,scale(j),period(j)
100 CONTINUE
c**----- end loop
C** construct the cone of influence
DO 110 i=1,(n+1)/2
coi(i) = coi1*dt*(DBLE(i)-1.D0)
coi(n-i+1) = coi(i)
110 CONTINUE
RETURN
END
C****************************************************************************
C WAVE_FUNCTION: computes the daughter wavelets for a particular
C wavelet function, with appropriate parameters.
C
C
C INPUTS:
C
C nk [INT] = the number of points in "kwave"
C
C dt [DOUBLE PRECISION] = amount of time between each Y value, i.e. the sampling time.
C
C mother [INT] = An integer giving the mother wavelet to use.
C 0='Morlet'
C 1='Paul'
C 2='DOG' (derivative of Gaussian)
C
C param [DOUBLE PRECISION] = mother wavelet parameter. If <0 then default is used.
C For 'Morlet' this is k0 (wavenumber), default is 6.
C For 'Paul' this is m (order), default is 4.
C For 'DOG' this is m (m-th derivative), default is 2.
C
C scale1 [DOUBLE PRECISION] = the wavelet scale used to construct the daughter.
C
C kwave [DOUBLE PRECISION(n)] = vector of wavenumbers, used to construct daughter.
C
C
C OUTPUTS:
C
C period1 [DOUBLE PRECISION] = the "Fourier" period (in time units) that corresponds
C to "scale1".
C
C coi1 [DOUBLE PRECISION] = the e-folding factor used for the cone of influence.
C
C daughter [DCMPLX(nk)] = real & imaginary parts of the wavelet function
C at "scale1" and "kwave".
C
C
C REQUIRES: FACTORIAL, CHISQR
C
C
C Reference: Tables 1 & 2 in
C Torrence, C. and G. P. Compo, 1998: A Practical Guide to
C Wavelet Analysis. *Bull. Amer. Meteor. Soc.*, 79, 61-78.
C
C Copyright (C) 1998, Christopher Torrence and Gilbert P. Compo
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made. This
C routine is provided as is without any express or implied warranties
C whatsoever.
SUBROUTINE WAVE_FUNCTION (nk,dt,mother,param,scale1,
& kwave,period1,coi1,daughter)
IMPLICIT none
INTEGER nk,mother
DOUBLE PRECISION dt,kwave(nk),param,scale1,period1,coi1
DOUBLE COMPLEX daughter(nk),norm
DOUBLE PRECISION expnt,sk,pi,fourier_factor
INTEGER k,m,factorial
DOUBLE PRECISION gamma
pi = 4.D0*ATAN(1.D0)
IF (mother.EQ.0) THEN
C******************************************* Morlet wavelet
IF (param.LT.0) param = 6.D0
norm = SQRT(2.D0*pi*scale1/dt)*(pi**(-0.25D0))
DO 10 k=1,nk/2+1
expnt = -0.5D0*(scale1*kwave(k) - param)**2
daughter(k) = DCMPLX(norm*EXP(expnt))
10 CONTINUE
DO 20 k=nk/2+2,nk
daughter(k) = DCMPLX(0.D0)
20 CONTINUE
fourier_factor = (4.D0*pi)/(param + SQRT(2.D0+param**2))
period1 = scale1*fourier_factor
coi1 = fourier_factor/SQRT(2.D0)
ELSE IF (mother.EQ.1) THEN
C******************************************* Paul wavelet
IF (param.LT.0) param = 4.D0
m = INT(param)
norm = SQRT(2.D0*pi*scale1/dt)*
& (2.D0**m/SQRT(DBLE(m*FACTORIAL(2*m-1))))
DO 30 k=1,nk/2+1
expnt = -scale1*kwave(k)
daughter(k) = DCMPLX(norm*(scale1*kwave(k))**m*EXP(expnt))
30 CONTINUE
DO 40 k=nk/2+2,nk
daughter(k) = DCMPLX(0.D0)
40 CONTINUE
fourier_factor = (4.D0*pi)/(2.D0*m + 1.D0)
period1 = scale1*fourier_factor
coi1 = fourier_factor*SQRT(2.D0)
ELSE IF (mother.EQ.2) THEN
C******************************************* DOG wavelet
IF (param.LT.0) param = 2.D0
m = INT(param)
norm = SQRT(2.D0*pi*scale1/dt)*SQRT(1.D0/GAMMA(m+0.5D0))
norm = -norm*(DCMPLX(0.D0,1.D0)**m)
DO 50 k=1,nk
sk = scale1*kwave(k)
daughter(k) = norm*(sk**m)*EXP(-0.5D0*sk**2)
50 CONTINUE
fourier_factor = 2.D0*pi*SQRT(2.D0/(2.D0*m+1.D0))
period1 = scale1*fourier_factor
coi1 = fourier_factor/SQRT(2.D0)
ELSE
stop
END IF
RETURN
END
C****************************************************************************
C FACTORIAL: compute the factorial (n!) of an integer n
C Copyright (C) 1998, Christopher Torrence
FUNCTION FACTORIAL(n)
IMPLICIT NONE
INTEGER factorial,n,i
factorial = 1
DO 10 i=1,n
factorial = factorial*i
10 CONTINUE
END
C****************************************************************************
C WAVE_SIGNIF: computes the significance levels for a wavelet transform.
C
C
C INPUTS:
C
C isigtest [INT] = 0, 1, or 2.
C
C If 0, then just do a regular chi-square test,
C i.e. Eqn (18) from Torrence & Compo.
C If 1, then do a "time-average" test, i.e. Eqn (23).
C In this case, DOF(j) should be set to NA, the number
C of local wavelet spectra that were averaged together
C at each scale. For the Global Wavelet Spectrum,
C this would be dof(j)=N-scale(j),
C where N is the number of points in your time series.
C If 2, then do a "scale-average" test, i.e. Eqns (25)-(28).
C In this case, "dof(1)" and "dof(2)" should be set to the
C smallest (S1) and largest (S2) scales that were averaged
C together, respectively.
C e.g. if you scale-averaged scales between 2 and 8,
C then dof(1)=2.0 and dof(2)=8.0
C
C
C n [INT] = the number of points in "y".
C
C y [DOUBLE PRECISION] = the time series of length "n".
C
C dt [DOUBLE PRECISION] = amount of time between each Y value, i.e. the sampling time.
C
C mother [INT] = An integer giving the mother wavelet to use.
C 0='Morlet'
C 1='Paul'
C 2='DOG' (derivative of Gaussian)
C
C param [DOUBLE PRECISION] = mother wavelet parameter.
C
C s0 [DOUBLE PRECISION] = the smallest scale of the wavelet.
C
C dj [DOUBLE PRECISION] = the spacing between discrete scales.
C
C jtot [INT] = the # of scales.
C
C scale [DOUBLE PRECISION(jtot)] = the wavelet scales that were used.
C
C period [DOUBLE PRECISION(jtot)] = the "Fourier" periods corresponding to "scale".
C
C lag1 [DOUBLE PRECISION] = lag 1 Autocorrelation, used for signif levels.
C Default is 0.0, which corresponds to white-noise.
C
C siglvl [DOUBLE PRECISION] = significance level to use. Default is 0.05 (the "5%" level)
C
C dof [DOUBLE PRECISION(jtot)] = degrees-of-freedom for signif test.
C IF SIGTEST=0, then the input dof is ignored.
C IF SIGTEST=1, then dof(j) = NA, the number of times averaged together.
C IF SIGTEST=2, then dof(1)=S1, dof(2)=S2, the range of scales averaged.
C
C
C OUTPUTS:
C
C dof [DOUBLE PRECISION(jtot)] = degrees-of-freedom that were actually used.
C IF SIGTEST=0, then dof(j) = 2 (or 1 for the 'DOG')
C IF SIGTEST=1, then dof(j) = degrees-of-freedom versus scale.
C IF SIGTEST=2, then dof(1)=degrees-of-freedom, dof(2...jtot)=0.0
C
C fft_theor [DOUBLE PRECISION(jtot)] = theoretical red-noise spectrum vs scale.
C IF SIGTEST=2, then fft_theor(1) = the average spectrum from S1-->S2
C fft_theor(2...jtot) = 0.0
C
C signif [DOUBLE PRECISION(jtot)] = significance levels vs scale.
C IF SIGTEST=2, then signif(1) = the significance level
C signif(2...jtot) = 0.0
C
C ymean [DOUBLE PRECISION] = the mean of the time series.
C
C variance [DOUBLE PRECISION] = the variance of the time series.
C
C Cdelta [DOUBLE PRECISION] = the constant "Cdelta" for the mother wavelet (Table 2).
C
C psi0[DOUBLE PRECISION] = the constant 'psi(0)' for the mother wavelet (Table 2)
C
C REQUIRES: CHISQR
C
C Copyright (C) 1998, Christopher Torrence and Gilbert P. Compo
C This software may be used, copied, or redistributed as long as it is not
C sold and this copyright notice is reproduced on each copy made. This
C routine is provided as is without any express or implied warranties
C whatsoever.
SUBROUTINE WAVE_SIGNIF (isigtest,n,y,dt,mother,param,dj,jtot,
& scale,period,lag1,siglvl,dof,fft_theor,signif,
& ymean,variance,Cdelta,psi0)
IMPLICIT none
INTEGER isigtest,n,mother,jtot
DOUBLE PRECISION y(n),dt,param,dj,scale(jtot),period(jtot)
DOUBLE PRECISION lag1,siglvl,dof(jtot),fft_theor(jtot),
& signif(jtot)
DOUBLE PRECISION ymean,variance,Cdelta,psi0
INTEGER i,j,m,status,javg1,javg2,navg
DOUBLE PRECISION pi,freq1,dofmin,gammafac,dj0,Savg,Smid
DOUBLE PRECISION fft_theor1
DOUBLE PRECISION chisqr,p,q,bound
pi = 4.D0*ATAN(1.D0)
IF (siglvl.LE.0.) siglvl = 0.05D0
IF (lag1.LE.0.D0) lag1 = 0.D0
Cdelta = -1.D0
gammafac = -1.D0
dj0 = -1.D0
psi0 = -1.D0
IF (mother.EQ.0) THEN
C******************************************* Morlet wavelet
dofmin = 2.D0
IF (param.EQ.6.D0) THEN
Cdelta = 0.776D0
gammafac = 2.32D0
dj0 = 0.60D0
psi0 = pi**(-0.25D0)
END IF
ELSE IF (mother.EQ.1) THEN
C******************************************* Paul wavelet
m = INT(param)
dofmin = 2.D0
IF (m.EQ.4) THEN
Cdelta = 1.132D0
gammafac = 1.17D0
dj0 = 1.5D0
psi0 = 1.079D0
END IF
ELSE IF (mother.EQ.2) THEN
C******************************************* DOG wavelet
m = INT(param)
dofmin = 1.D0
IF (m.EQ.2) THEN
Cdelta = 3.541D0
gammafac = 1.43D0
dj0 = 1.4D0
psi0 = 0.867D0
ELSE IF (m.EQ.6) THEN
Cdelta = 1.966D0
gammafac = 1.37D0
dj0 = 0.97D0
psi0 = 0.884D0
END IF
ELSE
stop
END IF
C** find the time-series variance
ymean = 0.D0
DO 10 i=1,n
ymean = ymean + y(i)
10 CONTINUE
ymean = ymean/n
variance = 0.D0
DO 15 i=1,n
variance = variance + (y(i) - ymean)**2
15 CONTINUE
variance = variance/(DBLE(n)) ! - 1.D0)
C** construct theoretical red(white)-noise power spectrum [Eqn(16)]
DO 20 j=1,jtot
freq1 = dt/period(j)
fft_theor(j) = variance*(1.D0-lag1**2)/
& (1.D0 - 2.D0*lag1*COS(freq1*2.D0*pi) + lag1**2)
20 CONTINUE
q = DBLE(siglvl)
p = 1d0 - q
IF (isigtest.EQ.0) THEN
C******************************************* no smoothing, dof=dofmin
C see Eqn(18)
DO 30 j=1,jtot
dof(j) = dofmin
CALL CDFCHI(2,p,q,chisqr,DBLE(dofmin),status,bound)
signif(j) = fft_theor(j)*chisqr/dofmin
30 CONTINUE
ELSE IF (isigtest.EQ.1) THEN
C*********************************** time-averaged, dof depend on scale
IF (gammafac.LE.0.D0) THEN
PRINT*,'**WAVE_SIGNIF: "gammafac" undefined for this wavelet'
RETURN
END IF
DO 40 j=1,jtot
IF (dof(j).LT.1.) dof(j) = 1.D0
C see Eqn(23)
dof(j) = dofmin*SQRT(1.D0+(dof(j)*dt/gammafac/scale(j))**2)
IF (dof(j).LT.dofmin) dof(j) = dofmin
CALL CDFCHI(2,p,q,chisqr,DBLE(dof(j)),status,bound)
signif(j) = fft_theor(j)*chisqr/dof(j)
40 CONTINUE
ELSE IF (isigtest.EQ.2) THEN
C*********************************** scale-averaged, dof depend on scale
IF (Cdelta.LE.0.) THEN
PRINT*,'**WAVE_SIGNIF: "Cdelta" and "dj0" '//
& 'undefined for this wavelet'
RETURN
END IF
javg1 = 0
javg2 = 0
DO 50 j=1,jtot
IF ((scale(j).GE.dof(1)).AND.(javg1.EQ.0)) javg1 = j
IF (scale(j).LE.dof(2)) javg2 = j
50 CONTINUE
IF ((javg1.EQ.0).OR.(javg2.EQ.0).OR.(javg1.GT.javg2)) THEN
PRINT*,'**WAVE_SIGNIF: Scales in "dof(1)" & "dof(2)" '//
& 'are out of range.'
RETURN
END IF
navg = javg2 - javg1 + 1
C see Eqn(25)
Savg = 0.D0
DO 60 j=javg1,javg2
Savg = Savg + 1.D0/scale(j)
60 CONTINUE
Savg = 1.D0/Savg
C see Eqn(27)
fft_theor1 = 0.D0
DO 70 j=javg1,javg2
fft_theor1 = fft_theor1 + fft_theor(j)/scale(j)
70 CONTINUE
fft_theor(1) = Savg*fft_theor1
C see Eqn(28)
Smid = EXP(0.5D0*(LOG(scale(javg1)) + LOG(scale(javg2))))
dof(1) = (dofmin*navg*Savg/Smid)*SQRT(1 + (navg*dj/dj0)**2)
C see Eqn(26)
CALL CDFCHI(2,p,q,chisqr,DBLE(dof(1)),status,bound)
signif(1)=(dj*dt/Cdelta/Savg)*fft_theor(1)*chisqr/dof(1)
DO 80 j=2,jtot
dof(j) = 0.D0
fft_theor(j) = 0.D0
signif(j) = 0.D0
80 CONTINUE
ELSE
stop
END IF
END