Is there a method to calculate the maximum gain of a biquad?

Question

Given an IIR expressed as a biquad:

$$ H(z) = {b_0 + b_1 z^{-1} + b_2 z^{-2} \over 1 + a_1 z^{-1} + a_2 z^{-2}} $$

Is there a way to calculate the maximum gain the biquad attains over the range $0$ and $\pi$?

One way is to calculate the turning points of $H(z) \bar H(z)$ and plug those back into $H(z) \bar H(z)$ but that would result in a horrible equation to solve in large powers of $z$ (probably up to order $7$) or differentiation by parts of a quotient based on terms of $\cos$ (depending on the two ways to tackle $H(z) \bar H(z) $).

Has this problem been simplified down to a ready made formula? Or is there an easier way?

Answer 1

in this previous answer, the squared-magnitude of a biquad frequency response is shown to be:

$$\begin{align} f(\phi) \triangleq |H(e^{j\omega})|^2 & = \left| \frac{b_0 + b_1 e^{-j\omega} + b_2 e^{-j2\omega}}{a_0 + a_1 e^{-j\omega} + a_2 e^{-j2\omega}} \right|^2 \\ & = \frac{B_0 + B_1 \phi + B_2 \phi^2}{A_0 + A_1 \phi + A_2 \phi^2} \\ \end{align}$$

where $ \phi \triangleq \sin^2\left(\frac{\omega}{2} \right) $

and

$\begin{align} B_0 & = \left(\frac{b_0+b_1+b_2}{2}\right)^2 \\ B_1 & = - (4b_0b_2 + b_1(b_0+b_2)) \\ B_2 & = 4b_0b_2 \\ A_0 & = \left(\frac{a_0+a_1+a_2}{2}\right)^2 \\ A_1 & = - (4a_0a_2 + a_1(a_0+a_2)) \\ A_2 & = 4a_0a_2 \\ \end{align}$

maxima and minima occur when the derivative is zero.

the derivative of the above is zero when

$$ \phi = \frac{A_2B_0 - A_0B_2 \pm \sqrt{A_0^2B_2^2 + A_2^2B_0^2 - 2A_0A_2B_0B_2 - A_0A_1B_1B_2 \\ + A_0A_2B_1^2 + A_1^2B_0B_2 - A_1A_2B_0B_1}}{A_1B_2 - A_2B_1} $$

then evaluating $f(\phi)$ at those two values of $\phi$ is another laborious pain in ass. and it won't be my ass.

Answer 2

Let me refer to the denominator coefficients as the pole coefficients, and the numerator coefficients as the zero coefficients. This is for those who visualize poles as able to peak magnitude, and zeros able to dip magnitude magnitude. For instance, if the biquad were all-pole (numerator = 1), the max magnitude could be found at the pole’s phase angle. Similarly the minimum for an all-zero (denominator = 1), the minimum would be at the zero’s phase angle. (Complex conjugate poles and zeros for a real filter, we only need look at one.)

Having both poles and zeros spoils that simplicity, as a nearby zero can pull down the pole response and allow it to be higher elsewhere. However, if we don’t find the maximum at the pole angle, it will be at 0 or pi. Similarly, the minimum magnitude will be the minimum of the responses at the zero angle, 0, or pi.

Here I’ll present some Python code, written for clarity (I hope). Getting the pole (for max) and zero (if interested in min) locations is a matter of solving for the roots of the denominator and numerator, respectively. We’ll take only the “+” case here, since we have complex conjugates for real filters:

import math
import cmath

# get location via quadratic formula
def quadraticToZplane(a, b, c):
    return (-b + cmath.sqrt(b * b - 4 * a * c)) / (2 * a)

To evaluate the magnitude response at 0, pi, and the pole angle for max magnitude, or the zero angle for min magnitude, the evalMagZplane function takes zero and pole locations, the location on the unit circle to evaluate, and a gain factor. Magnitude response in the z plane is the product of the zero distances to the evaluation point, divided by the product of the distances from the poles to the evaluation point, multiplied by the filter gain (the absolute value of a0 divided by b0).

def evalMagZplane(zeroLoc, poleLoc, unitLoc, gain):
    magNumer  = math.hypot(unitLoc.real - zeroLoc.real, unitLoc.imag - zeroLoc.imag)
    magNumer *= math.hypot(unitLoc.real - zeroLoc.real, unitLoc.imag + zeroLoc.imag)
    magDenom  = math.hypot(unitLoc.real - poleLoc.real, unitLoc.imag - poleLoc.imag)
    magDenom *= math.hypot(unitLoc.real - poleLoc.real, unitLoc.imag + poleLoc.imag)
    return magNumer / magDenom * gain

Here’s our solution to max and min, using the others; math.phase is equivalent to atan2(<imaginaryPart>, <realPart>):

def biquadMaxMin(zeros, poles):
    zeroLoc = quadraticToZplane(zeros[0], zeros[1], zeros[2])
    poleLoc = quadraticToZplane(poles[0], poles[1], poles[2])

    zeroAngle = cmath.phase(zeroLoc)
    poleAngle = cmath.phase(poleLoc)

    zeroUnitLoc = cmath.rect(1, zeroAngle)
    poleUnitLoc = cmath.rect(1, poleAngle)

    gain = abs(zeros[0] / poles[0])
    mag0 = evalMagZplane(zeroLoc, poleLoc, 1, gain)
    magPi = evalMagZplane(zeroLoc, poleLoc, -1, gain)
    magZero = evalMagZplane(zeroLoc, poleLoc, zeroUnitLoc, gain)
    magPole = evalMagZplane(zeroLoc, poleLoc, poleUnitLoc, gain)

    # max is at the pole angle, or 0 or pi
    magMax = max(mag0, magPole, magPi)
    # min is at the zero angle, or 0 or pi
    magMin = min(mag0, magZero, magPi)

    return (magMax, magMin)

A simple test, with biquad coefficients for a peaking filter set to +7 dB gain. So, max of 7 dB is expected, min of 0 dB (converting to dB in the print statement, rounding to make it prettier):

# +7 dB peaking filter
zeros = [ 1.0473669305387907, -1.4557390814245152, 0.8761559204661808 ]
poles = [ 1.0, -1.4557390814245152, 0.9235228510049719 ]

mm = biquadMaxMin(zeros, poles)
print 'max = {} dB, min = {} dB'.format(round(20 * math.log10(mm[0]), 2), round(20 * math.log10(mm[1]), 2))

max = 7.0 dB, min = -0.0 dB

I was improving some old plotting code on my website, to make it so that the maximum is the true max and not just the max of points I chose to plot. I found this page, it didn’t answer my question, so I thought about it a few minutes and came up with this. I expect that examination of the real zero and pole locations would indicate whether to solve for the pole angle, 0, or pi, but I'm not sure offhand whether it would be a computational improvement over solving all three and taking the maximum. I won’t be available till after the holidays, so I’ll leave that as an exercise for anyone interested.