pro fftpacket,series1,fftseries,smoother,apodize,min,max,relative=relative,broadcast=broadcast,spikes=spikes,nofft=nofft,zeroer=zeroer
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; NAME:
;	FFTPACKET 
; CALLING SEQUENCE:
;	fftpacket,series,fftseries,smoother,apodize,min,max,/relative,/broadcast,spikes=value
; PURPOSE:
;	Fourier transform of a time series 
; INPUTS:
;	series 		time series to be treated (number of points is a power of 2) 
;	smoother	Perform smoothing.  It is a dblarr(3)
;		0	No smoothing [0,*,*]
;		1	Binomial smoothing [1, sampling time in sec, cut-on in mHz]
;		2	Polynomial smoothing [2, degree of polynome, *]
;		3	Triangle smoothing [3,Width in minutes,*]
;		4	Spline interpolation
;	apodize		Apodization function
;		0	No apodization (default)
;		1	Welsh (parabolic)
;		2	Hanning (cosine)
;		3	Parzen (triangle)
;	min		min value for invalid data
;	max		max value for invalid data
; OPTIONAL KEYWORDS:
;	relative	after smoothing the relative difference is computed
;	broadcast	plot a few plots to see what is going on
;	nofft		do not compute fft. fftseries is replaced by the series before fft
; OUTPUTS:
;	fftseries 	power spectra (Npoints)
; OPTIONAL OUTPUT:
;	none
; EXAMPLE:
;	To be written
; LIMITATIONS:
;	None ?
; COMMONS:
;	None
; PROCEDURES USED:
;       fft, convol, poly_fit, moment, total
; MODIFICATION HISTORY:
;	Beta 1: WRITTEN, Thierry Appourchaux, September, 25, 1995.
;	Beta 2: Various improvements still to be checked...
;	Beta 3:	Everything is done in dbl precision, May 6, 1996.
;	Beta 4: Add smoothing by linear interpolation, May 7, 1996.
;	Beta 5: With /nofft we return the apodized and normalized time series
;		May 7, 1996
;	Beta 6: Removed a bug preventing to smooth with a polynomial where there
;		are no invalid data, May 9, 1996.
;	Beta 7: Now compute the number of zeros put by spikes, and include the 
;		invalid data inot the normalization factor, August 19, 1996
;------------------------------------------------------------------------------
;
;
;
	Npoints=N_elements(series1)
;
	if (keyword_set(effect)) then begin
		series1(0:Npoints-1)=0.d0
		series1(Npoints/2)=1.d0
	endif
	series=series1*1.d0
;
; Here we create an index of invalid data
;
	invalid=where((series LT min) OR (series GT max))
;
; Here we create an index of valid data
;
	valid=where((series GT min) AND (series LT max))
	if (valid(0) EQ -1) then begin
		print,'WHAT! No valid data?'
		stop
	endif
;
; We replace the invalid data by a linear interpolation of the valid data. This is for
; reducing edge effects around the invalid data.
;	
	


	if (invalid(0) NE -1) then begin
		bounce=dblarr(N_elements(valid)/1440.)
		Nbounce=N_elements(bounce)
		bounce=series(valid(findgen(Nbounce)*1440))
		u=dindgen(Npoints)
		x=valid(findgen(Nbounce)*1440)
		smoothed4=interpol(bounce,x,u)
		series(invalid)=smoothed4(invalid)
	endif


;	if (invalid(0) NE -1) then begin
;		meanout=moment(series(valid))
;		series(invalid)=meanout(0)
;	endif
;
;
;
;
	If (smoother(0) NE 0) then begin
;
; ********************* Start binomial smoothing ******************************
;
	If (smoother(0) EQ 1) then begin
		print,'data are being smoothed with the binomial stuff...'
;
;
;	plot,series
;	read,blurp
;
;
; define binomial parameters for smoothing
;
;
	sampling=smoother(1)
	cuton=smoother(2)
	fre=cuton*sampling/1000.
;
	result=[6.3447729d0,-108.40157d0,866.92492d0,-2966.7498d0]
;
; Compute number of pass for the binomial smoothing
;
	npass=floor(exp(result(0)+result(1)*fre+result(2)*fre^2+result(3)*fre^3))+1

;     	npass=4*npass
;
	print,'npass=',npass
;
	binomial=dblarr(5)
	binomial(0)=0.d0
	binomial(1)=0.25d0
	binomial(2)=0.5d0
	binomial(3)=0.25d0
	binomial(4)=0.d0
;
	bounce=dblarr(3*Npoints)
;
; mirror beginning
	bounce(0:Npoints-1)=reverse(series(0:Npoints-1))
; real array
	bounce(Npoints:2*Npoints-1)=series(0:Npoints-1)
; mirror end
	bounce(2*Npoints:3*Npoints-1)=reverse(series(0:Npoints-1))
;
; smooth the date using CONVOL and the binomial smoothing
	smoothed=series
	for i=0,npass do begin
			bufferi=bounce
			buffero=convol(bufferi,binomial)
			bounce=buffero
	endfor
;
; put the bounced array in smoothed
	smoothed(0:Npoints-1)=bounce(Npoints:2*Npoints-1)
	endif
;
; ********************* End binomial smoothing ******************************
;
;
; ****************** Start polynomial smoothing ******************************
;

	If (smoother(0) EQ 2) then begin
		print,'data are being smoothed with the polynomial...'
		deg=smoother(1)
		afit=poly_fit(valid*1.d0,series(valid),deg)
			
		smoothed=dblarr(Npoints)
		smoothed(0:Npoints-1)=afit(deg)
		for i=1,deg do begin
			j=deg-i
			smoothed=afit(j)+dindgen(Npoints)*smoothed
		endfor
	endif
;
; ****************** End polynomial smoothing *********************************
;

;
; ********************* Start triangle smoothing ******************************
;
	If (smoother(0) EQ 3) then begin
		print,'data are being smoothed with the triangle stuff...'
;
;
;	plot,series
;	read,blurp
;
; define boxcar parameters for smoothing
;
;
	width=floor(smoother(1))
;*1440
;
;
;
	bounce=dblarr(3*Npoints)
;
; mirror beginning
	bounce(0:Npoints-1)=reverse(series(0:Npoints-1))
; real array
	bounce(Npoints:2*Npoints-1)=series(0:Npoints-1)
; mirror end
	bounce(2*Npoints:3*Npoints-1)=reverse(series(0:Npoints-1))
;
; smooth the date using SMOOTH and the binomial smoothing
	smoothed=series

	buffero=smooth(bounce,width)	
	bounce=buffero

	buffero=smooth(bounce,width)
	bounce=buffero

;
; put the bounced array in smoothed
	smoothed(0:Npoints-1)=bounce(Npoints:2*Npoints-1)
	endif
;
; ********************* End triangle smoothing ******************************
;

;
; ********************* Start interpolation smoothing ******************************
;
	If (smoother(0) EQ 4) then begin
		print,'data are being smoothed with the interpolation stuff...'
;
;
;
;
	
;
; Do not redo the interpolation if there are invalid data
;
	if (invalid(0) EQ -1) then begin
		bounce=dblarr(N_elements(valid)/1440.)
		Nbounce=N_elements(bounce)
		bounce=series(valid(findgen(Nbounce)*1440))
		u=dindgen(Npoints)
		if (Nbounce eq 1) then begin
			Print,'Sorry we cannot perform smoothing by interpolation if you do not hav emore than 1440 points'
			stop
		endif
		x=valid(findgen(Nbounce)*1440)
		smoothed4=interpol(bounce,x,u)
	endif



	smoothed=smoothed4
	

	endif
;
; ********************* End triangle smoothing ******************************
;




;
;
; ************* Start relative or absolute variations *************************
;
; compute relative variation or absolute variation
	If (not keyword_set(relative)) then begin
		print,'Computing absolute variations...'
		out=series-smoothed
	endif else begin
		print,'Computing relative variations...'
		out=(series/smoothed-1.)
	endelse
	If (keyword_set(effect)) then begin
		out=smoothed
	endif
;
; ************* End relative or absolute variations ****************************
;

;
;
; ********************** No smoothing: raw data! *******************************
;
;
	endif else begin
;
	If (smoother(0) EQ 0) then begin
		print,'***WARNING***'
		print,'You chose NOT to smooth.  We compute the fft of the raw time series'
		if (invalid(0) NE -1) then begin
			meanout=moment(series(valid))
			series(invalid)=meanout(0)
		endif
		out=series
	endif else begin
		print,'This is not a valid number for smoother(0)=',smoother(0)
		stop
	endelse
	endelse
;
; The invalid data are put to ZERO (or the mean) whatever the treatment is (Smoothing or not)
;
	if (invalid(0) NE -1) then begin
		meanout=moment(out(valid(*)))
		out(invalid(*))=0.d0
	endif
	If (smoother(0) EQ 0) then begin
		if (invalid(0) NE -1) then begin
			meanout=moment(out(valid(*)))
			out(invalid(*))=meanout(0)
		endif
	endif

;
	if (keyword_set(broadcast)) then begin
		plot,out,title='output data after smoothing'
		read,blurp
	endif
;
; ******************* End of smoothing or no smoothing *************************
;
;
;
	If (smoother(0) NE 0) then begin
	if (keyword_set(spikes)) then begin

		zeroer=out
	
		bufferm=dblarr(100)
	
		sigmamax=N_elements(valid)/100-1
		sigmaarr=dblarr(sigmamax)
	
	
		for k=0,sigmamax-1 do begin
			bufferm(0:99)=out(valid(long(k)*100:99+long(k)*100))
			sigmaarr(k)=stdev(bufferm)
		endfor
;		plot_io,sigmaarr
		histoplot,alog(sigmaarr)
		read,blurp
;		sig=min(sigmaarr(where(sigmaarr GT 0.)))
		
		new=alog(sigmaarr(where(sigmaarr GT 0.d0)))

		Nnew=N_elements(new)

		Nn1=(Nnew+Nnew mod 2)/2


		
		sig1=exp(total(new(0:Nn1-1))/Nn1)

		sig2=exp(total(new(Nn1:Nnew-1))/(Nnew-Nn1))

		print,'Sigma mean1=',sig1,' ,sigma mean2=',sig2
		sig=sig1 < sig2
		print,sig
		zeroer=floor(abs(out/spikes/sig)) LT 1
		if (keyword_set(broadcast)) then begin
			plot,zeroer,title='zeroer',yrange=[0,2]
			read,blurp
		endif
		
		print,'Number of zeros=',N_elements(zeroer)-total(zeroer)
		out=out*zeroer
	endif
	endif






	apod=dblarr(Npoints)
	apod(0:Npoints-1)=1.d0
;
	if (apodize NE 0) then begin
		if (apodize EQ 1) then begin  ;welsh
			print,'Apdodization is Welsh'
			apod=1.d0-(2.d0*(dindgen(Npoints)-(Npoints-1.d0)/2.d0)/(Npoints+1.))^2
		endif
		if (apodize EQ 2) then begin  ;hanning
			print,'Apdodization is Hanning'
			apod=hanning(Npoints)
		endif
		if (apodize EQ 3) then begin  ;parzen
			print,'Apdodization is Parzen'
			apod=1.d0-abs((2.d0*(dindgen(Npoints)-(Npoints-1.d0)/2.d0)/(Npoints+1.)))
		endif
	endif
;
	if (keyword_set(broadcast)) then begin
		plot,apod,title='apodization function'
		read,blurp
	endif
	out=apod*out
;
; Now we compute the normalization factor from the apodized window function
; including the gaps
;
	window=dblarr(Npoints)
	window(1:Npoints-1)=1.d0
;
	if (invalid(0) NE -1) then begin
		window(invalid)=0.
	endif

	if (keyword_set(spikes)) then begin
		window=window*zeroer
		print,'Total number of zeros=',N_elements(window)-total(window)
	endif

;
	window=window*apod
;
	Norm=sqrt(Npoints/total(window^2))	

	print,'Norm=',Norm


;

;
; Now remove mean and renormalize
;
	

		meanout=moment(out)
		print,meanout

		print,total(out)/N_elements(out)
;
		out=out-meanout(0)
		print,meanout(0),total(out)/N_elements(out)
;
		out=Norm*out
		
	if (keyword_set(nofft)) then begin
		fftseries=out
		if (keyword_set(broadcast)) then begin
			plot,out,title='time series to be processed'
			read,blurp
		endif
	endif else begin

		
		if (keyword_set(broadcast)) then begin
			plot,out,title='time series to be processed'
			read,blurp
		endif

		series=out
		fftseries=abs(fft(out,-1))^2
;		if (keyword_set(effect)) then begin
;			fftseries=(1.-sqrt(fftseries/max(fftseries)))^2
;		endif
	endelse
;
end