Non-linear curve fitting: y = Ae^(-t/tau) + C

Hi,

A signal I am trying to process in real time has multiple exponential
decay curves embedded in it, at known points in time. The curves can be
offset vertically by arbitrary amounts.

Each curve has the form

y = Ae^(-t/tau) + C    (Equation 1)

where A is the initial amplitude, t is time, tau is the time constant,
and C is some arbitrary vertical offset.

Fitting a curve in the form

y = Ae^(-t/tau)     (Equation 2)

is simple, because the equation can be linearized and you can use
linear least squares to find A and tau. Equation 2 has the vertical
offset C and the exponential term makes it impossible to isolate A, tau
and C when deriving the least squares relations.

So my question is, does anyone have any C or C++ (or any other
language) code snippets that can find the coefficients A, tau and C in
Equation 1, given a set of (t,y) data points?

Any help would be greatly appreciated.

Thanks,
Markus.

in article 1148481913.865727.64470@j55g2000cwa.googlegroups.com, Markus
Svilans at msvilans@gmail.com wrote on 05/24/2006 10:45:

> A signal I am trying to process in real time has multiple exponential
> decay curves embedded in it, at known points in time. The curves can be
> offset vertically by arbitrary amounts.
> 
> Each curve has the form
> 
> y = Ae^(-t/tau) + C    (Equation 1)
> 
> where A is the initial amplitude, t is time, tau is the time constant,
> and C is some arbitrary vertical offset.
> 
> Fitting a curve in the form
> 
> y = Ae^(-t/tau)     (Equation 2)
> 
> is simple, because the equation can be linearized and you can use
> linear least squares to find A and tau. Equation 2 has the vertical
> offset C and the exponential term makes it impossible to isolate A, tau
> and C when deriving the least squares relations.
> 
> So my question is, does anyone have any C or C++ (or any other
> language) code snippets that can find the coefficients A, tau and C in
> Equation 1, given a set of (t,y) data points?

personally, it would seem to me that using least square fitting to get a
single A and C is easy (assuming tau is known) and getting tau (or 1/tau) is
hard.  if there are multiple overlapping components of Eq. 1 and you want
all of the A's and tau's and C's, i guess i would say that there is a single
C.

in any case, you take the partial derivative of the square error with
respect to each unknown parameter and set each of those to zero.  you get an
equal number of equations as unknown parameter (that, if solved, yields a
unique solution), but solving for some of the unknowns might be a bitch.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

robert bristow-johnson wrote:
> in article 1148481913.865727.64470@j55g2000cwa.googlegroups.com, Markus
> Svilans at msvilans@gmail.com wrote on 05/24/2006 10:45:
>
> > A signal I am trying to process in real time has multiple exponential
> > decay curves embedded in it, at known points in time. The curves can be
> > offset vertically by arbitrary amounts.
> >
> > Each curve has the form
> >
> > y = Ae^(-t/tau) + C    (Equation 1)
> >
> > where A is the initial amplitude, t is time, tau is the time constant,
> > and C is some arbitrary vertical offset.
> >
> > Fitting a curve in the form
> >
> > y = Ae^(-t/tau)     (Equation 2)
> >
> > is simple, because the equation can be linearized and you can use
> > linear least squares to find A and tau. Equation 2 has the vertical
> > offset C and the exponential term makes it impossible to isolate A, tau
> > and C when deriving the least squares relations.
> >
> > So my question is, does anyone have any C or C++ (or any other
> > language) code snippets that can find the coefficients A, tau and C in
> > Equation 1, given a set of (t,y) data points?
>
> personally, it would seem to me that using least square fitting to get a
> single A and C is easy (assuming tau is known) and getting tau (or 1/tau) is
> hard.  if there are multiple overlapping components of Eq. 1 and you want
> all of the A's and tau's and C's, i guess i would say that there is a single
> C.

All of A, tau and C are unknown.

I should clarify: I would like to be able to find the unknowns using
some sort of non-linear curve fitter, for a single decay curve in the
form of Equation 1. Once I know how to do that, I have a strategy in
mind for subtracting each decay curve from subsequent decay curves.

> in any case, you take the partial derivative of the square error with
> respect to each unknown parameter and set each of those to zero.  you get an
> equal number of equations as unknown parameter (that, if solved, yields a
> unique solution), but solving for some of the unknowns might be a bitch.


It is exactly the bitchy solving for unknowns that is my problem. It is
impossible to get least squares equations for Equation 1, making it
necessary to resort to an iterative non-linear curve fitter.

Thanks,
Markus.

Markus Svilans wrote:
> It is exactly the bitchy solving for unknowns that is my problem. It is
> impossible to get least squares equations for Equation 1, making it
> necessary to resort to an iterative non-linear curve fitter.
>
> Thanks,
> Markus.

Hello Markus,

What you want to use is the method of differential correction. LeGendre
used it to adjust planetary orbital parameters given new observations.

Take your original equation and perform implicit differentiation on it
with respect to one of your parameters. Then cross multiply by the
denominator portion of the differential.

You should get


dy = [exp(-t/tau)] da  + [ (A*t/tau)exp(-t/tau)] d tau  + dc


Think of this in matrix form


[dy] = [ exp(-t/tau)      (A*t/tau)exp(-t/tau)      1] [ da   d tau
dc ]transposed


The left hand side is a column vector (n terms deep) of your true
measured values minus your predicted values. The matrix on the right is
formed as a n by 3 matrix where each row has its own t value. And the
column vector on the right (just 3 terms) is what you are solving for.
So you start with seed values for A, tau, and C. Use these to predict y
values by using your model equation. Find the dy by finding the deltas
between the measured and predicted data. And then back solve for  da, d
tau, and dc. Since your system is likely (hopefully) over determined
(i.e., n>3), use a least squares method to find these. Then update your
A, tau, and C with the changes implied by da, d tau, and dC. Iterate
this last part using the new A, tau, and C as seed values until nothing
changes.

IHTH,

Clay

Markus Svilans wrote:
> Hi,
>
> A signal I am trying to process in real time has multiple exponential
> decay curves embedded in it, at known points in time. The curves can be
> offset vertically by arbitrary amounts.
>
> Each curve has the form
>
> y = Ae^(-t/tau) + C    (Equation 1)
>
> where A is the initial amplitude, t is time, tau is the time constant,
> and C is some arbitrary vertical offset.
>
> Fitting a curve in the form
>
> y = Ae^(-t/tau)     (Equation 2)
>
> is simple, because the equation can be linearized and you can use
> linear least squares to find A and tau. Equation 2 has the vertical
> offset C and the exponential term makes it impossible to isolate A, tau
> and C when deriving the least squares relations.
>
> So my question is, does anyone have any C or C++ (or any other
> language) code snippets that can find the coefficients A, tau and C in
> Equation 1, given a set of (t,y) data points?
>
> Any help would be greatly appreciated.
>
> Thanks,
> Markus.
This can be solved in an number of ways.
Excel has a LINEST function that can do non linear fitting.   For
sonething as simple as what you are doing this should be enough and you
probably already have Excel.  I think the LINEST function is an add in.
SciLab has a Levenberg-Marquardt function that can find the minimum
error between two functions or a function and some data.
Mathcad has a Minerr function.
Numerical Recipies in C has a Powell, DFP, and Levenberg-Marquardt
programs that can minimize or optimize errors.
One can also search for levmar and lmfit and find other C or C++
sources on the internet.

Once you get into non-linear optimization a whole world opens up.

Peter Nachtwey

Clay wrote:
> Markus Svilans wrote:
> > It is exactly the bitchy solving for unknowns that is my problem. It is
> > impossible to get least squares equations for Equation 1, making it
> > necessary to resort to an iterative non-linear curve fitter.
> >
> > Thanks,
> > Markus.
>
> Hello Markus,
>
> What you want to use is the method of differential correction. LeGendre
> used it to adjust planetary orbital parameters given new observations.
>
> Take your original equation and perform implicit differentiation on it
> with respect to one of your parameters. Then cross multiply by the
> denominator portion of the differential.
>
> You should get
>
>
> dy = [exp(-t/tau)] da  + [ (A*t/tau)exp(-t/tau)] d tau  + dc
>
>
> Think of this in matrix form

Hi Clay,

Thank you for your insight and clear explanation of what to do. This is
exactly the kind of solution I was looking for. The next thing is to
code it up and try it out.

I just found one minor correction to make in your math:

>
> [dy] = [ exp(-t/tau)      (A*t/tau)exp(-t/tau)      1] [ da   d tau
> dc ]transposed
>

The partial derivative of with respect to tau should be:

(A*t)/(tau^2) exp(-t/tau) d tau

The tau should be squared in the denominator out front. :)

Regards,
Markus.

ya know what, we were talking about something very similar:

 http://groups.google.com/group/comp.dsp/msg/381d1a02136b0346

and then it occured to me that a result from, believe or not- my
master's thesis in 1981, might be relevant.

restating it a little:

      y(t) = A*e^(-t/tau) + C

anyway the above is the solution to the following diff eq:

     y(t) - C = A*e^(-t/tau)

     y'(t) = (-A/tau)*e^(-t/tau)

     y'(t) = (C - y(t))/tau

defining the error function to be:

     e(t) = y'(t) - (C - y(t))/tau

            t
    E = integral{ e(u) du} = integral{ y'(u) - (C - y(u))/tau du}
           t0

taking the partial derivative of E w.r.t. C and 1/tau and setting both
to zero gets you (from my memory):

             integral{y(u)^2 du} - integral{ y(u) du}*(y(t)+y(t0))/2
      C  =  -----------------------------------------------------------
                   integral{ y(u) du} - (y(t)+y(t0))/2

the integrals all taken over the interval, t0 < u < t

i think you can figure out what 1/tau is and then maybe A.  it looks
like that somehow, A comes out in the wash.  you should be able to use
this estimate for C to take it out of the samples in the original
equation:

    y(t) - C = A*e^(-t/tau)

take the natural logarithm:

    log( y(t) - C ) = (-1/tau)*t + log(A)

and fit to a straight line over the interval, t0 < u < t, and determine
a best estimate for -1/tau and log(A).

i think that might be enough stuff to do it, but i want to leave the
labor to someone else.

r b-j

looks like the OP is with some big Canadian company:

http://www.aeroquestsurveys.com/english/team.php

i hope i didn't give that one away too cheap.

r b-j

made one mistake.  it should be:


             integral{y(u)^2 du} - integral{ y(u) du}*(y(t)+y(t0))/2
      C  =  -----------------------------------------------------------
                  integral{ y(u) du} - (t - t0)*(y(t)+y(t0))/2


i think that is better dimensionally, also (than the previous expression).

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

in article C09A9B0B.14D07%rbj@audioimagination.com, robert bristow-johnson
at rbj@audioimagination.com wrote on 05/24/2006 23:21:

> 
> made one mistake.  it should be:
> 
> 
>              integral{ y(u)^2 du} - integral{ y(u) du}*(y(t)+y(t0))/2
>       C  =  -----------------------------------------------------------
>                   integral{ y(u) du} - (t-t0)*(y(t)+y(t0))/2
> 

(all integrals over the interval:  t0 < u < t .

made another mistake in presentation:

the error function is:

    e(t) = y'(t) - (C - y(t))/tau


but the total square error to minimize is:

             t                      t
    E = integral{ e(u)^2 du} = integral{ (y'(u) - (C - y(u))/tau)^2 du}
             t0                     t0

taking derivatives of E with respect to C and 1/tau and setting to zero,
results in the estimate for C above.


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."