How to properly take derivative of discrete data numerically in matlab?

Question

Taking derivative of discrete data requires some fitting and then using the b-form of the polynomial.

How to make a 'good' fit and properly take derivative? I am confused what is the right behavior of the derivative. Things change quite rapidly with small adjustments of the fitting csaps coefficient 0.09 in the example below.

The test.txt data is attached. Test data test.txt click

Matlab code.

T = load('test.txt'); x = T(:,1); y=T(:,2)
pp22 = csaps(x,y,0.09);%attempt to create 'good' fit in b-form. ... Notworking as intended
pp22val = ppval(pp22,x);
figure(1012)
plot(x,pp22val,x,y) %quality of the fit
p_der202 = fnder(pp22,2); % second derivative
y_prime22 = ppval(p_der202,x); % values

figure(101)
plot(x,y_prime22) %second derivative

% test discrete centered derivative
for i = 4:length(x)-4
prime23(i) = (-y(i+2)+16*y(i+1)-30*y(i)+16*y(i-1)-y(i-2))/(12*((x(i+1)-x(i))^2));
end
figure(102)
plot(x(2:length(x)-3),prime23)

Data and fit

Matlab fnder

Discrete centered

Answer 1

Splines tend to consider that control points are not noisy, they can induce some overshoot. Since your time sampling is almost regular, I would suggest to keep the denominator constant (median oof the differences). Plus, your data looks very smooth. Resultingly, you can use a lot of standard linear filtering tools. To name a few:

• finite difference approximations of the second derivative, a procedure is here, your operator is numbered (13),
• filter the signal with a Gaussian kernel, or a discrete approximation, and differentiate,
• use Savitzky-Golay filters, implemented in Matlab under sgolayfilt and sgolay, which give you access to the derivatives quite directly.

The last two options seem appropriate to me. What is important the the choice of the scale under which the derivatives are meaningful.

I did a try, adapting Matlab code. On its right end, the derivative seems blocky (piecewise constant), suggesting a close to piecewise linear signal, hence the peaks in your second derivative. You will have to play with parameters.

N = 3;                 % Order of polynomial fit
F = 41;                % Window length
[b,g] = sgolay(N,F);   % Calculate S-G coefficients

T = load('SE_Math_1594190_test.txt'); x = T(:,1); y=T(:,2)

dx = median(diff(x));

SG0 = zeros(size(x));
SG1 = zeros(size(x));
SG2 = zeros(size(x));
HalfWin  = ((F+1)/2) -1;
for n = (F+1)/2:length(x)-(F+1)/2,
  % Zero-th derivative (smoothing only)
  SG0(n) =   dot(g(:,1), y(n - HalfWin: n + HalfWin));

  % 1st differential
  SG1(n) =   dot(g(:,2), y(n - HalfWin: n + HalfWin));

  % 2nd differential
  SG2(n) = 2*dot(g(:,3)', y(n - HalfWin: n + HalfWin))';
end

SG1 = SG1/dx;         % Turn differential into derivative
SG2 = SG2/(dx*dx);    % and into 2nd derivative

% Scale the "diff" results
DiffD1 = gradient(y)/ dx; 
DiffD2 = gradient(gradient(y)) / (dx*dx);

subplot(3,1,1);
plot([y(1:length(SG0)), SG0])
legend('Data','S-G Smoothed data')
axis  tight
subplot(3, 1, 2);
plot([DiffD1,SG1])
legend('Grad-generated 1st-derivative', ...
'S-G Smoothed 1st-derivative')
axis  tight
subplot(3, 1, 3);
plot([DiffD2,SG2])
legend('Grad-generated 2nd-derivative', ...
'S-G Smoothed 2nd-derivative')
axis  tight