PRO Periodogram,t,x,nu,power,freq_range=freq_range,fap=fap,power_sign=power_sign ; This subroutine performs a periodogram analysis of a signal x(t). ; ; INPUT PARAMETERS: ; x is a 1D array of the signal data points; ; t is a 1D array of the times in s when the data points were measured; ; ; freq_range is an (optional) string keyword: ; 'fft' - the frequency range coincides with the natural frequencies ; given by FFT; - DEFAULT ; 'fft_10' - the frequency range is 10 times denser than the natural ; frequencies given by FFT ; This option can ; fap is an (optional) keyward prescribing the value of the ; false alarm probability; - DEFAULT fap=0.01 ; ; OUTPUT PARAMETERS: ; power - is a 1D array giving the spectral density as ; a funtion of the frequency nu; The power is normalised to the total variance. ; of the input signal. The size of the array depends upon the keyward freq_range ; nu - is a 1D array giving the values of the frequency of the spectrum; ; The size of the array depends upon the keyward freq_range. ; power_sign - the power threshold corresponding to the given false alarm probability fap. ; ; CALLING SEQUENCE: ; Before the use, the subroutine must be compiled by the command: ; .r periodogram ; IDL> Periodogram,t,x,nu,power,power_sign=power_sign,fap=.1,freq_range='fft' ; IDL> plot,nu,power,thick=3,yrange=[0,20.] ; IDL> oplot,[0.,max(nu)],[power_sign,power_sign] ; ; HISTORY: ; 15/07/2007 - was created from per_sub.pro by V.M.Nakariakov (V.Nakariakov@warwick.ac.uk) ; 25/02/2008 - the mismatch between the output array sizes fixed. ; 27/02/2008 - the parameter /DOUBLE is used in variance ; ; defaults IF n_elements(freq_range) EQ 0 THEN freq_range = 'fft' IF n_elements(fap) EQ 0 THEN fap = 0.01 sz = size(x) npt = sz[1] sz = size(t) IF (npt NE sz[1]) THEN BEGIN message,'Input array size mismatch' GOTO,nafik ENDIF Tdur = t(npt-1)-t(0) ; (Scargle periodogram is time-shift invariant) ti = Tdur/(npt-1) ; the cadence time CASE freq_range OF 'fft': begin npw = npt/2 ; the number of points in output spectrum omega = dblarr(npw) ; the frequency power = dblarr(npw) ; the spectral density tau = dblarr(npw) ; the Scargle's parameter freq = findgen(npw)/(2*npw*ti) ; old version: findgen(npw+1)/(2*npw*ti) omega = 2.D0*!DPI*freq end 'fft_10':begin N_ext = 10 npw = npt/2 *N_ext ; the number of points in output spectrum omega = fltarr(npw) ; the frequency power = fltarr(npw) ; the spectral density tau = fltarr(npw) ; the Scargle's parameter freq = findgen(npw)/(2*npw*ti) ; old version: findgen(npw+1)/(2*npw*ti) omega = 2.D0*!DPI*freq ;help,omega end ENDCASE for i = 1,npw-1 do begin ssin2 = 0.D0 scos2 = 0.D0 for j=0, npt-1 do begin ssin2 = ssin2 + sin(2.D0*omega[i]*t[j]) scos2 = scos2 + cos(2.D0*omega[i]*t[j]) endfor tau(i) = 1.D0/(2.D0*omega[i])*atan(ssin2/scos2) ;print,tau sxcos = 0.D0 sxsin = 0.D0 scos2 = 0.D0 ssin2 = 0.D0 for j=0, npt-1 do begin sxcos = sxcos + x[j]* cos(omega[i]*(t[j]-tau[i])) sxsin = sxsin + x[j]* sin(omega[i]*(t[j]-tau[i])) scos2 = scos2 + cos(omega[i]*(t[j]-tau[i]))*cos(omega[i]*(t[j]-tau[i])) ssin2 = ssin2 + sin(omega[i]*(t[j]-tau[i]))*sin(omega[i]*(t[j]-tau[i])) endfor power(i) = 1.D0/2.D0 * (sxcos^2/scos2 + sxsin^2/ssin2) endfor nu = omega/(2.D0*!DPI) power = power/variance(x, /DOUBLE) ; Scargle's false alarm probability: horne = long(-6.362+1.193*npt+0.00098*npt^2.) power_sign = -alog( 1.D0 - ((1.D0-fap)^(1./horne)) ) nafik: end