# Frequency response of numerical derivative

Analytical derivative of a function is equivalent to convolution of that function with $$s$$ in Laplace domain. Numerical derivatives are limited in bandwidth due to finite sampling rate, so they are not synonymous with convolving the signal with with $$s$$ term. At higher frequencies one would expect attenuation of the numerically differentiated signal from one that was computed analytically. Recently, I found that there are some differences at the low frequency limit as well which I cannot explain.

Attached is a plot of a signal sampled from a normal distribution (blue) and it's first derivative in time (red). As expected, at high frequencies the derivative signal begins to attenuate. But why does it not cross $$\omega$$ = 1 rad/s or 0.16 Hz as would be the case if the solution was obtained analytically?

Here's the code I am running in MATLAB

sr = 100000;
y = randn(1,sr);
dydt = y;
for i = 2:length(y)-1
dydt(i) = (y(i+1)-y(i-1))*sr*2;
end
hold on, plot(abs(fft(y)));
plot(abs(fft(dydt)));
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')

• How did you filter the signal? Finite differencing or something else? – fibonatic Sep 30 '19 at 3:11
• The low frequency 'mis'-behavior happens regardless of type of differentiation used. In this case I used a symmetric derivative though: en.wikipedia.org/wiki/Symmetric_derivative . – Roman Vas Sep 30 '19 at 3:39
• Might you be running into rounding errors? For example what happens when lower the sample frequency significantly? – fibonatic Sep 30 '19 at 3:46
• It does appear to be sampling rate dependent. The 'corner' frequency of derivative (red) scales with sampling rate. Could you explain how the rounding error contributes to the observed behavior? Thanks – Roman Vas Sep 30 '19 at 3:55

Even with @MattL.'s fix you are discarding typically non-zero parts of the discrete-time derivative by not including its first and last sample, which destroys its autocorrelation properties near the end points, typically resulting in the low-frequency plateau in the frequency spectrum as you have observed. We can add a bit of a zero-valued safety buffer at the start and at the end of the signal to ensure that what we are discarding will be zero-valued:

sr = 100000;
y = randn(1,sr);
y(1) = 0;
y(2) = 0;
y(end-1) = 0;
y(end) = 0;
dydt = zeros(1,sr);
for i = 2:length(y)-1
dydt(i) = (y(i+1)-y(i-1))*sr*2;
end
hold on, plot(abs(fft(y)));
plot(abs(fft(dydt)));
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')


The result is as desired:

Figure 1. Result using safety buffers at the start and at the end of the signal.

Another way is to treat the signal as periodic and to wrap around the subscripts:

sr = 100000;
y = randn(1,sr);
dydt = zeros(1,sr);
dydt(1) = (y(1+1)-y(end))*sr*2;
for i = 2:length(y)-1
dydt(i) = (y(i+1)-y(i-1))*sr*2;
end
dydt(end) = (y(1)-y(end-1))*sr*2;
hold on, plot(abs(fft(y)));
plot(abs(fft(dydt)));
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')


This will also give the desired result:

Figure 2. Results when treating the signal as periodic.

• More generally, given any signal it looks like we have to make sure it's starting and ending values are the same before evaluating fft. It's easiest to just set them to 0 granted the signal in question is long enough. Thanks for your input! – Roman Vas Oct 1 '19 at 20:11

When I try this, the results looks as expected. So you really have to explain in more detail what exactly it is that you're doing, because it doesn't seem to be a property of the discrete-time derivative.

EDIT: Now that I see your code, I'm convinced that the problem will disappear if you use

dydt = zeros(1,sr);

to initialize the derivative vector.

• That's very strange! Here's the code I am running in MATLAB: sr = 100000; y = randn(1,sr); dydt = y; for i = 2:length(y)-1 dydt(i) = (y(i+1)-y(i-1))*sr*2; end hold on plot(abs(fft(y))); plot(abs(fft(dydt))); set(gca, 'YScale', 'log') set(gca, 'XScale', 'log') – Roman Vas Sep 30 '19 at 13:33
• @RomanVas It is probably because you do not overwrite the first element of dydt, which essentially adds a impulse to the derivative approximation which has a flat spectrum, which dominates at the low frequencies. – fibonatic Sep 30 '19 at 14:36
• @fibonatic The issue persists even if I don't include the first and last points of dydt in my evaluation of the Fourier transform. – Roman Vas Sep 30 '19 at 14:58
• @Matt L. Nope! That doesn't fix a thing. What is going on?!? I am so confused by this issue. – Roman Vas Sep 30 '19 at 21:29
• Yep this is only a partial fix, see my answer for the rest. – Olli Niemitalo Oct 1 '19 at 5:55

Not an answer but too long for a comment:

It looks like you may have two problems: First you seem to plot the negative frequencies as the top half of the spectrum. You should only plot the output of the FFT from $$0$$ to $$N/2$$.

Second: did you a try a symmetric differentiator? The simplest approach to in the discrete domain is simply. $$y_0[n] = x[n]-x[n-1]$$

That induces more or less a half sample delay. More precise would be something like $$y_1[n] = x[n+1/2]-x[n-1/2]$$

but that's a lot harder to implement. Is't basically the simple differentiator cascaded with a half sample delay of -0.5. When designing fractional delays, you can run into the problems that you show.