Rapid fitting of exponentials

Tue, 04/01/2008 - 04:09 pm

//  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 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More