2
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ i haven't done this before except for specified LPF/HPF/BPF, but not for the general biquad. however, if you want to tackle that problem, check out the equation at the bottom of this answer, compute the derivative and set it to zero. there will likely be more than one value of $\omega$ (or $\phi$) that does that. in your case $a_0=1$. for the completely general problem, i don't think there is an easier way. please publish the result when you are done. $\endgroup$ – robert bristow-johnson Jun 2 '16 at 16:41
  • $\begingroup$ @robertbristow-johnson : Comment deleted. $\endgroup$ – Peter K. Jun 2 '16 at 20:36
  • $\begingroup$ @robertbristow-johnson, yeah I thought as much, hopefully a lot of terms will cancel. I'll probably do it differently to that equation you posted. If an expression in $z$ ends up being a quartic or less I'll post the solution, otherwise I doubt it will be useful to anyone. $\endgroup$ – keith Jun 2 '16 at 21:26
  • $\begingroup$ keith, you won't succeed doing it much differently. $\endgroup$ – robert bristow-johnson Jun 2 '16 at 23:58
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You're right, expanding $H(z)\bar{H(z)}$ directly in terms $z$ is not pretty, this method is more concise. For anyone else the explicit coefficients of the quadratic polynomial problem $p_0 + p_1 \phi + p_2 \phi^2 = 0$ are: $p_0 = {A_0 B_1 - A_1 B_0 \over 2}$, $p_1 = A_0 B_2 - A_2 B_0$, $p_2 = {A_1 B_2 - A_2 B_1 \over 2}$. $\endgroup$ – keith Jun 3 '16 at 11:43
  • $\begingroup$ i haven't figgered out what you are doing with the $p_n$ coefficients. i don't recognize the relationship with the $A_n$ and $B_n$ coefficients at all. $\endgroup$ – robert bristow-johnson Jun 3 '16 at 15:15
  • $\begingroup$ Your expression for $\phi$ is the solution to the quadratic roots of the polynomial in $p_n$ (some people use $a$, $b$, $c$ as per wiki article). The expression in that form when run on a PC can cause catastrophic cancellation. Out of habit I tend to pass the polynomial coefficients to a quadratic root solver that deals with underflow/overflow and avoids catastrophic cancellation, so I decoded your expression for $\phi$ into the coefficients. I assumed others would find it useful. $\endgroup$ – keith Jun 3 '16 at 15:54
  • $\begingroup$ i get it now ("duh"). my quadratic was a factor of two bigger (and i normally normalize the leading coefficient, "$p_2$"), which doesn't change the roots at all. rather than decoding the $\phi$ expression (that sounds really hard), i would think that just applying the quotient rule to $f'(\phi)$ would be easier. (and the $\phi^3$ term does drop out.) $\endgroup$ – robert bristow-johnson Jun 3 '16 at 16:14
  • 1
    $\begingroup$ Just a remark: this formula won't give you maxima of $|H(e^{j\omega})|$ at $\omega=0$ or $\omega=\pi$. So any maxima obtained from the quadratic equation must be compared to the frequency response values at DC and at Nyquist. $\endgroup$ – Matt L. Jun 4 '16 at 19:58
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ PS—There may be a more elegant solution in z-plane, which I didn't have time to consider. The main point that I started out to convey was that the maximum would only be at one of three places. Then it grew into showing that... $\endgroup$ – earlevel Dec 12 '16 at 10:08

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.