; ; Computes and plots the frequency response of various filters. ; ; We want lots of overlapping band-passed filters, whose width increases ; linearly with time scale (from a width of 10 for the 0-10 year timescale, ; to a width of 50 for the 50-100 year timescale, noting that only an 80-yr ; run of data is analysed). ; ; Define filters to apply ; allmin=[findgen(26),findgen(8)*2+26,replicate(40,12)] allmax=[findgen(26)*2+10,findgen(8)*5+65,110.,replicate(120,11)] allts=[replicate(0.,26+8+1),findgen(11)*0.1] nband=n_elements(allmin) print,allmin print print,allmax print print,allts ; ; Prepare for plotting ; loadct,39 multi_plot,nrow=8,ncol=3,layout='large' if !d.name eq 'X' then begin wintim,ysize=850,xsize=700 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=11 endelse def_1color,10,color='blue' def_1color,11,color='red' def_1color,12,color='mlgrey' ; perlist=findgen(148)+3. nper=n_elements(perlist) nt=1000 ; ; Consider each filter separately ; for iband = 0 , nband-1 do begin t1=allmin(iband) t2=allmax(iband) ts=allts(iband) ; ; Generate sine waves at different frequencies and filter these ; fresp=fltarr(nper) for iper = 0 , nper-1 do begin ; y=sin(findgen(nt)*2*!pi/perlist(iper)) ; if t1 gt 0 then begin filter_cru,t1,tsin=y,tshigh=yh1,tslow=yl1 endif else begin yh1=y*0. yl1=y endelse ; if t2 ne 999 then begin filter_cru,t2,tsin=y,tshigh=yh2,tslow=yl2 endif else begin yh2=y yl2=y*0. endelse ; yfilt=yl1-yl2+ts*yl2 ; dummy=moment(y,sdev=sd1) dummy=moment(yfilt,sdev=sd2) fresp(iper)=sd2/sd1 ; endfor ; pause plot,perlist,fresp,$ /xstyle,xrange=[0,150],xtitle='Period (years)',$ /ystyle,yrange=[-0.1,1.1],ytitle='Filter response',$ title=string([t1,t2,ts],format='(2I5,F5.1)') oplot,!x.crange,[0,0],linestyle=1 oplot,!x.crange,[1,1],linestyle=1 oplot,[t1,t1],!y.crange,linestyle=1 oplot,[t2,t2],!y.crange,linestyle=1 ; endfor ; end