Rapid fitting of exponentials
meverest
// 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
// (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
// 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
Forum
Support
Gallery
Igor Pro 9
Learn More
Igor XOP Toolkit
Learn More
Igor NIDAQ Tools MX
Learn More
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