4
$\begingroup$

Say I want to apply a zero-phase Butterworth bandpass filter, implemented in MATLAB using butter and filtfilt. So, for an order-$4n$ filter with passband $[f_{\mathrm{min}\,},f_{\mathrm{max}}]$, I believe the transfer function $H$ is given by $$ H(2\pi if) \ = \ \frac{1}{1+\left( \frac{f}{f_\mathrm{max}-f_\mathrm{min}} + \frac{f^{-1}}{f_\mathrm{max}^{-1}-f_\mathrm{min}^{-1}} \right)^{\!2n}}\,. $$

It seems that for any given digital signal $x(t)$ (other than the constant zero signal), there is a critical number $N_x$ such that:

  • for any $n < N_x$ the results are fairly sensible;
  • but for $n \geq N_x$, the results change considerably and are unreliable.

(1) What is the cause of this behaviour? Is it due to numerical imprecisions that grow large when the stability of the filter is poor?

(2) Is there a straightforward way to determine what this critical number $N_x$ is?

Let me illustrate the phenomenon with a couple of examples:

EXAMPLE 1

Consider the signal $x(t)=\sin(2\pi t)$, $t \in \{0,0.01,0.02\ldots,99.99,100\}$. (This signal has 10001 points altogether.)

Suppose I apply the Butterworth filter with passband $[f_{\mathrm{min}\,},f_{\mathrm{max}}]=[0.5,3]$. Of course, I expect my filtered signal $x_{\mathrm{filt}}(t)$ to be basically the same as $x(t)$, except:

  • there may be edge effects near the start and/or end;
  • for $n=1$ or $2$, the amplitude is likely to be slightly reduced, as the frequency of the sinusoid $x(t)$ is not very close to the centre frequency $\sqrt{1.5}$ of the filter.

For each $n \leq 5$, the results are as expected. E.g., for $n=5$ (corresponding to filter order 20), the result is as follows:

However, when I now take $n=6$ (filter order 24), I suddenly get a drastically different result:

A wavelet transform of the original and filtered signal for $n=6$ looks as follows:

For even larger $n$, the results simply blow up. For $n=7$, the filtered signal is equal to NaN for about the first two thirds of the duration, and then starts taking values of the order of $\,-10^{305}$, going up to a maximum value of about $4 \times 10^{231}$.

So in this example, $N_x=6$.

EXAMPLE 2

This is exactly the same as Example 1, except with the sampling frequency considerably reduced.

Consider the signal $x(t)=\sin(2\pi t)$, $t \in \{0,0.08,0.16\ldots,99.92,100\}$. (This signal has 1251 points altogether.)

Once again, suppose I apply the Butterworth filter with passband $[f_{\mathrm{min}\,},f_{\mathrm{max}}]=[0.5,3]$. For $n \leq 16$, I get sensible results. Here is the filtered signal with $n=16$:

I think the "early warning sign" here is that the slight inaccuracy of the filtered signal persists throughout the whole duration, not just at the start and end; for $n \leq 15$, the visible inaccuracy is only at the edges.

For $n=17$, I get quite a bad result:

In this case, Fourier and wavelet transforms are as follows:

Once again, for $n=18$, the results blow up:

(Please note that the values of the filtered signal are extremely large throughout, not just for $t<35$. At $t=70$, the value is about $-6.5 \times 10^{11}$.)

So in this example, $N_x=17$.

REMARK: In many (probably most) other examples, there is no intermediary $n$-value between the $n$-values giving consistent results and the $n$-values for which the filtered signal blows up to extreme values. In this regard (but only this regard), the above two examples seem somewhat exceptional.

$\endgroup$
  • $\begingroup$ Okay, it seems a very similar question has been asked at dsp.stackexchange.com/questions/17980/…. The accepted answer there suggests that the problem can be overcome by using [z, p, k] and constructing a filter from there, rather than using the more direct filtering code with [b, a]. $\endgroup$ – Julian Newman Jan 29 '18 at 22:49
  • $\begingroup$ But it still perhaps remains an interesting question as to what the largest order is for which the [b, a] approach will give good results (i.e. I don't think my question is completely redundant). After all, the [b, a] approach is used (presumably because it is less effort). $\endgroup$ – Julian Newman Jan 29 '18 at 22:51
  • 4
    $\begingroup$ Usually we use cascades of second-order sections instead of doing it all in one section. $\endgroup$ – endolith Jan 29 '18 at 23:01
  • 1
    $\begingroup$ As described in en.wikipedia.org/wiki/Digital_biquad_filter, I guess. $\endgroup$ – Julian Newman Feb 3 '18 at 18:49
  • 1
    $\begingroup$ Yep, that's the standard in both digital and analog. $\endgroup$ – endolith Feb 4 '18 at 20:39
7
$\begingroup$

One cause is that higher order Butterworth filters have poles closer to the unit circle. This nearby infinite gain point increases the likelihood of numerical instabilities. (e.g. rounding/arithmetic/quantization noise may move a pole to the “wrong” side of the unit circle.)

Where the numerical noise will blow up depends on your executable code’s exact arithmetic sequence, and the quantization and number format used. Often, dividing up a high-order filter into a sequence of biquad stages produces a more stable sequence of arithmetic operations.

$\endgroup$
3
$\begingroup$

The cause is that in general, relating the roots of a polynomial to it’s coefficients is numerically ill conditioned. A well known example is Wilkerson’s polynomial.

https://en.m.wikipedia.org/wiki/Wilkinson%27s_polynomial

Small perturbations have relatively large effects.

The safe order also depends on numerical precision. Analog Butterworth filters are also difficult as order increases and the order one cascades biquad factors are important to avoid clipping a stage.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.