; See "http://paos.colorado.edu/research/wavelets/" ; Written January 1998 by C. Torrence ; ; Reads in the combined Hugershoff and Age-banded ALL reconstruction ; and performs evolutionary wavelet analysis on it. ; donh=1 ; 0=NH 1=ALL if donh eq 0 then begin regtit='NH' endif else begin regtit='ALL' endelse ; ; Get the combined reconstruction ; restore,filename='combtemp'+regtit+'_calibrated.idlsave' ; ; Remove missing data from start & end (end in 1960 due to decline) ; kl=where((yrmxd ge 1402) and (yrmxd le 1960),n) sst=prednh(kl) ;------------------------------------------------------ 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 = 1. ; annual data time = FINDGEN(n)*dt + 1402.0 ; construct time array xrange = [1402,1960] ; plotting range pad = 1 s0 = 2.*dt ; this says start at a scale of 1 year 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' mothper=6 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, $ param=mothper,$ 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 loadct,39 multi_plot,nrow=3,layout='large' if !d.name eq 'X' then begin !P.FONT = -1 WINDOW,0,XSIZE=600,YSIZE=800 endif else begin !P.FONT = 0 endelse !X.STYLE = 1 !Y.STYLE = 1 !x.omargin=[0,2] def_1color,19,color='white' def_1color,20,color='vlgrey' def_1color,21,color='lgrey' def_1color,22,color='mgrey' def_1color,23,color='mdgrey' def_1color,24,color='black' ;--- Plot time series ; pos1 = [0.1,0.75,0.7,0.95] ; PLOT,time,sst,XRANGE=xrange, $ ; XTITLE='Time (year)',YTITLE='(!Uo!NC wrt 1961-90)', $ ; TITLE='a) ALL reconstructed Temperature', $ ; POSITION=pos1 ; IF (N_ELEMENTS(recon_sst) GT 1) THEN OPLOT,time,recon_sst,COLOR=144*cfac ; 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 = [256,2.] ; years levels = [0.,0.05,0.1,0.2,0.4,0.8,100.] colors = [19,20,21,22,23,24] period2 = FIX(ALOG(period)/ALOG(2)) ; integer powers of 2 in period ytickv = 2.^(period2(UNIQ(period2))) ; unique powers of 2 kl=where((ytickv le yrange(0)) and (ytickv ge yrange(1))) ytickv=ytickv(kl) ; pos2 = [pos1(0),0.35,pos1(2),0.65] scale_vert,levels=levels,c_colors=colors,/noextremes,$ ytitle='Wavelet power spectrum (!Uo!NC!U2!N)' CONTOUR,power,time,period,$ ;POSITION=pos2, $ XRANGE=xrange,YRANGE=yrange,/ystyle,/ylog, $ 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.05,0.1,0.2,0.4,0.8!Uo!NC!U2!N)' CONTOUR,power,time,period,/overplot,$ LEVELS=levels ; significance contour, levels at -99 (fake) and 1 (significant) CONTOUR,signif,time,period,/OVERPLOT,LEVEL=1,THICK=12, $ C_LABEL=0,color=!p.background CONTOUR,signif,time,period,/OVERPLOT,LEVEL=1,THICK=6 ; cone-of-influence, anything "below" is dubious x = [time(0),time,MAX(time)] y = [MAX(period),coi,MAX(period)] color = !p.background POLYFILL,x,y,COLOR=color,NOCLIP=0 POLYFILL,x,y,COLOR=color,NOCLIP=0 PLOTS,time,coi,NOCLIP=0,THICK=1 scale_vert,/off ;--- Plot global wavelet spectrum ; pos3 = [0.74,pos2(1),0.95,pos2(3)] ; blank = REPLICATE(' ',29) ; PLOT,global_ws,period,/NOERASE,POSITION=pos3, $ ; xrange=[0,1],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,0.2,8,'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 ; END