pro varadi,file,rmax,M,K,filters,seed,S_lag,eigen_values ;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; NAME: ; varadi ; CALLING SEQUENCE: ; varadi,file,rmax,M,K,filters,seed,S_lag,eigen_values ; PURPOSE: ; Compute the RLSSA from Varadi et al, ApJ, 1999, 526, 1052 ; INPUTS: ; file filename ; rmax number of largest eigenvalues for reconstruction ; M number of random lags ; K maximum lag ; filters reconstructed filters of size fltarr(K) ; seed dummy argument ; S_lag return set of lags fltarr(M) ; eigen_values return eigenvalues fltarr(M) ; OPTIONAL KEYWORDS: ; none ; OUTPUTS: ; none ; OPTIONAL OUTPUT: ; none ; EXAMPLE: ; thomson,file ; LIMITATIONS: ; None yet ; COMMONS: ; None ; PROCEDURES USED: ; fft, readfits ; MODIFICATION HISTORY: ; Beta 1: WRITTEN, Thierry Appourchaux, November 14, 2000 ; Beta 2: Return also eigenvalues, November 17, 2000 ;------------------------------------------------------------------------------ ; ; Read file print,'r_max (10 is a good number)=',rmax print,'M (number of different lags, 250 is OK)',M print,'K (maximum lag, 25000 is OK)',K a=readfits(file,/silent) ; a=randomn(seed,131072)+cos(2.*!pi*findgen(131072)/3.) m_2=moment(a,/double) a=(a-m_2(0))/sqrt(m_2(1)) ; a=randomn(seed,131072)*1.d0 ; ; N_a=N_elements(a) ; k_lag=K ; ; Here we assume that the time series is sampled at 60 s ; xfreq=findgen(N_a)*1.e06/60./N_a a_rev=reverse(a) ; ; Pad the time series with zero ; a_new=dblarr(2*N_a) a_new_rev=dblarr(2*N_a) a_new(0:N_a-1)=a a_new_rev(0:N_a-1)=a_rev ; ; Then compute the correlation using fft (55) ; Please note the use of the different normalization ; It was calibrated using the real formula (53) ; q=2.d0*N_a*fft(fft(a_new,-1,/double)*fft(a_new_rev,-1,/double),1,/double) ; ; Shift to put the max at 0 ; q=shift(double(q),-(N_a-1)) q=q(0:k_lag-1) print,q(0:10) ; ; Calibration ; ; ; q_clas=fltarr(k_lag) ; for k=0,k_lag-1 do begin ; q_clas(k)=total(a(0:N_a-1-k)*a(k:N_a-1)) ; endfor ; ; Equation (57) ; q=q/double(N_a-findgen(k_lag)) ; ; normalize, equation (58) ; c=q/q(0) ; ; Set randomly the lag and makes sure that they are all different ; C_diff=dblarr(M,M) S=floor(k_lag*randomu(seed,M))+1 for i=0,M-1 do begin for j=i,M-1 do begin ;print,i,j C_diff(i,j)=abs(S(i)-S(j)) if (i ne j) then begin repeat begin nodiff=1 if (C_diff(i,j) eq 0.) then begin S(j)=floor(k_lag*randomu(seed))+1 C_diff(i,j)=abs(S(i)-S(j)) ;print,C_diff(i,j),i,j nodiff=0 endif endrep until nodiff endif endfor endfor ; print,S S_lag=S ; ; Set the random matrix ; C_rand=dblarr(M,M) for i=0,M-1 do begin for j=i,M-1 do begin C_rand(i,j)=c(C_diff(i,j)) C_rand(j,i)=C_rand(i,j) endfor endfor C_buffer=C_rand Print,'Eigenvalues...' ; ; Compute eigenvalues and eigenvectors ; evecs=1.d0 eigenvalues=eigenql(C_buffer,eigenvectors=evecs,/double) ; ; Set the number of eigenvalues ; r=rmax ; ; Set filter vector to 0 ; f1=dblarr(k_lag) ; ; Compute filter ; Print,'Filter...' for mm=0,r-1 do begin for i=0,M-1 do begin nn=S(i)-S index=where(nn ge 0.) f1(nn(index))=f1(nn(index))+evecs(i,mm)*evecs(index,mm)/M endfor endfor f1(0)=0. filters=f1 eigen_values=eigenvalues end