pro varadi_csa,filex,filey,rmax,M,K,filters,seed ;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; NAME: ; varadi ; CALLING SEQUENCE: ; varadi_csa,filex,filey,rmax,M,K,filters,seed ; PURPOSE: ; Compute the RLSSA from Varadi et al, ApJ, 1999, 526, 1052 ; INPUTS: ; filex filename for first series ; filey filename for second series ; 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 ; 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 24, 2000 ;------------------------------------------------------------------------------ ; ; Read files 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_x=readfits(filex,/silent) m_2=moment(a_x,/double) a_x=(a_x-m_2(0))/sqrt(m_2(1)) a_y=readfits(filey,/silent) m_2=moment(a_y,/double) a_y=(a_y-m_2(0))/sqrt(m_2(1)) !p.multi=[0,1,2] ; ; Check if x and y are the same, if yes flag=1 ; diff=total(a_x-a_y) flag=0 if (diff eq 0.) then flag=1 print,diff,flag ; a=randomn(seed,131072)*1.d0 ; ; N_a=N_elements(a_x) ; 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_x=reverse(a_x) a_rev_y=reverse(a_y) ; ; Pad the time series with zero ; a_new_x=dblarr(2*N_a) a_new_rev_x=dblarr(2*N_a) a_new_x(0:N_a-1)=a_x a_new_rev_x(0:N_a-1)=a_rev_x ; a_new_y=dblarr(2*N_a) a_new_rev_y=dblarr(2*N_a) a_new_y(0:N_a-1)=a_y a_new_rev_y(0:N_a-1)=a_rev_y ; ; Then compute the correlation using fft (3.3) ; Please note the use of the different normalization ; It was calibrated using the real formula (3.1) ; q=2.d0*N_a*fft(fft(a_new_x,-1,/double)*fft(a_new_rev_y,-1,/double),1,/double) q_x=2.d0*N_a*fft(fft(a_new_x,-1,/double)*fft(a_new_rev_x,-1,/double),1,/double) q_y=2.d0*N_a*fft(fft(a_new_y,-1,/double)*fft(a_new_rev_y,-1,/double),1,/double) ; ; Shift to put the max at 0 ; q=shift(double(q),-(N_a-1)) q_x=shift(double(q_x),-(N_a-1)) q_y=shift(double(q_y),-(N_a-1)) qt=dblarr(2L*k_lag-1) qt_x=dblarr(2L*k_lag-1) qt_y=dblarr(2L*k_lag-1) qt(0:k_lag-1)=q(0:k_lag-1) qt(k_lag:2L*k_lag-2)=q(2L*N_a-1-(k_lag-2):2L*N_a-1) ; qt_x(0:k_lag-1)=q_x(0:k_lag-1) qt_x(k_lag:2L*k_lag-2)=q_x(2L*N_a-1-(k_lag-2):2L*N_a-1) ; qt_y(0:k_lag-1)=q_y(0:k_lag-1) qt_y(k_lag:2L*k_lag-2)=q_y(2L*N_a-1-(k_lag-2):2L*N_a-1) ; ; 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 (3.5) corrected... ; norm=N_a-abs(findgen(2L*k_lag-1)-(k_lag-1)) qt=qt/double(norm) qt_x=qt_x/double(norm) qt_y=qt_y/double(norm) ; ; normalize, equation (58) ; if (flag eq 1) then begin c=qt/qt(0) c_x=qt_x/qt_x(0) c_y=qt_y/qt_y(0) endif else begin c=qt c_x=qt_x c_y=qt_y endelse ; ; Set randomly the lag and makes sure that they are all different ; C_diff=dblarr(M,M) C_diff_x=dblarr(M,M) C_diff_y=dblarr(M,M) S_x=floor(k_lag*randomu(seed,M))+1 S_y=floor(k_lag*randomu(seed,M))+1 if (flag eq 1) then S_y=S_x for i=0,M-1 do begin for j=i+1,M-1 do begin ;print,i,j C_diff_x(i,j)=S_x(i)-S_x(j) C_diff_y(i,j)=S_y(i)-S_y(j) repeat begin nodiff=1 if ((C_diff_y(i,j) eq 0.) or (C_diff_x(i,j) eq 0.)) then begin print,C_diff_x(i,j),C_diff_y(i,j),i,j,S_x(i),S_x(j),S_y(i),S_y(j) S_x(j)=floor(k_lag*randomu(seed))+1 S_y(j)=floor(k_lag*randomu(seed))+1 if (flag eq 1) then S_y(j)=S_x(j) C_diff_x(i,j)=S_x(i)-S_x(j) C_diff_y(i,j)=S_y(i)-S_y(j) print,C_diff_x(i,j),C_diff_y(i,j),C_diff(i,j),i,j nodiff=0 endif endrep until nodiff endfor endfor for i=0,M-1 do begin for j=0,M-1 do begin C_diff(i,j)=S_y(i)-S_x(j) endfor endfor ; ; Set the random matrix ; C_rand=dblarr(M,M) C_rand_x=dblarr(M,M) C_rand_y=dblarr(M,M) for i=0,M-1 do begin for j=i,M-1 do begin C_rand_x(i,j)=c_x((2L*k_lag-1)*(1-sign(C_diff_x(i,j)))*0.5+C_diff_x(i,j)) C_rand_y(i,j)=c_y((2L*k_lag-1)*(1-sign(C_diff_y(i,j)))*0.5+C_diff_y(i,j)) C_rand_x(j,i)=C_rand_x(i,j) C_rand_y(j,i)=C_rand_y(i,j) endfor endfor for i=0,M-1 do begin for j=0,M-1 do begin C_rand(i,j)=c((2L*k_lag-1)*(1-sign(C_diff(i,j)))*0.5+C_diff(i,j)) endfor endfor Print,'Eigenvalues...' ; ; Compute SVD according to Eq (2.3) ; svdc,C_rand,lambda,evecs_x,evecs_y,/double ; ; Compute sigma_x according to (2.5) ; pseudo_x=transpose(evecs_x) ## C_rand_x ## evecs_x lambda_x=dblarr(M) for i=0,M-1 do begin lambda_x(i)=pseudo_x(i,i) endfor ; ; Compute sigma_x according to (2.5) modified for y ; lambda_y=dblarr(M) pseudo_y=transpose(evecs_y) ## C_rand_y ## evecs_y for i=0,M-1 do begin lambda_y(i)=pseudo_y(i,i) endfor ; ; Set the number of eigenvalues (or highest correlation coefficient) ; r=rmax ; ; Compute B according to (3.9) ; ratio=lambda/lambda_x ; ; Compute correlation coefficient according to (3.10) ; r_corr=lambda/sqrt(lambda_x)/sqrt(lambda_y) ; ; Sort data ; if (flag eq 1) then begin index_c=reverse(sort(lambda)) endif else begin index_c=reverse(sort(r_corr)) endelse evecs=1.d0 eigenvalues=eigenql(C_rand_x,eigenvectors=evecs,/double) ; ; Set filter vector to 0 ; f1=dblarr(k_lag) ; ; Compute filter ; evecs_x=transpose(evecs_x) evecs_y=transpose(evecs_y) Print,'Filter...' for mm=0,r-1 do begin for i=0,M-1 do begin nn=S_x(i)-S_y index=where(nn ge 0.) if (index(0) ne -1) then begin f1(nn(index))=f1(nn(index))+evecs_x(i,index_c(mm))*evecs_y(index,index_c(mm))*ratio(index_c(mm))/M endif endfor endfor f1(0)=0. filters=f1 end