# 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

Forum

Support

Gallery