apulSoft Blog

Jun 1st, 2024 - filters math dsp

Matched Digital Allpass IIR Filters

This is inspired by the papers of Martin Vicanek, who created new versions of the typical digital IIR filters that match the magnitude response of analog prototypes better than the commonly used bilinear transform method. Instead of transforming poles and zeros, his methods use mathematical constraints to construct digital filters that match the ideal analog magnitude curves.

Matching magnitudes for allpass filters is obviously pointless. By definition they have magnitude 1 at all frequencies, regardless of being calculated in analog or digital space. This means phase values have to be matched instead. It is best done in complex space to avoid converting the phase to angles using atan2(), which is problematic to work with mathematically.

As allpass filters have magnitude 1 everywhere, their complex frequency response always lies on the complex unit circle $e^(-ix)$. Some relations between phase angles and complex locations $c$:

	Phase	c
	0	1 + i * 0
	-pi/2	0 + i * (-1)
	-pi	-1 + i * 0

1-Pole Case

A 1-pole digital allpass filter has the complex frequency response:

$H_1(z)=(c_0+z^-1)/(1+c_0*z^-1)$ with $z = e^(-iω)$ and $ω = 2pif$

There is only one variable, so only one phase point can be matched. The only reasonable choice is asking for phase $-pi//2$ at the center frequency $f_c$ like the analog version has.

$w_c = 2pif_c$, $H_1(ω_c) = -i$

The free algebra application Maxima can easily solve this for $c_0$. As $c_0$ is known to be real, the imaginary part can be discarded.

trigsimp(realpart(solve((c_0+%e^(-%i*ω_c))/(1+c_0*%e^(-%i*ω_c))=-%i, c_0)));

$c_0 = (-cos(ω_c))/(sin(ω_c)+1)$

Unsurprisingly, this is equivalent to any other reasonable method to create a 1-pole digital allpass IIR filter.

2-Pole Case

The 2-pole case is much more interesting, as the quality factor $Q$ comes into play, this gives two parameters to work with and two complex frequency response values can be matched. To keep the math manageable, the points chosen are the center frequency (phase $-pi$) and a new point at a lower frequency where the phase is $-pi//2$.
The first task is finding this lower frequency on an analog prototype filter.

For simpler terms, the damping factor $ζ$(zeta) is used instead of $Q$.

$ζ = 1/(2Q)$

An analog allpass IIR prototype filter with center frequency 1 has the complex frequency response:

$H_2(s)=(s^2-2ζs+1)/(s^2+2ζs+1)$ with $s=i*f$

$H_2(f)=(1-f^2+i(-2ζf))/(1-f^2+i*2ζf)$

This can be solved for the $-pi//2$ point using Maxima again:

$H_2(f)=-i$

solve((1-f^2+%i*(-2*ζ*f))/(1-f^2+%i*(2*ζ*f))=%i,f);

$f_h=-zeta+-sqrt(ζ^2+1)$

To get $0 < f_h < 1$, one solution is selected: $f_h=-ζ+sqrt(ζ^2+1)$

With this result, the digital IIR filter can now be matched. It has the complex frequency response:

$H_2(z)=(c_0+c_1*z^-1+z^-2)/(1+c_1*z^-1+c_0*z^-2)$

Two matching points are defined:

Center frequency: $w_c = 2pif_c$, phase =$-pi$ => $H_2(ω_c) = -1$

-Half-pi frequency: $w_h = 2pif_cf_h$, phase =$-pi/2$ => $H_2(ω_h) = -i$

Using Maxima, this 2-variable system can be solved.

realpart(solve([(c_0+c_1*%e^(-%i*ω_c)+%e^(-2*%i*ω_c))/(1+c_1*%e^(-%i*ω_c)+c_0*%e^(-2*%i*ω_c))=-1,(c_0+c_1*%e^(-%i*ω_h)+%e^(-2*%i*ω_h))/(1+c_1*%e^(-%i*ω_h)+c_0*%e^(-2*%i*ω_h))=-%i
        ], [c_0, c_1]));

This produces very extensive results, but repeated trigsimp(trigexpand(%)) and some extra help from WolframAlpha finally lead to usable results for $c_0$ and $c_1$.

$c_0 = (2sin(ω_h))/(sin(ω_h)+cos(ω_h)-cos(ω_c)) - 1$

$c_1 = (-2cos(ω_c)sin(ω_h))/(sin(ω_h)+cos(ω_h)-cos(ω_c))$

Improved 1-Pole Matching

As can be seen when adjusting the interactive curves, the $-pi//2$ phase angle is not ideal as a matching point. Especially high-$f_c$/high-$Q$ situations miss the target curve in a way that looks like it could be improved by matching at a different phase angle.

$H_1(s) = H_1(f_a) = (1-s)/(1+s) = (1 - i*f_a)/(1 + i*f_a) = e^(-i*a)$

Looking at the real part only:

$(1-f_a^2)/(1+f_a^2)=cos(a)$

$f_a = +-sqrt((1 - cos(a))/(1 + cos(a)))$

Determining the digital coefficient:

$omega_a = 2pif_af_c$

$H_1(z) = H_1(w_a) = (c_0 + z^-1)/(1 + c_0 * z^-1) = (c_0 + e^(-i*omega_a))/(1 + c_0 * e^(-i*omega_a)) = e^(-i*a)$

$c_0=(cos omega_a - cos a)/(cos(a+omega_a)-1)$

For the one-pole matching, I chose $a = pi//4$ to get an extra octave.

Improved 2-Pole Matching

For the two pole case, the new goal is to match at two frequencies: The center frequency where the phase is $-pi$, and a flexible frequency $f_a$ where the phase angle is $a$.

To determine $f_a$ from $a$, this analog prototype equation can be solved.

$H_2(s) = H_2(s) = H_2(f_a) = (1-f_a^2+i(-2\zetaf_a))/(1-f_a^2+i*2\zetaf_a)=e^(-i*a)$

If $\zeta, f$ and $a$ are real, it is equivalent to:

$(f^4 -2 (2 \zeta^2 + 1) f^2 + 1)/(f^4 + 2(2\zeta^2 - 1) f^2 + 1) + i*(4 \zeta f (f^2 - 1))/(f^4 + 2(2\zeta^2 - 1) f^2 + 1) = cos(a) - i * sin(a)$

Restricting the possible results to $0 <= f_a < pi//2$, setting the real parts equal leads to:

$R = (cos a * (\zeta^2 - 1) + 2) * cos a$

$B = -2 * \zeta^2 * cos a + 2 * z * sqrt(R) + cos a - 2$

$f_a = sqrt(B / (cos a - 1))$

Now the digital complex frequency response can be matched at the two points.

$H(z) = H(w_c) = -1$

$omega_a = 2pif_af_c$

$H(z) = H(omega_a) = e^(-i*a)$

This can be solved using Maxima.

solve([(c_0+c_1*%e^(-%i*c)+%e^(-2*%i*w_c))/(1+c_1*%e^(-%i*c)+c_0*%e^(-2*%i*c))=-1,(c_0+c_1*%e^(-%i*w_a)+%e^(-2*%i*w_a))/(1+c_1*%e^(-%i*w_a)+c_0*%e^(-2*%i*w_a))=%e^(-%i*a)], [c_0, c_1]);

The page-filling result can be simplified to:

$C = sin\frac{a}{2} * (cos w_a - cos w_c),\ D = sin w_a * cos\frac{a}{2}$

$c_0 = (D - C) / (C + D)$

$c_1 = (-2D*cos w_c) / (C + D)$

One question remains: How to choose the phase angle $a$ to target. As seen in the first attempt, $pi//2$ is too high. Very small values tend to produce numerical problems as $f_a$ becomes tiny. Through experimentation I arrived at:

$a = sqrt(zeta/(2f_c))$, clamped to $0.01<=a<=1.5$ to stay out of numerical trouble.

Interactive Demo

The plot shows the phase response of four allpass filters at 44.1kHz samplerate.

Blue: ideal analog curve
Red: bilinear transform to digital
Green: digital matching at fc and at pi/2 phase.
Orange: digital matching at fc and at phase a.

One Pole (green and red are equal):

Two Poles:

Interactive plots created using Function Plot.

C++ Implementation

The results have been implemented in c++. Higher order allpass filters can be split into 1-pole and 2-pole sections. The _r2 functions are the improved versions.

/*
  ap_allpass.h Copyright (c) 2024 Adrian Pflugshaupt. 
  Blog-article: https://apulsoft.ch/blog/matched-allpass/

  This work is licensed under the terms of the MIT license.  
  For a copy, see <https://opensource.org/licenses/MIT>. 
*/

#pragma once
#include <algorithm> // std::min
#include <cmath>

namespace ap_allpass {

///        b0 + b1 * z^-1 + b2 * z^-2
/// H(z) = --------------------------
///         1 + a1 * z^-1 + a2 * z^-2
///
struct BiquadDigitalCF {
    double a1 = 0, a2 = 0;
    double b0 = 1, b1 = 0, b2 = 0;
};

/// Create digital IIR allpass coefficients for a 1-pole allpass 
/// with normalized center freq f_c (0.5 = nyquist) 
inline BiquadDigitalCF createDigitalAllpass1P(double f_c) {
    auto w_c = 2 * M_PI * std::min(0.499, f_c);

    auto c_0 = -std::cos(w_c) / (std::sin(w_c) + 1);

    BiquadDigitalCF res;
    res.b0 = c_0;
    res.b1 = 1;
    res.b2 = 0;
    res.a1 = c_0;
    res.a2 = 0;
    return res;
}

/// Create digital IIR allpass coefficients for a 2-pole allpass 
/// with normalized center freq f_c (0.5 = nyquist) and quality factor Q
inline BiquadDigitalCF createDigitalAllpass2P(double f_c, double Q) {
    auto w_c = 2 * M_PI * std::min(0.499, f_c);

    // figure out -pi/2 phase point f_h from analog prototype
    auto zeta = 1 / (2 * Q);
    auto f_h = -zeta + std::sqrt(zeta * zeta + 1);

    // create digital allpass through two points
    auto w_h = 2 * M_PI * std::min(0.498, f_c * f_h);
    auto cos_w_c = std::cos(w_c);
    auto sin_w_h = std::sin(w_h), cos_w_h = std::cos(w_h);

    auto bot = 1 / (sin_w_h + cos_w_h - cos_w_c);
    auto c_0 = (2 * sin_w_h) * bot - 1;
    auto c_1 = (-2 * cos_w_c * sin_w_h) * bot;

    BiquadDigitalCF res;
    res.b0 = c_0;
    res.b1 = c_1;
    res.b2 = 1;
    res.a1 = c_1;
    res.a2 = c_0;
    return res;
}

/// Create digital IIR allpass coefficients for a 1-pole allpass 
/// with normalized center freq f_c (0.5 = nyquist) 
inline BiquadDigitalCF createDigitalAllpass1P_r2(double f_c) {
    const auto a = M_PI / 4; // always fit at pi/4 phase value
    const auto cos_a = std::cos(a);
    const auto f_a = std::sqrt((1 - cos_a) / (1 + cos_a));

    auto w_a = 2 * M_PI * std::min(0.499, f_c * f_a);
    auto c_0 = (std::cos(w_a) - cos_a) / (std::cos(a + w_a) - 1);

    res.b0 = c_0;
    res.b1 = 1;
    res.b2 = 0;
    res.a1 = c_0;
    res.a2 = 0;
    return res;
}

/// Create digital IIR allpass coefficients for a 2-pole allpass 
/// with normalized center freq f_c (0.5 = nyquist) and quality factor Q
inline BiquadDigitalCF createDigitalAllpass2P_r2(double f_c, double Q) {
    auto w_c = 2 * M_PI * std::min(0.4995, f_c);

    // zeta damping factor
    auto zeta = 1 / (2 * Q); 
    auto zetasq = zeta * zeta;

    // match phases at f_c and the point where the phase is -a.
    auto a = std::clamp(std::sqrt(zeta / (2 * f_c)), 0.01, 1.5);
    auto cos_a = std::cos(a);
    auto A = cos_a + 1;
    auto sin_a_h = std::sqrt(0.5 - 0.5 * cos_a);
    auto cos_a_h = std::sqrt(0.5 + 0.5 * cos_a);

    auto R = (A * (zetasq - 1) + 2) * A;
    auto B = -2 * zetasq * A + 2 * zeta * std::sqrt(R) + A - 2;
    auto f_a = std::sqrt(B / (cos_a - 1));

    // create digital allpass through two points
    auto w_a = 2 * M_PI * std::min(0.499, f_c * f_a);

    auto cos_w_c = std::cos(w_c);
    auto sin_w_a = std::sin(w_a), cos_w_a = std::cos(w_a);

    auto C = sin_a_h * (cos_w_a - cos_w_c);
    auto D = sin_w_a * cos_a_h;

    auto bot = -1 / (C + D);
    auto c_0 = (C - D) * bot;
    auto c_1 = 2 * cos_w_c * D * bot;

    res.b0 = c_0;
    res.b1 = c_1;
    res.b2 = 1;
    res.a1 = c_1;
    res.a2 = c_0;
    return res;
}

} // namespace ap_allpass