Blogs

Evaluate Noise Performance of Discrete-Time Differentiators

Neil RobertsonMarch 28, 20222 comments

When it comes to noise, all differentiators are not created equal.  Figure 1 shows the magnitude response of two differentiators.  They both have a useful bandwidth of a little less than π/8 radians (based on maximum magnitude response error of 2%).  Suppose we apply a signal with Gaussian noise to each of these differentiators.  The sinusoidal signal with noise is shown in the top of Figure 2.  Signal frequency is π/12.5 radians.  The output of the so-called central-difference differentiator is shown in the middle plot, where we see that the signal has been attenuated more than the noise, resulting in obvious degradation.  The output of the “41 tap” differentiator is shown in the bottom plot.  Clearly, the noise performance of this latter differentiator is superior.

Referring back to Figure 1, it’s apparent the central-difference differentiator is passing much more noise energy at frequencies above π/8 radians than the 41 tap one.  The spectrum at the input of the differentiators is shown in Figure 3 (top).  The spectrum at the output of the central-difference differentiator is shown in Figure 3 (middle), where we see that the signal is attenuated about 12 dB, while the noise is barely reduced over much of the spectrum.  You could say the central-difference differentiator is an “open barn door” for noise.  The spectrum at the output of the “41 tap” differentiator is shown in Figure 3 (bottom):  the signal is also attenuated about 12 dB, but the noise power is much lower than that of the central-difference differentiator.

It would be nice to be able to quantify the difference in performance between differentiators.  In this article, I’ll introduce a Gaussian noise performance measure for differentiators that I’ll call Differentiator Noise Power Ratio.

This article is available in PDF format for easy printing

Figure 1.  Magnitude responses of a central-difference differentiator and a 41-tap

differentiator with equal usable bandwidth. 


Figure 2.  Differentiator Noise Performance.

Top:  differentiator input signal with added Gaussian noise.

Middle:  output of central-difference differentiator (note change in vertical scale).

Bottom:  output of “41-tap” differentiator.


Figure 3.  Differentiator spectra.

Top:  input sinusoid with noise.   Middle:  Central-Difference output.   Bottom:  “41 tap” output.


Discrete-Time Differentiator Review

For continuous-time signals, the differentiation property of the Continuous Fourier Transform is:

$$\frac{d}{dt}x(t) \iff j\Omega X(\Omega) \qquad (1) $$

where Ω is continuous radian frequency.  There is no differentiation property of the Discrete Fourier Transform (DFT).  But given discrete signal x(n), it appears from Equation 1 that the DFT of its derivative should be jω X(ω), where ω= 2πf/fs is discrete radian frequency, and fs is the sample rate.

It turns out we can use a Finite Impulse Response (FIR) filter with real coefficients b(n) as an approximate differentiator:

$$ b(n) \circledast x(n) \iff H(\omega)X(\omega) \approx j\omega X(\omega) \qquad (2) $$

where $\circledast$ indicates discrete convolution, and the bi-directional arrow now represents the DFT.  The approximation sign means that the derivative is not exact.  What are the coefficients b(n)?  From Equation 2, we have:

$$ b(n) \iff H(\omega) \approx j\omega \qquad (3) $$

So, b(n) are coefficients whose DFT H(ω) is approximately jω.  Normally, the goal is that H(ω) approximate jω over some bandwidth less than 0 to π (0 to fs/2).  Since jω is pure imaginary, we know from the symmetry properties of the DFT that b(n) must be odd-symmetric.  But odd-symmetric time functions are non-causal, so we typically only require odd symmetry with respect to the center of b(n).

A simple, often-used differentiator is the central-difference differentiator, who’s difference equation is:

$$ y(n)= 0.5x(n)-0.5x(n-2) \qquad (4) $$

The coefficients and block diagram are shown in Figure 4 a and b.  The z-domain transfer function is:

$$ H(z) = 0.5 - 0.5z^{-2} \qquad (5) $$

Letting z = ejω, H(ω) is given by:

$$ H(\omega)=je^{-j\omega} sin(\omega) \qquad (6) $$

where the complex exponential term occurs because the center tap is delayed 1 sample from that of an odd-symmetric time function.  If need be, we can avoid this delay as follows.  The coefficients and block diagram of an odd-symmetric version of the differentiator are shown in Figures 4c and 4d.This filter has output y(n) time-aligned with the delayed input x(n-1).  The transfer function with respect to x(n-1) is:

$$ H(z) = 0.5z - 0.5z^{-1} \qquad (7) $$

resulting in

$$ H(\omega)= jsin(\omega) \qquad (8) $$

Note that the magnitude of H(ω) is the same for either differentiator:

$$ |H(\omega)| = |sin(\omega)| \qquad (9) $$

Recall that our goal is that H(ω) approximate jω, or at least that |H(ω)| approximate ω.  Equation 9 is plotted in Figure 1 as the central-difference response, along with |H| = ω (dashed line).  We see that the response approximates ω for small values of ω:

$$ |H(\omega)| \approx \omega,\; \omega<< \pi \quad (10) $$

Figure 4.  Central-Difference Differentiator

a) Coefficients.   b) Block diagram.

c) Odd-symmetric Coefficients.   d) Odd-symmetric block diagram.


Differentiator Noise Power Ratio

Now we’ll develop our noise performance metric for differentiators.  As an example, we’ll use an FIR differentiator designed using the method described in [1].  I provide a Matlab function diff_synth.m for this purpose in Appendix A of the PDF version.  Suppose we want a differentiator with bandwidth approaching π/2 radians.  The desired magnitude response is shown in Figure 5 (top).  For number of taps = 25, diff_synth.m returns the coefficients shown in the bottom of Figure 5 (We’ll go into more detail about diff_synth.m later in this article).

The magnitude response of this filter is plotted in Figure 6 (top), along with the ideal response.  As you can see, the bandwidth of the filter is less than our original goal.  The response percent error is given by:

$$ \%\, error = 100 \frac{|H(k)|-\omega(k)}{\omega(k)} \qquad (11) $$

The percent error is plotted in Figure 6 (middle).  We can define the useful bandwidth of the filter based on some maximum allowed error.  For this example, if we choose maximum error of 2%, we get a bandwidth of ωmax ~= 0.39π as shown.

Figure 6 (bottom) plots |H(k)|, along with the response of an ideal differentiator with bandwidth ωmax.  The actual response obviously passes more noise power than the ideal response.  I define the quotient of Gaussian noise power passed by these two responses as the Differentiator Noise Power Ratio.

Figure 5.  Differentiator Example.   Top:  Desired response.

Bottom:  Coefficients of 25-tap differentiator.


Figure 6.  Differentiator magnitude response and error.

Top:  Magnitude response and ideal response.

Middle:  Magnitude response error, showing frequency ωmax for 2% error.

Bottom:  Ideal response over 0 to ωmax and actual response.


Now we’ll develop a formula for Differentiator Noise Power Ratio.  Figure 7 shows a discrete-time differentiator with Gaussian noise input having a noise power density of Nb W/bin.



Figure 7.  Discrete-time differentiator with input noise density Nb W/bin.


At the output of the differentiator, we have:

$$ output\,noise\,power\,density= N_b|H(k)|^2 \quad W/bin \qquad (12) $$


The average output noise power is the sum of the noise power density over k:

$$ N_{out}= N_b \sum_{k=0}^{M/2-1}|H(k)|^2 \quad W \qquad (13) $$


where M/2 is the number of bins in the one-sided frequency response H(k).  We define the usable bandwidth of the differentiator as ωmax.  If we now consider the ideal differentiator with bandwidth ωmax (see Figure 6, bottom), the average output noise power is just:

$$ N_{out\,ideal}= N_b \sum_{k=0}^{k_{max}}\omega(k)^2 \quad W \qquad (14) $$


where kmax = Mωmax/(2π).  Taking the quotient of equations 13 and 14, we define the Differentiator Noise Power Ratio R as:

$$ R= \frac{N_{out}}{N_{out\,ideal}}= \frac{\sum_{k=0}^{M/2-1}|H(k)|^2} {\sum_{k=0}^{k_{max}}\omega(k)^2} \qquad (15) $$


For an ideal differentiator, R = 1.  In Appendix B of the PDF version, I derive a simpler expression for R based on ωmax and b(n):

$$ R = \frac{N_{out}}{N_{out\,ideal}} = \frac{3\pi}{\omega_{max}^3} \sum b(n)^2 \qquad (16) $$


We can state that The Signal to Noise Ratio (SNR) at the output of a given differentiator is reduced from that of an ideal differentiator of equal ωmax by a factor of R.  Equations 15 and 16 require that the passband slope of the differentiator magnitude response be 1.0.  It may sometimes be necessary to scale the coefficients b(n) to achieve this.  Keep in mind that the Differentiator Noise Power Ratio calculation is based on an assumption of Gaussian (white) noise.

What about the case of a differentiator cascaded with an FIR lowpass?  We can just convolve the coefficients of the differentiator with those of the lowpass, and use the resulting coefficients in Equation 16, making sure to include any change to ωmax.


Simple Differentiators

Table 1 lists Differentiator Noise Power Ratio R for several simple differentiators, where R was computed from Equation 16, and ωmax is based on 2% response error.  The magnitude responses of the differentiators are plotted in Figure 8.  Note that you don’t get to choose ωmax:  it is fixed for each differentiator. 

Using a different response error criterion would result in different values for R; for example, for a 5% error criterion, the central-difference differentiator would have ωmax = 0.175π and R = 28.3 (14.5 dB).  Also, for a narrowband signal, the differentiator response need only approximate linearity over a narrow frequency range.  In this case, the signal frequency could be yet higher, further improving noise performance. (Note that response slope would be less than 1, which could be compensated with a cascaded gain stage if desired).

Figure 9 plots RdB vs fmax/fs, where fmax/fs = ωmax/(2π).  Given the existence of the ωmax3 term in the denominator of Equation 16, it is not surprising that performance improves with increasing fmax.  Keep in mind though, that if the signal spectrum occupies less than this full bandwidth, you will not obtain the calculated noise performance.  As shown, the central difference differentiator has the worst performance, although we should note that the first difference differentiator would perform poorly vs. high frequency noise or interference, given its high-pass magnitude response.


Notes:

1.  ωmax = 2π fmax/fs

2.  RdB = 10*log10(R)

3.  Maximum response error = 2%, except for Lyons, which has 2.5% error at low frequencies.


Figure 8.  Magnitude responses of simple differentiators.


Figure 9.  Differentiator Noise Power Ratio RdB vs. fmax/fs for simple differentiators.

Maximum response error = 2%.


Custom Differentiators

The simple differentiators we discussed above are often appropriate, but they do not perform well in the presence of noise, and their usable bandwidth is arbitrary.  We can design better-performing differentiators using the following equation, which lets us choose a goal for ωc [1]:

$$ b(k)= \frac{\omega_c cos(\omega_c k)}{\pi k} - \frac{sin(\omega_c k)}{\pi k^2} \qquad (17)  $$

where -N/2 <= k <= N/2 and N = ntaps – 1.  The Matlab function diff_synth.m listed in Appendix A of the PDF version uses this equation, along with a Hanning window function, to design custom differentiators.  Note that the differentiators’ usable bandwidth ωmax is less than ωc.

As an example, suppose we want to design a 25-tap differentiator with ωc = π/2, as shown in the top of Figure 5.  The Matlab code is:

    ntaps= 25;
    wc= pi/2;              % radians cutoff frequency goal
    b= diff_synth(ntaps,wc,2*pi);

The resulting coefficients are shown in Figure 5 (bottom).  The magnitude response and response percent error (Equation 11) are plotted in Figure 6.  They are computed as follows:

    [h,w]= freqz(b,1,256);
    H= abs(h);              % magnitude response
    e_pct= 100*(H-w)./w;    % percent error

We can use diff_synth.m to design several custom differentiators with different lengths and different values of ωmax (or equivalently, fmax).  Table 2 lists parameters of five differentiators.  RdB is plotted vs fmax/fs for these differentiators in Figure 10.  Comparing Figure 10 to Figure 9, we see that the custom differentiators provide lower RdB, at the expense of more taps.  Also, as you might expect, to obtain a given value of RdB, fewer taps are required for differentiators with higher fmax/fs vs. lower fmax/fs.

Going back to our first example (Figures 1 - 3), we compared a central difference differentiator (R = 110, Table 1) and a custom “41 tap” differentiator (R = 3.05, Table 2), both having the same ωmax of 0.111π.  The difference in noise output power is 10*log10(110/3.05) = 15.6 dB.  That is, the central-difference differentiator has 15.6 dB higher noise power output (15.6 dB worse SNR), given Gaussian noise at the input. 

Finally, another option for designing differentiators is the Matlab function firpm, with ‘ftype’ = ‘differentiator’.  This function uses the Parks-McClellan algorithm, and should be more efficient than the method used in diff_synth.m.

Notes:

1.  ωmax = 2π fmax/fs

2.  RdB = 10*log10(R)

3.  Maximum response error = 2%

Figure 10.  Differentiator Noise Power Ratio RdB vs. fmax/fs for some custom

differentiators.  Maximum response error = 2%.


For Appendices, see the PDF version of this post.


References

1.  Lyons, Richard, Understanding Digital Signal Processing, 3rd Ed., Prentice Hall, 2011, section 7.1.3.

2.  Lyons, section 7.1.2.

3.  Lyons, p 810.

4.  Mitra, Sanjit K., Digital Signal Processing, 2nd Ed., McGraw-Hill, 2001, p 141.

5.  Spiegel, Murray R., et. al., Shaum’s Outlines Mathematical Handbook of Formulas and Tables, McGraw-Hill, 2018, equation 21.10.


                        Neil Robertson      March, 2022        Revised 3/29/22, 4/17/22



[ - ]
Comment by yatesMay 11, 2022

Hi Neil,

Thank you for a thought-provoking article. 

Can you give a mathematical expression involving the magnitude responses in Figure 1 that shows these two transfer functions have "equal usable bandwidth?" It sure doesn't _look_ like it to me, but I'm not sure how you define "usable bandwidth."

--Randy

[ - ]
Comment by neiroberMay 12, 2022

Hi Randy,

I base usable bandwidth on the percent error of the responses (I allowed 2% maximum error in my examples). In figure 1, each response has 2% error at slightly less than pi/8 radians.  

An example of the error is given in the center plot of figure 6.  I don't explicitly evaluate H(k):  I use the freqz function to compute it.  Later in the article, I give the following code to compute percent error:

[h,w]= freqz(b,1,256);

H= abs(h);             % magnitude response

e_pct= 100*(H-w)./w;   % percent error

I'm assuming here that I must meet the error critereon over the entire band from 0 to w_max.  For a narrowband signal, you only need to evaluate the error (with respect to linear) over the bandwidth of the signal.  In that case, the maximum usable center frequency of the differentiator can be higher.

regards,

Neil



To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.

Sign up

I agree with the terms of use and privacy policy.

Try our occasional but popular newsletter. VERY easy to unsubscribe.
or Sign in