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.
[z, p, k]
and constructing a filter from there, rather than using the more direct filtering code with[b, a]
. $\endgroup$[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$