# 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*.

** Analytical derivative 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.

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:

- The imaginary value notation
`j`

indicates that the phase angle is shifted by π/2 radians. - The response magnitude is proportional to the frequency:
`2*pi m/M`

as indicated by index`m`

. - 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.

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: