filters - Compute time constant of discrete exponential decay pulses

I want to retrieve (within a FPGA) the time constant $\tau$ of incoming pulses of the form $$x[n]=Ae^{-n/\tau}u[n]+C$$ In addition to the offset, inputs are subject to noise.

To measure the integral of the pulses, I used a trapezoidal filter (2 successive differentiator+accumulator) whose properties allow for DC offset removal and noise attenuation. Starting from that, I was wondering if there is were a smart way to compute an approximation of the time constant.

A few articles mention Fourier Transform or Corrected Successive Integration technics but I could not find any free material.

Any tip to put me in the right direction would be appreciated.

Answer

Take four sample values.

$$ y_1 = x[p] $$ $$ y_2 = x[q] $$ $$ y_3 = x[p+d] $$ $$ y_4 = x[q+d] $$

Subtract the third from the first. Notice that the $C$s cancel out.

$$ y_1 - y_3 = x[p] - x[p+d] = Ae^{-p/\tau}-Ae^{-(p+d)/\tau} = Ae^{-p/\tau} \left( 1 - e^{-d/\tau} \right) $$

Similarly, subtract the fourth from the second.

$$ y_2 - y_4 = Ae^{-q/\tau} \left( 1 - e^{-d/\tau} \right) $$

Take the quotient of these two differences. Notice that the $A$s cancel out.

$$ \frac{y_1 - y_3}{y_2 - y_4} = e^{(q-p)/\tau} $$

Take the log of each side.

$$ \ln \left( \frac{y_1 - y_3}{y_2 - y_4} \right) = \frac{q-p}{\tau} $$

Solve for $\tau$.

$$ \tau = \frac{q-p}{ \ln \left( \frac{y_1 - y_3}{y_2 - y_4} \right) } $$

I did this more thoroughly in another answer, but this should give you start. You want the differences to be large compared to the noise level, meaning d should be large. $(q-p)$ being larger will help too. I am guessing that four evenly spaced samples are as good as anything. Something to investigate.

If noise is a problem you can spread the calculation across more points. You can either do an average of neighbor points for each $y$ value as long as you use the same weighting or employ more points directly in the quotient as in my referenced answer.

The other answer will also show you how to find the best fit values of $A$ and $C$.

---

Using more points to mitigate noise.

Here is a way you can fully utilize all your samples to reduce the effects of noise I believe as much as possible.

Select the number of samples ($N$) to be a multiple of four. Partition the samples into four quarters. Then set the $y$ values to be the sum of the samples of the corresponding quarters. You don't need to take an average (an extra division) because that washes out in the taking of the quotient.

$$ q-p = N/4 $$

In addition to taking one measure of the entire curve, you might want to try a sliding window taking many smaller measures. The values you get should be fairly constant. If not, your assumption of exponential decay may not be correct.

---

Sample Python code:

import numpy as np

#================================================
def main():

#---- Set Parameters


        A   = 1.234
        C   = 5.678
        tau = 24.08016032

        N = 100

#---- Construct the pulse

        x = np.zeros( N )


        for n in range( N ):
            x[n] = A * np.exp( -n / tau ) + C

        print x

#---- [Add noise here]

#---- Example Using a Single Point

        p = 10

        q = 20
        d = 20

        y1 = x[p]
        y2 = x[q]
        y3 = x[p+d]
        y4 = x[q+d]

        S = ( y1 - y3 ) / ( y2 - y4 )


        #++[Check S here]

        tau_calc1 = ( q - p ) / np.log( S )

        print tau_calc1

#---- Example by Quarters

        y1 = np.sum( x[p:p+10] )
        y2 = np.sum( x[q:q+10] )

        y3 = np.sum( x[p+d:p+d+10] )
        y4 = np.sum( x[q+d:q+d+10] )

        S = ( y1 - y3 ) / ( y2 - y4 )

        #++[Check S here]

        tau_calc2 = ( q - p ) / np.log( S )

        print tau_calc2


#================================================
main()