Dear Luke,

Maybe fit of overlapping Hermittian polynomials to data (by pchip() in matlab) could allow you to compute continuous 1st order derivatives and detect acute rates of change.

Acute changes of the signal could also be detected by histogram statistics (matlab uses Otsu's method to perform such a segmentation).

But if you need a global analytic expression for the whole signal (not local overlapping ones) and the formula of a prototype is not available, I think, in general, you'll have to project the signal onto a kernel of functions, e.g. wavelets, gaussian mixtures, etc.

In case that you have a formula for the prototype that fits the signal, you could detect the optimal parameters of this model via matlab fit() function. The explicit formula of the parametric model should be predefined via 'fittype'.

Best regards.

If you want to remove noise from your data, you use smooth() first and then fit.

You'll need to be careful with noise removal, as the rise and fall information may well have frequency components in the noise. Differentiation will also amplify noise. Sometimes a short (3 or 5 point) median filter is best if the noise is visible as discrete spikes in the data stream.

I would suggest the use of an AR model. Then look at the poles. You are looking for the most important real pole (inside the unit circle, of course). Convert this pole from z to s using the bilinear or the Euler s = (1-z^(-1))/T transformation. Then the time constant is the inverse of the pole in plane s.

Hi,

there is a very useful and self-explanatory curve-fitting tool (provided you have the toolbox installed). Just type "cftool" in command window of MATLAB. Have the data stored in a variable then it'll appear in the tool's option. You can choose from a variety of fitting options, and the fit implicitly takes care of noise. Hope this helps!

Hi Luke, You need Curve Fitting Toolbox™ . You can perform exploratory data analysis, preprocess and post-process data, compare candidate models, and remove outliers using this tool. You can also conduct regression analysis using the library of linear and nonlinear models provided or specify your own custom equations. The library provides optimized solver parameters and starting conditions to improve the quality of your fits. It also supports nonparametric modeling techniques, such as splines, interpolation, and smoothing. You can do many other things such as

variety of postprocessing methods for plotting, interpolation, and extrapolation; estimating confidence intervals; and calculating integrals and derivatives. It is a great tool and you will love it.

If your signal frequency is large compared to your sample frequency, wouldn't it be easiest to just check the Fourier Transform?

No one has described the quick techniques to calculate time constants.

For a quick assessment locate the peak, time and amplitude, estimate the point at which the waveform reaches 0.368*the peak amplitude, and the difference in time defines the time constant. Good results depend on a high oversampling ratio, at least 25:1 relative to the time constant, depending on the SNR.

The next quantitative improvement would be a non-linear least squares. Try a generalized least squares fit to the sample points defined by:

{tn, ln(((yn - min(YN)+a)/apk))} where tn are the time points, yn are the associated measurements, min(YN) is the minimum of all the measurements and a is the magic number which makes the vector quantity evaluated in the natural log entirely positive.

Fourier approaches will run into problems with frequency resolution - with the advantage that the result may indicate that you need a higher order model.

If you know that the model for each event is an exponential increase followed by an exponential decrease and you know when the events occur, then nonlinear least-squares fitting of each event is definitely the way to go. There's several ways to do this in MATLAB, but the easiest is probably to just define a cost function (measurement - (prediction from model))^2 and find the time constant that minimizes it for each event using fminsearch. You could get the initial guess by doing linear least-squares on ln(initial value - final value); this value might be good enough for your purposes that you don't need to do the full nonlinear least-squares.

I cannot imagine simple Matlab-type smoothing or filtering will do the job for you. You seem to have very detailed knowledge of the typical shape of your function which rises about the noise. So a generative statistical model would be the way to go. The question now remains what. Assuming your shape doesn't change too much you could apply a matched filter. The matched filter would provide you with the positions of maximal match with the template (which you could have generated by ensemble averaging of several traces) and you could use this knowledge to clean the signal.

I would use Fourier Transforms or bilateral Fourier transform (Laplace) to pick out the decay rate. So, take a length of time for the signal history and Fourier transform, you should be able to pick out or measure the relaxation time in the Fourier codomain.

why not to use match filtering.

If you have an expected shape of the signal try to use filter function and locate the pick of the output.

In case of problem send me mat file

First, assuming that the signal is sampled at 1sps (sample per second), according to the Nyquist theorem the maximum frequency content that can be determined is 0.5 Hz. Second, the sampling rate appears to be quite slow, so that the sampled process/event must also be quite slow as well. I would suggest a couple of strategies.

1) If you require a rough estimation, try increasing the sampling rate (10-12 times higher). Depending on the Signal-to-noise ratio, a simple FIR filter can be used to filter out the noise to reveal the process response; simple averaging using 3, 5, or 7 terms can help to reduce the noise content. If the response resembles that of an exponential decay, simply identify the beginning of the transient response and the beginning of the steady state response (within 1% of the final value). Divide the time difference between the beginning of the steady state and the beginning of the transient response times by 5 to give an estimate of the time constant; some authors use the value just below 2% and dividing by 4.

2) Whether is possible to increase the sampling rate or not, you may try wavelet denoising first. To get you started quickly, load the measured data into a 1-D variable in MATLAB workspace. Start the wavelet GUI by typing wavemenu. Select the “SWT denoising 1-d” button from the GUI. Load your data into the Stationary Wavelet Denoising interface and try out different combination of wavelet and parameters. If the procedure works for you, you can even create your own code for automating the analysis process. Again if you identify a wavelet denoising process that works for your data, export the results (denoised data) back to MATLAB workspace. Start the system identification GUI by typing ident. Load your data and try out different time domain modeling options. You need one more vector corresponding to the input excitation signal; since you are only interested in finding the time constant you can generate a vector of ones (unit step input), the same size as the denoised signal.

There are many other methods that you can try out, for either denoising or identifying the dynamics of time domain data (for instance correlation, cross correlation, DFT, etc). However, it is necessary to have much more information about the process and the nature of the noisy signal. I hope this helps.

If you noise it high frequency than 1Hz then just a simple filter using a series resistor with a capacitor will work as a low pass filter. 16K ohm and 10uF will work. The cutoff frequency is 0.99Hz

Your case looks similar to ECG signal or EEG signal.processing....

Global level of curve fitting looses minute information which are critical for further analysis and vice verse for low level curve fitting using Cubic splines....

For filtering

http://174.142.16.66/csc/manuscript/Journals/SPIJ/volume3/Issue5/SPIJ-48.pdff

for seperation of peaks for time constant identification

1) using a template matching as proposed by Michael Bono & Iead Rezek

or

2)http://www.physik.uni-freiburg.de/~severin/ECG_QRS_Detection.pdf

Thanks everyone for your help. I am trying some of the methods and approaches suggested which I can implement and if I run in to problems I will post further.