MATLAB: What is the best way to smooth and compute the derivatives of noisy data

curve fittingdata smoothingderivative calculatormaximum slope

Hello all,

I have a set of experimental data (temperature vs time), and have no idea of the function which satisfies the data behaviour. As a part of analysis, I need to estimate the first derivative using central differences and obtain the maximum gradient.

Because my data is too noisy, I need to filter it before taking the derivative. To keep the precision of data and minimize any distortion, I tried to remove the outliers from my data using a Savitzky–Golay filter. It shows great results, but my data is not quite smoothed as it can be seen in a picture of Savitzky–Golay filter. To compute the exact inflection point further curve fitting is required, as it can be seen in the first derivative picture.

Any suggestions or better approaches , would be highly appreciated as I’m new to Matlab and signal processing. I heard about the cubic spline for better smoothness, but I have no idea about it and not sure if it can solve my problem.

Many thanks in advance

Best Answer

The "best" way? That will be completely dependent on your data, on what you know about your data, and on what you are willing to assume about the process that generated the data.

Looking at your data, since I don't have your data, nor do I know anything about the data or the process, I can only guess. Thus, can I assume that what looks like noise really is noise? I had one group of clients who were not willing to make that assumption. Any smoothing at all was a bad thing, because there was potentially information in the high frequencty stuff that they needed to know about. Luckily for them, the noise was not overwhelming.

In general though, numerical differentiation is a noise amplifying process. So you need to kill off the noise to whatever extent possible. You can do so by use of a spline fit (IF you use the right spline model. An interpolating spline would be generally a bad idea for you), or you could use a regression model of some sort, IF you have a viable model for this process.

Again, I don't have your data. So, at best, I can only make some data up for an example. So, let me see, this should look at least a little like what you have shown:

x = linspace(0.5,2.5,100);
y = 19 + exp((x - 0.5)/2)/10 + randn(size(x))/100;
plot(x,y,'.')

The picture you show has a fair amount of data, yet it is not that complex of a relationship. As such, a low order polynomial model might be sufficient. Then you could just differentiate the polynomial. For example...

P3 = polyfit(x,y,3)
plot(x,y,'b.',x,polyval(P3,x),'r-')

P3d = polyder(P3)
plot(x,polyval(P3d,x),'-')

A problem with such an approach is it does not build information that you may have about the underlying functional relationship. That is, do you know that temperature will ALWAYS increase with time? If so, then the first derivative must always be positive, yet a polynomial fit gives you no such assurance. In fact, you might notice that the true underlying function behind my data has an everywhere increasing derivative, yet a cubic polynomial model had a negative second derivative. Worse, what if you use too high an order for a polynomial? That is easy to do.

P7 = polyfit(x,y,7);
plot(x,y,'b.',x,polyval(P7,x),'r-')

The same fact applies to a Savitsky-Golay filter. A Savitsky-Golay derivative estimation on noisy data will need to have a long window width, and low order implicit polynomial model. However, remember that Savitsky-Golay does not produce a differentiable or continuous function/derivative estimate. So, using my own movingslope code, that builds polynomial model on a moving window, then returns the derivative, we see:

plot(x,movingslope(y,20,1),'.')

Anyway, if we look at the bottom end of the curve, the 7th degree polynomial started to get a little squirrely. So, I might try to use a regression spline, constrained to have the expected shape. Lets see,

slm = slmengine(x,y,'plot','on','increasing','on','concaveup','on','jerk','positive');

So here, a cubic spline, constrained to have first, second, and third derivatives entirely positive over the domain.

Better yet, of course, is to use a nonlinear model, IF you know a good choice of model in advance.

But really, what matters is as I said in the beginning. What do you know? Because once you can provide information to the modeling process, you will GREATLY improve the derivative estimation process. Tools like Savitsky-Golay assume almost nothing about the system. A polynomial model assumes relatively little, but somewhat more.

Related Solutions

MATLAB: How to interpolate data using quartic or higher order polynomials? Similar to the cubic spline interpolator csapi

You don't really want to do so. Yes, you think you might. But you don't.

High order polynomials are easy enough to build, polyfit can generate them. Just use one order less than the number of data points. So two points exactly determine a line, 5 points a quartic polynomial, etc.

That is the good news. The bad news is high order polynomials will almost always cause you to be forced to post a new anguished question, like "Why does my high order polynomial interpolant do insane things?"

In fact, depending on your data, even something relatively low order like a cubic polynomial can cause problems. Those problems typically arise from two things. The first is poorly scaled data. So if your independent variable (I'll call it x) varies over the interval [0,1000], then cubing those numbers will result in values on the order of 1e15 for a 5th order polynomial. Now, when you add and subtract numbers that vary by 15 powers of 10, expect to have numerical problems in double precision arithmetic.

High order polynomials say 15 or so, almost always tend to have numerical problems, because now you are raising your numbers to seriously high powers. The linear algebra needed to compute those coefficients often starts to have problems now, because again, it works in double precision arithmetic. So unless you are very careful about issues like scaling and centering of your data, expect problems.

Next, polynomials always SEEM like a good idea. (Hey, I can think of a lot of things I've done that SEEMED like a good idea at the time.) After all, a Taylor series is just a polynomial, and they can represent almost anything. But the fact is, polynomials are squirrelly things. They squirm around, introducing lots of added extrema to data that never did anything bad.

The classic example is seen here:

x = -10:10;
y = 1./(1+x.^2);
p = polyfit(x,y,20);
Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as
described in HELP POLYFIT. 
> In polyfit at 75 
ezplot(@(x) polyval(p,x),[-10,10])
hold on
plot(x,y,'ro')

It turns out that the squiggles you see are not due to the warning message, although scaling the data would have allowed it to work without generating a warning message.

And NEVER use a polynomial to extrapolate. Ok, linear extrapolation tends to be moderately safe. Stop there though.

Set('soapboxmode','off')

MATLAB: Find tangent to set of 2D points

Your curve is somewhat noisy. I might think to point out that my interparc tool from the file exchange can be made to give estimates of the derivatives of the curve at any point along the curve. So it would return two numbers at a point that are essentially dx/dt and dy/dt. The ratio is the slope of the tangent line at that point.

The above works even for non-functional closed curves as you have, since interparc works in a parametric form. It computes a simple cumulative linear chordal arclength along the curve, then uses spline models x(t), y(t).

The problem is, your data is somewhat noisy, based on the curve we are shown. So that tangent line would be pointing all over the place.

Anyway, if you were to use interparc in that way, then it would be just as simple to do the work yourself. You do not provide your data. But, suppose you have data that is in the form of two vectors x and y, sorted in order along such a path in the (x,y) plane?

t = cumsum([0;diff(x(:)).^2 + diff(y(:)).^2)]);

So, t is the cumulative piecewise linear distance between each point, one to the next along your curve.

splx = spline(t,x);
sply = spline(t,y);

Differentiate each spline, then evaluate the ratio of the derivatives at the original points. That is the slope of the tangent line, at that point. If you have the curve fitting toolbox, you can use fnder to differentiate the splines.

tanslope = ppval(fnder(sply),t)./ppval(fnder(splx),t);

If no curvefitting toolbox, the derivative function is trivial to generate for the pp form that spline returns.

Again, the issue is that you may wish to produce a smoother estimate of the tangent slope. If you provide your data, and you need such a smoothed tangent line, just attach the data to a .mat file as a comment.

Best Answer

Related Solutions

MATLAB: How to interpolate data using quartic or higher order polynomials? Similar to the cubic spline interpolator csapi

MATLAB: Find tangent to set of 2D points

Related Question