PRO compute_fft, dat, init_ind, trend_coef,main_title=main_title, frequency=frequency, key_lowpass=key_lowpass, key_highpass= key_highpass,key_smooth=key_smooth ; This IDL program allows finding the dominant frequency signal(s) that are ; present in a time series data set. ; The part between BEGIN - END is an add-on which removes the trend of the data ; set and does the spectral analysis on a detrended data ; The second BEGIN - END helps empasize on specific frequencies. ; ; key_smooth can be set to 2, 3, 4,... to smooth out the data with 2, 3, 4, ... moving average. ; nl= N_ELEMENTS(dat) IF NOT KEYWORD_SET(frequency) THEN frequency= 30./(60.*24) ; per day sample= 1/ frequency ; Number of data per day x= FINDGEN(nl)*frequency y= dat ; ------------------------ BEGIN ---------------------------- ; --------------- Extract and reprocess the data ------------ ; Just take a sample of the data ;; init_ind= 800 xn=x[init_ind:*] & yn= y[init_ind:*] nln= nl- init_ind x= xn & y= yn & nl= nln ; Remove any trend from the data ;; trend_coef= .02 y= y- trend_coef* x ; ------------------------- END ----------------------------- ; ------------------------- BEGIN -------------------------- ; ---------- Empasize on specific frequencies -------------- period= 1. ; Day c03= .0 & c05= .00 & c065= 0.0 & c1= 0 & c2= 0 & c10= 0 & c30= 0 y1= c03*SIN(2.*!PI*x/(period*0.3))+ $ c05* SIN(2.*!pi*x/(period*0.5))+ $ c065* SIN(2.*!pi*x/(period*0.65))+ $ c1* SIN(2.*!PI*x/(period))+ $ c2* SIN(2.*!PI*x/(2.* period))+ $ c10* SIN(2.*!pi*x/(period*10))+ $ c30*SIN(2.*!pi*x/(period*30)) ; The above data will have 3 frequencies of period, period*10, and period*36. ; ; The first component will have pic at 1.E-1 (= 1/0.1 days) ; The second at 1.E-2 (= 1/0.01 days) ; The third at 0.00277= 2.77E-3 (= 360. days) y= y+y1 ;& y= y1 ; -------------------------- END --------------------------- ; ------------- SMOOTH ----------- IF KEYWORD_SET(key_smooth) THEN BEGIN IF key_smooth NE 1 THEN y= SMOOTH(y,key_smooth) ENDIF ; --------- ; ------------- LOWPASS or HIGHPASS filter ------------- fy= FINDGEN(nl) IF (nl MOD 2) THEN fy(nl/2:*)= -REVERSE(fy(0:nl/2)) ELSE fy(nl/2:*)= -REVERSE(fy(1:nl/2)) cut= 3./100* nl cut= 3./100* nl filter=1./(1+(fy/cut)^10) lowpass=FFT(FFT(y,1)*filter,-1) highpass=fft(fft(y,1)*(1.0-filter),-1) IF KEYWORD_SET(key_highpass) THEN BEGIN y= highpass ; lowpass ; ENDIF IF KEYWORD_SET(key_lowpass) THEN BEGIN y= lowpass ; ENDIF ; ------------------------------------------------------ !p.multi= [0,1,2] PLOT, x, y, xtitle='Day nb',title=main_title;, xrange=[0,100] v= FFT(y) v_h= FFT(HANNING(nl)*y) f= FINDGEN(nl/2+ 1) / (nl*(1./sample)) ;PLOT, f, ABS(v_h(0:nl/2))^2, ytitle= 'Power spectrum' $ PLOT, f, ABS(v(0:nl/2))^2, ytitle= 'Power spectrum of u(k)' $ , title= 'With Hanning window (black solid) !Cand without window (red solid)' $ , /ylog, xtitle= 'Frequency in cycles/day', /xlog, xstyle=1,/nodata, xrange=[1./nl, nl/2.] OPLOT, f, ABS(v(0:nl/2))^2, color=fsc_color('red') OPLOT, f, ABS(v_h(0:nl/2))^2 PRINT, [1./nl, nl/2.], ALOG10([1./nl, nl/2.]) !p.multi= 0 END