;****************************************************** 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(cross_wavelet)^2)/(power1*power2 > 1E-9) ; wave_phase = wave_phase + 360.*(wave_phase LT 0.) END