# 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.

• so by using your DC offset removal filter, does that mean that $C$ is zero? – robert bristow-johnson Dec 26 '18 at 20:37
• The filter is just to measure the integral of the pulse, and one of its property is to remove the offset – Dams0622 Dec 26 '18 at 22:19
• my question remains. – robert bristow-johnson Dec 26 '18 at 22:34
• and may i suggest modifying your equation to be completely discrete-time? like this: $$x[n]=Ae^{-n/\tau}u[n]+C$$ and use small-case $x[\cdot]$ for functions of time and leave the capital $X[\cdot]$ for frequency domain. – robert bristow-johnson Dec 26 '18 at 22:37
• looks like Dawg pointed you to another answer. i am still "on the road". see if that other answer helps, but if you want one that takes many more than 4 samples and computes $\tau$ using a least-square fit, then you might need something else. but Dawg's answer is a good one, i think. – robert bristow-johnson Dec 27 '18 at 2:40

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

#---- 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()


• Thank you Dawg. Your answer and the reference post are a lot to digest for me, but I think that get the point. I am going to test the second part and take the measure with all the samples. The noise is a real problem, especially for small pulses. I will provide the results on real pulses as soon as I get back to my laptop – Dams0622 Dec 27 '18 at 11:43
• @Dams0622,You're welcome. You will get better results at the steeper part of your curves. This is because the differences you will be taking will be larger compared to your noise level. When the curve gets close to flat the difference gets close to zero and value becomes almost all noise, a very poor SNR. The nice thing about this solution is that it is fairly efficient. It would be nice to see some results. – Cedron Dawg Dec 27 '18 at 13:54
• @Dams0622, I have added some sample code to show that the formulas work in the noiseless case. As you can see, the calculations are pretty straightforward. You may get a better read by finding the best fit $A$ and $C$ then casting away the outliers and recalculating. – Cedron Dawg Dec 28 '18 at 5:45
• I reproduced your script in Excel (my only tool available at the moment) and it is working really well. I will be able to produce results with real pulses in a few days – Dams0622 Dec 28 '18 at 11:43