pro test_varadi_csa,amplitude
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; NAME:
;	test_varadi_csa
; CALLING SEQUENCE:
;	test_varadi_csa,amplitude
; PURPOSE:
;	Test for detecting a pure sinewave embedded into noise with 2 different signal
; INPUTS:
;	amplitude: mean amplitude of the sine waves (sqrt(2.) or less is a good number)
; OPTIONAL KEYWORDS:
; 	none
; OUTPUTS:
;	none
; OPTIONAL OUTPUT:
;	none
; LIMITATIONS:
;	None yet
; COMMONS:
;	None
; PROCEDURES USED:
;       fft,writefits
; MODIFICATION HISTORY:
;	Beta 1: WRITTEN, Thierry Appourchaux, February 6, 2001
;------------------------------------------------------------------------------
	M=250
	K=25000
	seed=systime(1)
	filters=dblarr(K)
	filters_plot=dblarr(K)
	values=dblarr(M)
	eigen_val1=dblarr(M)
	values1=dblarr(M)
	eigen_val=dblarr(M)
	b=fltarr(131072)
	N_a=N_elements(b)
;
; Put 11 set of sine waves
;
	freq1=findgen(11)*500.+800.
;
; Create sum of sine waves
;
	for kk=0,10 do begin
		b=b+0.01*amplitude*sqrt(2.)*cos(2.*!pi*findgen(N_a)/(1./(1.e-06*freq1(kk))/60.))
	endfor
;
; Plot fft of sine waves, initial noise and noise + waves
;
	!p.multi=[0,1,2]
	a=randomn(seed,131072)
	c=randomn(seed,131072)
	writefits,'test_noise',a
	writefits,'test_signal1',a+b
	writefits,'test_signal2',c+b
	window,0,title='Noise and signal'
	wset,0
	plot,abs(fft(a,-1)),xrange=[0.,131072/2.],charsize=2.0
	plot,abs(fft(b,-1)),xrange=[0.,131072/2.],charsize=2.0
	window,2,title='Noisy data'
	wset,2
	!p.multi=[0,1,2]
	plot,abs(fft(a+b,-1)),xrange=[0.,131072/2.],charsize=2.0
	plot,abs(fft(c+b,-1)),xrange=[0.,131072/2.],charsize=2.0
;
; Define the filtering arrays
;
	filter_new=dblarr(N_a)
	filter_shift=dblarr(N_a)
;
; Create array for lags
;
	S_lag=fltarr(1000)
	window,1,title='Reconstructed filter'
	wset,1
	for i=0,99 do begin
		print,'User reassuring index=',i
		varadi_csa,'test_signal1','test_signal2',10,M,K,filters,seed
		values1=values1+eigen_val1
		filters_plot=filters_plot+filters
		filter_new(*)=0.
		filter_new(0:K-1)=filters_plot/(i+1)
		filter_shift(0:K-2)=filters_plot(1:K-1)/(i+1)
		filter_new=filter_new+reverse(filter_shift)
		;plot,filter_new
		uu=abs(fft(filter_new,-1))^2
		plot,uu/max(uu),xrange=[0.,131072/2.],charsize=2.0,title=strtrim(i,1)
		;oplot,(uu_b+0.5)/1.5,color=200
	endfor
end