I once had the problem that i had to differentiate noisy data - which is a problem. The differentiation process basically behaves like a high-pass filter so that the high frequency noise becomes dominant. The problem gets worse when you not only have to make the first but also higher derivatives.
Lets look at an example: we take a simple sin function and add a small amount of noise. Figure 1 shows the resulting curve. The noise is hardly visible, so Fig.2 shows a magnified version.
 |
| Fig.1 - a sine curve with some noise |
 |
| Fig.2. Here the noise curve (circles) is compared with the analytical sine curve. |
 |
| Fig.3 second derivative when the gradient function is used. |
How can one improve the derivative? One way is to make a fourier filtering afterwards. Another way is to apply a "smoothed" or "filtered" gradient. The basic idea here is to include more than 2 points in the differentiation process.
I'm not very good a theory. There is a very sophisticated explanation on the page of
Pavel Holoborodko
("Numerical Methods / Numerical Differentiation"). Here i found two Algorithms I implemented for octave.
The first one is the "low noise Lanczos Gradient":
function [x2,y2] = lanczos_gradient(x1,y1,DOrder);
%
% -----------------------------------------------
% usage :
% [x2,y2] = lanczos_gradient(x1,y1,o);
%
% parameters :
%
% x1 : x values (x2 will be cropped by 2 times o)
% y1 : function values to be list to differenciated
% o : order of Lanczos-Differentiation (odd number)
%
% IMPORTANT : x1-values have to be aequidistant !
%
% -----------------------------------------------
%
NumPoints = length(y1);
h = x1(2)-x1(1);
y2=[];
m = (DOrder-1)/2;
for i=m+1:NumPoints-m
NewValue = 0;
for k = 1:m
%printf("\n NumPoints %d und %d ",i+k,i-k)
NewValue = NewValue + 3/h * k * (y1(i+k) - y1(i-k)) / (m * (m+1) * (2*m+1));
endfor
y2 = [y2,NewValue];
endfor
x2=x1(m+1:NumPoints-m);
endfunction
The second one is the "Smooth noise-robust differentiator":
function [x2,y2] = filtered_gradient(x1,y1,DOrder);
%
% -----------------------------------------------
% usage :
% [x2,y2] = filtered_gradient(x1,y1,O);
%
% parameters :
%
% x1 : x values (x2 will be cropped by 2 times O)
% y1 : function values to be list to differenciated
% O : order of the Filter (an odd integer number)
%
% IMPORTANT : x1-values have to be aequidistant !
%
% -----------------------------------------------
%
NumPoints = length(y1);
h = x1(2)-x1(1);
y2=[];
m = (DOrder-1)/2;
for i=m+1:NumPoints-m
NewValue = 0;
for k = 1:m
%printf("\n NumPoints %d und %d ",i+k,i-k)
NewValue = NewValue + 3/h * k * (y1(i+k) - y1(i-k)) / (m * (m+1) * (2*m+1));
endfor
y2 = [y2,NewValue];
endfor
x2=x1(m+1:NumPoints-m);
endfunction
For details how these two work, please have a look on Pavels page. If you're not sooo much into theory (as i do) let's just apply these two on data. For the example above, the result is shown in Fig. 4
 |
| Fig.4 the two Algorithms give derivatives very close the the noise-free case |
As can be seen, the curves look much better.
Please keep in mind: we can apply these algorithms to our data. And yes, the results look very good (especially when we know the analytical result). But be carefull: by "smoothing" noise out, we alway reduce the amount of information in our data. So it's a tradeoff and one has to be very careful.
Have fun playing around with this.
Here's the code to produce the pictures:
clear all;
more off;
LO = 11; % the 'order' of smoothing
x = linspace(0,8,1000);
noise = 0.0001*(rand(1,1000)-.5);
y = sin(x)+noise;
% create the filtered gradients
[xll,yll] = filtered_gradient(x,y,LO);
[xl,yl]=lanczos_gradient(x,y,LO);
[xll2,yll2] = filtered_gradient(xll,yll,LO);
[xl2,yl2] = lanczos_gradient(xl,yl,LO);
%create the simple gradient
ysimp = gradient(y)/(x(2)-x(1));
ysimp2 = gradient(ysimp)/(x(2)-x(1));
figure(1)
plot(x,y, "linewidth",2)
%title "sin(x) + small amount of noise"
xlabel "x"
ylabel "y"
ax(1) = gca();
figure(2)
plot(x,sin(x),"linewidth",2,
x,y,"o","markersize",6,"linewidth",2)
legend("sin(x)","noisy sin(x)")
xlim([1.5,1.65])
ylim([0.996,1.001])
xlabel "x"
ylabel "y"
ax(2) = gca();
figure(3)
plot(x,ysimp2)
xlabel "x"
ylabel "y''"
ax(3) = gca();
figure(4)
plot( x, -y, "linewidth",2,
xl2, yl2,"linewidth",2,
xll2, yll2,"linewidth",2)
xlabel "x"
ylabel "y''"
legend("-sin(x)","lanczos gradient","filtered gradient")
ax(4)=gca();
set(ax,"linewidth",2,"tickdir","out","ticklength",[0.002,0.002])
set(ax([1,3,4]),"ylim",[-1.2,1.2])
%create the plots
figure(1)
print("sinnoise.png","-dpng","-F:Helvetica:11","-S900,550")
figure(2)
print("sinnoisedetail.png","-dpng","-F:Helvetica:11","-S900,550")
figure(3)
print("sinnoisesimple_gradient.png","-dpng","-F:Helvetica:11","-S900,550")
figure(4)
print("sinnoisegradients.png","-dpng","-F:Helvetica:14","-S900,550")
