Fractional Butterworth Filters: Optimal Approximations and Practical Challenges
Table Of Contents
Introduction
In a previous article, we discussed the Butterworth transform, which stretches the frequency axis according to the substitution: $$w \leftarrow w^n$$
This method operates under a constraint: \(n\) must be a positive integer.
In this article, we aim to explore a new method that allows for fractional orders, meaning \(n\) can be a real number rather than just an integer. There is no perfect solution for this, only approximations of such filters exist. Indeed, the elementary components of a transfer function are poles (which contribute -6 dB/oct slopes) and zeros (which contribute +6 dB/oct slopes). When their individual magnitudes are summed, we are left with a combined slope of +/-6 dB/oct at best. Our goal is therefore to find an efficient implementation that minimizes the number of poles and zeros used to approximate a fractional Butterworth filter.
The purpose of this approach is to achieve a continuous slope parameter, as shown in the graph below, which can offer greater flexibility for the end user.
|
|
\(\omega^n\) Approximation: The Tilt Filter Magnitude
Working with magnitude provides greater flexibility than directly manipulating the transfer function. As we previously discussed in this article, this approach allows us to adjust the magnitude without worrying about the phase. Our goal, therefore, is to transform the theoretical magnitude of the Butterworth filter by approximating \(\omega^n\) with a function that ensures the resulting magnitude can be used to generate an associated transfer function.
To approximate \(\omega^n\), one intuitive approach is to use the magnitude of a tilt filter. In a log-log plot, a tilt appears as a simple linear slope that varies continuously, making it an excellent approximation of \(\omega^n\) over a given range.
A tilt filter is commonly implemented by alternating poles and zeros. The greater the spacing between poles and zeros, the steeper the slope. Also, increasing the density of poles and zeros improves the approximation accuracy. To ensure a continuous transition of poles and zeros while \(n\) changes, we gradually shift the zeros toward \(s = 0\) in the complex plane, where they accumulate as the slope increases.
It is also worth noting that when \(n\) is an integer, the poles and zeros cancel out, resulting in exactly \(\omega^n\).
This approximation is given by:
$$ \boxed{ \omega^{2n} \approx \omega^{2n_{\text{int}}} \prod_{k=0}^{m-1} \frac{\left(1 + p_k^2\right) \left(\omega^2 + z_k^2\right)}{\left(1 + z_k^2\right) \left(\omega^2 + p_k^2\right)} } $$
where \(p_k\) and \(z_k\) are respectivly the poles and zeros of the corresponding tilt filter transfert function:
$$ \begin{aligned} z_k & = -2^{d\left(k - \frac{2m+1}{4} - n_{\text{frac}}\right)} \newline p_k & = -2^{d\left(k - \frac{2m+1}{4}\right)} \end{aligned} $$
for: $$k \in {0,1, \dots, m-1}$$
with:
- \(n \in \mathbb{R}\) → The order of the simulated tilt filter
- \(n_{\text{int}} = \lfloor n \rfloor\) → The integer part of n
- \(n_{\text{frac}} = n - n_{\text{int}}\) → The fractional part of n
- \(m \in \mathbb{N}^+\) → The true order of the approximated tilt filter
- \(d\) → A scaling factor that controls the density of \(z_k\) \(Z_{k}\), defining the precision of the approximation.
In order to simplify further expressions, we define:
$$ \begin{aligned} Z_{n}\left(\omega\right)=\omega^{2n_{int}}&\prod_{k=0}^{m-1}\frac{w^{2}+z_{k}^{2}}{1+z_{k}^{2}} \newline P\left(\omega\right)=&\prod_{k=0}^{m-1}\frac{1+p_{k}^{2}}{w^{2}+p_{k}^{2}} \end{aligned} $$
Such the approximation can be written as:
$$ \boxed{ \omega^{2n}\approx\frac{Z_{n}\left(\omega\right)}{P\left(\omega\right)} } $$
As you can see, \(p_k\), and therefore \(P(\omega)\), is not a function of n. This is a deliberate choice that enables certain factorizations, helping to reduce the resulting order when applying the fractional Butterworth transform to second-order prototype filters, as we will see in the next section.
The Python code below illustrates the approximation of \(\omega^n\) in a log-log space. You can think of it as a smoothed staircase, where the poles and zeros represent each step.
|
|
Injecting the Approximation into Prototype Filters
The next step is to integrate the approximation into first and second order low pass and high shelf filters. This is achieved using the following substitutions: \(\omega^{2n}\leftarrow\frac{Z_{n}\left(\omega\right)}{P\left(\omega\right)}\) and \(\omega^{4n}\leftarrow\frac{Z_{2n}\left(\omega\right)}{P\left(\omega\right)}\):
For a 1st-order low pass filter, we obtain:
$$ G_{LP1}\left(\omega^{n}\right)^{2}=\frac{1}{\omega^{2n}+1}\approx\frac{P\left(\omega\right)}{P\left(\omega\right)+Z_{n}\left(\omega\right)} $$
For a 2nd-order low pass filter, we obtain:
$$ G_{LP2}\left(\omega^{n},Q\right)^{2}=\frac{1}{\omega^{4n}+\left(\frac{1}{Q^{2}}-2\right)\omega^{2n}+1}\approx\frac{P\left(\omega\right)}{Z_{2n}\left(\omega\right)+\left(\frac{1}{Q^{2}}-2\right)Z_{n}\left(\omega\right)+P\left(\omega\right)} $$
For a 1st-order high shelf filter, we obtain:
$$ G_{HS1}\left(\omega^{n},g\right)^{2}=\frac{g\omega^{2n}+1}{\frac{1}{g}\omega^{2n}+1}\approx\frac{gZ_{n}\left(\omega\right)+P\left(\omega\right)}{\frac{1}{g}Z_{n}\left(\omega\right)+P\left(\omega\right)} $$
For a 2nd-order high shelf filter, we obtain:
$$ G_{HS2}\left(\omega^{n},g,Q_{up},Q_{down}\right)^{2}=\frac{g\omega^{4n}+\sqrt{g}\left(\frac{1}{Q_{down}^{2}}-2\right)\omega^{2n}+1}{\frac{1}{g}\omega^{4n}+\frac{1}{\sqrt{g}}\left(\frac{1}{Q_{up}^{2}}-2\right)\omega^{2n}+1}\approx\frac{gZ_{2n}\left(\omega\right)+\sqrt{g}\left(\frac{1}{Q_{down}^{2}}-2\right)Z_{n}\left(\omega\right)+P\left(\omega\right)}{\frac{1}{g}Z_{2n}\left(\omega\right)+\frac{1}{\sqrt{g}}\left(\frac{1}{Q_{up}^{2}}-2\right)Z_{n}\left(\omega\right)+P\left(\omega\right)} $$
The following code demonstrates these equations by comparing the desired magnitude with its approximation for each of the four filter prototypes.
|
|
Recovering the Transfer Function
To determine the transfer function from the magnitude \(G(\omega)\), we need to find the poles and zeros of \(G(\omega)\). This requires a root solver, as the order of the polynomial is typically too high, it’s known there is no closed-form solution for polynomials of degree \(\geq 5\).
To compute only half of the poles of \(G(\omega)\) needed to construct the transfer function of our desired prototype filter, we solve polynomials in terms of \(\omega^2\), reducing the computational complexity.
The code below demonstrates this polynomial resolution using np.roots(), after what we compute the transfer function of all common filter types. Poles and zeros that are too close and cancel each other out are removed from the plot for clarity.
|
|
The Problem of Sorting and Pairing Poles/Zeros
So, everything seems great, right? Well, not quite. Poles and zeros have no inherent ordering, which is a fundamental issue inherited from mathematics. For instance, if we vary \(n\) and \(Q\) and then return to their initial values (essentially completing a loop path), the order of the poles and zeros may shuffle. This issue does not arise for 1st-order low pass prototypes, as they depend on a single parameter \(n\), ensuring a single path.
The second challenge is root pairing. Even if we manage to order the poles and zeros correctly, their pairings may change, just as they do with the Bandpass Transform.
I haven’t explored a full solution yet, but I suspect that the best approach would involve strategically placing the poles of the tilt approximation, based on the cutoff frequencies of the decomposed prototype filters. This would require spreading the approximation’s poles and zeros dynamically, introducing new poles/zeros pairs that initially cancel each other out.
For the Bandpass Transform, we know that the root pairing problem makes it impossible to maintain continuous poles and zeros without reordering the pairs. This issue could potentially be addressed by using fourth-order filter sections instead of biquads.
Additionally, we need to be cautious about poles with very low cutoff frequencies (especially in high-pass filters), as they may lead to numerical error propagation.
Special Note on the Bandpass Transformation
Initially, I was confident that we could factorize the approximation similarly to second-order filters by using a shared denominator \(P(\omega)\). However, there doesn’t seem to be a straightforward solution. The main issue is that the bandpass transformation is applied after the Butterworth transform. When examining how poles and zeros are positioned in a BP-transformed filter, it becomes clear that it is difficult (if not impossible) to maintain a single, continuous path where the zeros slide smoothly as \(n\) changes on the complex s-plane. So in the end, we can simply apply the Bandpass Transform for bandpass, notch and peak filters, which unfortunatly doubles the filter order.
So, What Do We Get?
This approach works, but pole/zero continuity is handled only for 1-st order lowpass prototype filter. As mentioned earlier, I suspect there are viable solutions for the other prototypes, but more research is needed. Nevertheless, this approximation method provides excellent results in terms of CPU, as the \(\omega^n\) approximation is already optimal in term of poles/zeros usage. For our usage, we successfully implemented true fractional Butterworth low pass, high pass, band pass, and notch filters using just four extra biquad filters, in C++.