Fit Function for Data with Second Derivative Sometimes Flat

Hi,

Disclaimer: I'm a physicist, not a real coder!

I have beautiful data showing motion of each part of a stretchy particle:

Now I need to differentiate it. I currently try two different methods, diff for a ballpark figure, and polyfit-->polyder for an accurate measure.

An example for one curve:

%% velocity and acceleration by polyfit
 
[fitCoM,erCoM,muCoM] = polyfit(t,yCoM,7); %fit polynomial to CoM
vfitCoM = polyder(fitCoM)./1000; %derivative to poly fit
afitCoM = polyder(vfitCoM); %derivative to poly fit

%% velocity and acceleration by finite differences

vCoM = diff(yCoM).*framerate./1000; %measured vel of CoM (m/s)
aobsCoM = diff(smooth(vCoM,9,'moving')).*framerate; %measured acc of CoM (m/s^2)

I only show the polyfit part for acceleration (shown as 'observed', the 'calculated' part is a different model).

Minimal working code and data is attached.

This looks reasonable, right?

For that data ('Difftest1.csv'), yes. However, when I apply it to simpler data ('Difftest2.csv') under gravity only, the second derivative should be flat, but I get crazy results:

I've tried a few different modifications to my polyfits (degree, centring and scaling), but none make a major difference here.

Is there a better fit function I could apply consistently for differentiation, when some of my data has a constant second derivative and some doesn't?

Thanks :)

Answer 1

[Star Strider]

It’s difficult for me to follow your code.

Rather than using polyfit to smooth your data, I would use the smoothdata function (introduced in R2017a), specifically using the 'movmean' or 'sgolay' methods. (I believe the 'sgolay' method requires the Signal Processing Toolbox, although that is not in the smoothdata documentation.) Taking the derivative of a noisy signal amplifies the noise, so smooth it first. Then, use the gradient function to take the numerical derivative.

[Frederick Wells]

I have tried a few combinations of gradient(smoothdata(...)), with different methods etc.

The problem with smoothdata is that I always see artefacts approximately one window-length from each end. These translate into bigger and bigger issues on differentiation.

Eg. for

vCoM = gradient(smoothdata(yCoM,'lowess')).*framerate./1000; %measured vel of CoM (m/s)
vLQ = gradient(smoothdata(yLQ,'lowess')).*framerate./1000; %measured vel of lower quartile (m/s)
vUQ = gradient(smoothdata(yUQ,'lowess')).*framerate./1000; %measured vel of upper quartile (m/s)

aCoM = gradient(smoothdata(vCoM,'lowess')).*framerate./1000; %measured vel of CoM (m/s)
aLQ = gradient(smoothdata(vLQ,'lowess')).*framerate./1000; %measured vel of lower quartile (m/s)
aUQ = gradient(smoothdata(vUQ,'lowess')).*framerate./1000; %measured vel of upper quartile (m/s)

The acceleration plots should all be flat at -9.8.

They are relatively accurate numerically, but shape always seems to be distorted by smoothdata...

[Star Strider]

The smoothdata function is by definition a lowpass filter, using different methods to achieve the result. Filtering always removes details, so that may be what you are seeing. You can experiment with various frequency-selective filters (the Signal Processing Toolbox makes this relatively straightforward), or use wavelet denoising techniques (Wavelet Toolbox) although I do not have extensive experience with wavelets. Those may be preferable to smoothdata for what you are doing.