Prosig Signal Processing Blog

By Dr Colin Mercer, Technical Director, Prosig

Differentiation and integration are often used to compute between accelerations, velocities and displacements. These relatively straightforward mathematical operations may cause surprises when carried out numerically. The available frequency range over which reliable results are obtained is surprisingly small. This brief note illustrates the problems which may occur and indicates techniques to resolve the problems.
The basic method of differentiating a signal numerically is to estimate the local slope of the curve. The three simplest schemes use differences between values. If Δt is the time step between adjacent samples the three most well known schemes are

forward differences dx_j+0.5 = (x_j+1 - x_j) / Δt,

backward differences dx_j-0.5 = (x_j - x_j-1) / Δt,

central differences dx_j = (x_j+1 - x_j-1) / (2Δt).

The forward difference and backward difference forms have a problem in that they actually estimate the slope at the midpoint between samples. Hence the use of the subscripts j+0.5 and j-0.5 respectively. Thus in practice the simplest practical scheme is to use the central difference form.

What is often not appreciated is that any differentiation process is the equivalent of applying a filter to the data. The ideal differentiator^* has a frequency characteristic of D(iω) = (0 +iω). In modulus and phase form this is (ω, 90°). Thus on a normalised frequency axis the modulus should be a straight line from 0 to p and the phase should be at 90º.

[footnote - ^* this is known as a Hilbert Transformer but for the present purposes this is not significant].

To find the actual frequency response function use a broadband Gaussian Random Signal, differentiate it, form the Cross Spectral Density of the differentiated signal and the original signal and divide this by the Auto Spectral Density of the original signal.

How do we proceed if we need accurate integration or differentiation?There are 2 approaches. One is to use integration or differentiation in the frequency domain where we manipulate by using (1/iω) or (iω) respectively. This is referred to as “Omega Arithmetic”. Basically we FFT the signal, use Omega Arithmetic to integrate or differentiate and then transform the result back to the time domain. It is of course essential with the integration in particular that we first reduce the signal to a zero mean. Carrying out this operation gives the correct results and, as would be expected, the measured transfer functions for this method matches the ideal ones. We repeat however that eliminating the low frequencies is essential.The second approach is to use the Time Domain Interpolation function. This upsamples the signal by a factor up to 100 times. Typically upsampling by a factor of say 16 is generally sufficient. Upsampling is preferred to initially oversampling the original data. With upsampling good results are obtained up to 0.4 or 0.45 times the original sample rate (0.8 to 0.9 times Nyquist).For integration the low frequency problem is often referred to as “1/f” noise. There is no real solution except to high pass filter and remove the very low frequencies. If it is necessary to determine the integration at very low frequencies then the approach is to decimate the signal, again with a technique which maintains frequency content true (eg the Decimate module in the Time Series section).