Frequency Analysis for Derivative Estimates

Numerical Methods for Derivative Estimation

We have seen the hazards of using the Central Differences method for estimating derivative values in applications other than interpolation of highly accurate, precise data in tables. In this section we want to consider how to distinguish undesirable derivative estimator behaviors from the desirable ones. A frequency analysis is one way to obtain this kind of information.

In a frequency analysis, the input data are restricted to be pure sine or cosine waves, and the output responses of the estimator are examined. The estimator should produce fully predictable results for these special cases, since derivatives of sines and cosines are themselves sines and cosines. Adjusting the input wave frequency adjusts the derivative values, which can be made as large or as small as desired. There are well-developed frequency analysis methods that can assist the investigation.

 

Frequencies in real data

As a first exploratory step, let's take the arbitrary mathematical function used in the previous section, and calculate an FFT of that data sequence, just to get an idea of what frequencies are present — and that the derivative estimator formula will see. The following plots show the "mathematically exact" derivative curve and then its corresponding frequency spectrum.

ideal response curve

Analytical derivative curve

magnitude spectrum curve

Frequency spectrum for the analytical curve

In the spectrum plot, the higher the spectrum curve, the more of that corresponding frequency is present in the function sequence. The further to the right, the higher the frequency. The point at position 1 describes the degree to which the curve "looks like a constant" (frequency zero) function. The point at position 2 describes the degree to which the curve "looks like a wave with one cycle spanned by 100 samples." The point at position 2 describes the degree to which the curve "looks like a wave with two cycles spanned by 100 samples." Continuing on with this, the value at index k describes the degree to which the curve "looks like a sine wave with k-1 cycles spanned by 100 samples."

You might have noticed the peculiar symmetry of the spectrum plot. This is an artifact of having function values at discrete intervals. This phenomenon is known as aliasing and it is always a potential problem. Basically, it means that the second half of the spectrum plot (terms 51 to 100) have no independent meaning, and can be disregarded. The theoretical limit, halfway between frequency 0 and the sampling rate at which samples are taken, is known as the Nyquist Limit.

So, let us now repeat the FFT analysis, this time using the same function data but corrupted with about 1% random noise.

noisy spectrum curve

At first glance, it might appear that there is no difference. But look more closely near the bottom axis. You will see very small random variations. In low frequencies at the left end of the spectrum, relevant behaviors strongly dominate, and noise effects are impossible to distinguish. In the higher frequencies, noise effects are easily distinguishable, because nothing else of relevance is present there.

The main observation to take away from this is the fact that there is nothing of any significance beyond the 10th term of the frequency spectrum, which is at 20% of the Nyquist Limit. In general, it is reasonable to assume that any frequencies beyond this range are harmful to the derivative analysis. That is, variations in the input sequence that do not span at least 10 samples can be considered some kind of noise or data corruption.

 

Frequency characteristic of a derivative operation

What is the response of a theoretically perfect derivative operation when it is applied to a sine or cosine wave? Clearly, sine wave and cosine waves are really the same except for a different phase angle (or somewhat equivalently, a time shift). We can analyze how a differentiator formula responds to both, at the same time, by representing waves using the "complex exponential" function form with both sine and cosine terms:

     e j w t    ≡   cos(w t) + j sin(w t)

where j is the customary notation for the "imaginary component" used in signal analysis, w is the frequency expressed in radians per interval, and the function input t is the position within the waveform. We know that derivatives of a waveform with frequency zero are zero, so we do not need to analyze that frequency. We know that derivatives of the waveform repeat, so we do not need to analyze more than one cycle of the lowest frequency of interest. Responses that include the imaginary multiplier j will be related to the sin function part, and real-valued responses will be related to the cosine function part.

We can express the samples of the sine and cosine waves that we analyze in terms of a large number of discrete "phase angle" steps. The lowest frequency we have chosen for the analysis will have 100 positions analyzed in the course of one waveform — a larger number M could be chosen, but that wouldn't make much difference and adds computational load. Since the function input values can now be indexed by an integer value, we can change notations and use index m to designate the positions along the waveform.

     e j ((f 2*pi m) / M)  

What response should we expect according to a perfect analytical calculation of the derivative for pure waveforms at each frequency? Basic derivative calculus tells us that expected responses are:

cos'((f 2*π m)/M) = - ((2*π m)/M) sin((f 2*π m)/M)

sin'((f 2*π m)/M) =   ((2*π m)/M) cos((f 2*π m)/M)

Or using the complex exponential notation for the derivative, an equivalent and slightly more general expression is

     e' j ((f 2*pi m) / M) = j ((2*pi m) / M) e j ((f 2*pi m) / M)

Taking the ratio of the output to input at each frequency, we can see that the apparent gain factor of the derivative operation at each frequency is

   gain(n) = j 2*pi m / M 

It is clear from this relationship that:

  1. The imaginary value notation j indicates that the phase angle is shifted by π/2 radians.
  2. The response magnitude is proportional to the frequency: 2*pi m/M as indicated by index m.
  3. The gain response for "negative frequencies" with index -m can be seen to be the negative of the response with m positive. It is sufficient to analyze the positive frequencies only, and take into account the negative symmetry later.

 

Comparison: central differences estimator vs. ideal

So, we know how to calculate the frequency response of a proposed filter, by plugging the "complex frequency expression" into the "filter equation." We also know that the "ideal frequency response" for a mathematically perfect differentiator is a "straight line" with a response proportional to frequency. But what response do we want from a practical estimator? We would want to track the "straight line" response accurately in the band from zero up to 20% of the Nyquist, where relevant effects are found. Beyond that, the response should be as small as possible so that irrelevant noise effects are excluded.

Here is a plot of the frequency response characteristic of the Central Differences estimator, shown in blue, superimposed on the desired response characteristic just developed for an ideal estimator formula, shown in green.

estimator comparison

 

Notice how the Central Differences estimator applies an increasingly high gain to the noise effects at upper frequencies, where the desired response is approximately zero. It should now be clear where the extreme sensitivity of the Central Difference estimator comes from.

 

What can we do about this?

One simplistic solution is based on the old vaudeville joke:

Patient: Doctor, please help me! My arm hurts when I raise it up like this!

Doctor: Well then, don't raise your arm up like that!

It turns out that this kind of avoidance strategy can be relatively effective. More about this idea in the next section.

 


Footnotes:

 

Contents of the "Numerical Estimation of Derivatives" section of the ridgerat-tech.us website, including this page, are licensed under a Creative Commons Attribution-ShareAlike 4.0 International License unless otherwise specifically noted. For complete information about the terms of this license, see http://creativecommons.org/licenses/by-sa/4.0/. The license allows usage and derivative works for individual or commercial purposes under certain restrictions. For permissions beyond the scope of this license, see the contact information page for this site.