Rapid fitting of exponentials

 

//  The following function finds A, Tau, and B from an exponential decay by
//	(1) taking the FFT, and 
//	(2) using the complicated function for alpha found in Rev. Sci. Instrum. 79, 023108 (2008)
Function FFTFit(RDTWave, A, Tau, B)
	wave RDTWave
	Variable &A, &Tau, &B
	
	FFT/OUT=1/DEST=RDTWave_FFT RDTWave
	variable tm = pnt2x(RDTWave, (numpnts(RDTWave)-1)) - pnt2x(RDTWave, 0)
	variable omega = pnt2x(RDTWave_FFT,1)*2*pi
	variable alpha = ( 1/ deltax(RDTWave) )* ln(cos(2*pi/numpnts(RDTWave))+(real(RDTWave_FFT[1])/imag(RDTWave_FFT[1]))*sin(2*pi/numpnts(RDTWave) )  ) 
	Tau = 1/alpha
	A = (alpha^2 + omega^2)*(Real(RDTWave_FFT[1])*deltax(RDTWave))/( alpha*(1 - exp(-alpha*tm))  )
	B = (RDTWave_FFT[0]/numpnts(RDTWave) ) - ( A*(1 - exp(-alpha*tm) )/(alpha*tm) )
	
end




//  The following function finds A, Tau, and B from an exponential decay by
//	the method of corrected successive integration 
//	Rev. Sci. Instrum. 79, 023108 (2008)
//
Function DSIFit(RDTWave, A, Tau, B)
	wave RDTWave
	Variable &A, &Tau, &B
	
	variable Status = 0
	make /O/D/N=3 CoefMat=nan
	make /O/D/N=3 DataMat=nan
	make/O/D/N=(3,3) ArrayMat=nan
	make /O/D/N=(numpnts(RDTWave)) twave 
	CopyScales RDTWave twave
	twave = x
	// S_I
	Integrate   RDTWave /D=RDT_Int
	variable /D S_I = sum(RDT_Int)
	// S_II
	variable /D S_II = MatrixDot(RDT_Int, RDT_Int)
	// S_f
	variable /D S_f =sum(RDTWave)
	// S_tI
	variable /D S_tI = MatrixDot(twave, RDT_Int)
	// S_fI
	variable /D S_fI = MatrixDot(RDTWave, RDT_Int)
	// S_ft
	variable /D S_ft = MatrixDot(RDTWave, twave)
	// S_t
	variable /D S_t = deltax(RDTWave)*(numpnts(RDTWave)-1)*(numpnts(RDTWave))/2
	// S_tt
	variable /D S_tt = (deltax(RDTWave)^2)*(numpnts(RDTWave)-1)*(numpnts(RDTWave))*(2*(numpnts(RDTWave)-1)+1)/6
	
	ArrayMat = {{numpnts(RDTWave), S_I, S_t},{S_I, S_II, S_tI},{S_t, S_tI, S_tt}}
	DataMat={S_f, S_fI, S_ft}
 
	MatrixOp /O CoefMat =Inv(ArrayMat) x  DataMat 
	Tau = Deltax(RDTWave)/ln(1 - coefmat[1]*Deltax(RDTWave))  //  Rectangular Integration
	B = -CoefMat[2]/CoefMat[1]
	A =  exp(-Deltax(RDTWave)/Tau)*CoefMat[0] - B
	
	Return Status
end

Cool, Mike.

What are A and B? Amplitude and vertical offset? Oh- a little experimentation shows that's correct.

I tried it on this:

Make/D junk
setscale/I x 0,1,junk
junk += gnoise(.1) // pretty noisy

•printFFTFit(junk)
A= 1.37807 Tau= 0.0749016 B= -0.0283905
•printDSIFit(junk)
A= 0.984629 Tau= 0.102131 B= -0.0311092

Using Igor's built-in fit I got
y0 =-0.030051 ± 0.0128
A =0.99762 ± 0.0575
tau =0.10196 ± 0.00998

The FFT fit seems to be insensitive to where the exponential starts (that is, setscale/I x 1,2,junk gives the same answer), but DSIFit gives quite different answers.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

April 3, 2008 at 12:30 pm - Permalink

Forum

Support

Gallery

Igor Pro 10

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More