Biquad formulas

For fixed filters, we can plug biquad coefficients into our programs. But often, we need to calculate them on the fly, to user settings or changes in sample rate. As a companion to the biquad calculator, here are the formulas used, in JavaScript; refer to our article on biquads for implementation diagrams:

function calcBiquad(type, Fc, Fs, Q, peakGain) {
    var a0,a1,a2,b1,b2,norm;
    
    var V = Math.pow(10, Math.abs(peakGain) / 20);
    var K = Math.tan(Math.PI * Fc / Fs);
    switch (type) {
        case "lowpass":
            norm = 1 / (1 + K / Q + K * K);
            a0 = K * K * norm;
            a1 = 2 * a0;
            a2 = a0;
            b1 = 2 * (K * K - 1) * norm;
            b2 = (1 - K / Q + K * K) * norm;
            break;
        
        case "highpass":
            norm = 1 / (1 + K / Q + K * K);
            a0 = 1 * norm;
            a1 = -2 * a0;
            a2 = a0;
            b1 = 2 * (K * K - 1) * norm;
            b2 = (1 - K / Q + K * K) * norm;
            break;
        
        case "bandpass":
            norm = 1 / (1 + K / Q + K * K);
            a0 = K / Q * norm;
            a1 = 0;
            a2 = -a0;
            b1 = 2 * (K * K - 1) * norm;
            b2 = (1 - K / Q + K * K) * norm;
            break;
        
        case "notch":
            norm = 1 / (1 + K / Q + K * K);
            a0 = (1 + K * K) * norm;
            a1 = 2 * (K * K - 1) * norm;
            a2 = a0;
            b1 = a1;
            b2 = (1 - K / Q + K * K) * norm;
            break;
        
        case "peak":
            if (peakGain >= 0) {    // boost
                norm = 1 / (1 + 1/Q * K + K * K);
                a0 = (1 + V/Q * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - V/Q * K + K * K) * norm;
                b1 = a1;
                b2 = (1 - 1/Q * K + K * K) * norm;
            }
            else {    // cut
                norm = 1 / (1 + V/Q * K + K * K);
                a0 = (1 + 1/Q * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - 1/Q * K + K * K) * norm;
                b1 = a1;
                b2 = (1 - V/Q * K + K * K) * norm;
            }
            break;
        case "lowShelf":
            if (peakGain >= 0) {    // boost
                norm = 1 / (1 + Math.SQRT2 * K + K * K);
                a0 = (1 + Math.sqrt(2*V) * K + V * K * K) * norm;
                a1 = 2 * (V * K * K - 1) * norm;
                a2 = (1 - Math.sqrt(2*V) * K + V * K * K) * norm;
                b1 = 2 * (K * K - 1) * norm;
                b2 = (1 - Math.SQRT2 * K + K * K) * norm;
            }
            else {    // cut
                norm = 1 / (1 + Math.sqrt(2*V) * K + V * K * K);
                a0 = (1 + Math.SQRT2 * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - Math.SQRT2 * K + K * K) * norm;
                b1 = 2 * (V * K * K - 1) * norm;
                b2 = (1 - Math.sqrt(2*V) * K + V * K * K) * norm;
            }
            break;
        case "highShelf":
            if (peakGain >= 0) {    // boost
                norm = 1 / (1 + Math.SQRT2 * K + K * K);
                a0 = (V + Math.sqrt(2*V) * K + K * K) * norm;
                a1 = 2 * (K * K - V) * norm;
                a2 = (V - Math.sqrt(2*V) * K + K * K) * norm;
                b1 = 2 * (K * K - 1) * norm;
                b2 = (1 - Math.SQRT2 * K + K * K) * norm;
            }
            else {    // cut
                norm = 1 / (V + Math.sqrt(2*V) * K + K * K);
                a0 = (1 + Math.SQRT2 * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - Math.SQRT2 * K + K * K) * norm;
                b1 = 2 * (K * K - V) * norm;
                b2 = (V - Math.sqrt(2*V) * K + K * K) * norm;
            }
            break;
    }

    return [ a0, a1, a2, b1, b2 ];
}
This entry was posted in Biquads, Digital Audio, Filters, IIR Filters. Bookmark the permalink.

32 Responses to Biquad formulas

  1. Miles Egan says:

    Your entire series of posts on filters and resampling has been very helpful. Thanks!

  2. weili says:

    useful

  3. weili says:

    Hi Nigel,
    Do you have plan to write a tutorial how to calculate the FIR coefficient based on cut-off frequency for lowpass or bandpass filters?

    • Nigel Redmon says:

      I’ll keep that in mind. I’ve been busy—I have a couple of partially written articles on other subjects that I need to get back to and finish up, but I’ll think about an article on FIRs. Thanks for the feedback.

  4. Marika Picardal says:

    Hello Nigel,

    I am planning to implement a Biquad IIR filter on an FPGA platform. Right now, I’m searching for a simple yet useful Audio Application for the biquad filter. I would like to ask for your suggestions.

    Thanks,

    Marika

  5. Malcolm Rest says:

    Nigel hi,

    Thanks for the great site.

    What is the applicational difference between a notch and negative peak with hi-Q?

    Could you share/upload a calculation for Low/High Shelf EQ with “Q factor” for their sharpness?

    • Nigel Redmon says:

      Hi Malcolm,

      The application difference between a notch and a negative peak is that the notch, which comes to a very sharp point, is usually used to target a specific frequency; most often it’s used to remove hum from the AC power line. Power line frequency (60 Hz in the US, 50 Hz in the UK) is regulated precisely (if not, power grids would have trouble combining power from multiple generators), so surgical precision can be used to remove it and its harmonics. Moving notches are also used as an effect (as in a phase-shifter).

      Negative peak filters are more often used to de-emphasize a frequency area—perhaps a recording that was made in a room with a bad resonance.

    • Nigel Redmon says:

      Oh, and the Q for shelves…that’s a good idea, I’ll keep that in mind as I think about the next filter article…

      • Malcolm Rest says:

        Dear Nigel hi,

        Let me contribute to this great web site, not just “nag” you all the time ;-)

        I have calculated and tested the needed code change for “Q_Factor” for Low_Shelf & High_Shelf.

        It is done by simply dividing the second mathematical expression of each coefficient, besides a1 & b1, by Q.

        Here is an example based on your “Low Shelf”:

        case “lowShelf”:
        if (peakGain >= 0) { // boost
        norm = 1 / (1 + Math.SQRT2 * K /Q+ K * K);
        a0 = (1 + Math.sqrt(2*V) * K /Q+ V * K * K) * norm;
        a1 = 2 * (V * K * K – 1) * norm;
        a2 = (1 – Math.sqrt(2*V) * K /Q+ V * K * K) * norm;
        b1 = 2 * (K * K – 1) * norm;
        b2 = (1 – Math.SQRT2 * K /Q+ K * K) * norm;
        }
        else { // cut
        norm = 1 / (1 + Math.sqrt(2*V) * K/Q + V * K * K);
        a0 = (1 + Math.SQRT2 * K /Q+ K * K) * norm;
        a1 = 2 * (K * K – 1) * norm;
        a2 = (1 – Math.SQRT2 * K /Q+ K * K) * norm;
        b1 = 2 * (V * K * K – 1) * norm;
        b2 = (1 – Math.sqrt(2*V) * K /Q+ V * K * K) * norm;
        }

        Hope this is useful to some…

        Malcolm

        • Nigel Redmon says:

          Thanks Malcolm!

        • Nigel Redmon says:

          Just getting back around to this…After my recent C++ article, it just occurred to me that I didn’t work in Q into the shelves. I did, today, then remembered your comment…yours was different than mine, so I tried it out, and see a problem. Your version won’t work right for “flat” (Q = 0.7071) response.

          I did a quickie implementation, replacing all sqrt(2) with 1/Q. That should be correct, but I didn’t go through the BLT for a couple of reasons: first, lack of time, second, doing it this way (tweaking the pole/zero relationships on both sides of the shelf frequency) gives an overshoot on both sides. From what I’ve seen, the “good sounding” way to do it is to have the overshoot on the unity side of the shelf. So I need to look closer at the best “musical” implementation, but for now here’s the version with overshoot on both sides for Q > 0.7071 (in C++):

          case bq_type_lowshelf:  // with Q
          if (peakGain >= 0) {    // boost
            norm = 1 / (1 + 1/Q * K + K * K);
            a0 = (1 + sqrt(V)/Q * K + V * K * K) * norm;
            a1 = 2 * (V * K * K - 1) * norm;
            a2 = (1 - sqrt(V)/Q * K + V * K * K) * norm;
            b1 = 2 * (K * K - 1) * norm;
            b2 = (1 - 1/Q * K + K * K) * norm;
          }
          else {    // cut
            norm = 1 / (1 + sqrt(V)/Q * K + V * K * K);
            a0 = (1 + 1/Q * K + K * K) * norm;
            a1 = 2 * (K * K - 1) * norm;
            a2 = (1 - 1/Q * K + K * K) * norm;
            b1 = 2 * (V * K * K - 1) * norm;
            b2 = (1 - sqrt(V)/Q * K + V * K * K) * norm;
          }
          break;
          
  6. Malcolm rest says:

    Nigel hi,

    I believe there is a typo at “High Shelf/cut” calculation.
    “Norm” seems to appear twice.

    Thanks for the post !

    Malcolm

    • Nigel Redmon says:

      Hi Malcolm. As you noted later, yes, “norm” was recalculated. In fact, I grabbed the calculation from Udo Zölzer’s Digital Audio Signal Processing, essentially, for the biquad calculator code, and reproduced that code here. But you made me think…why should it be? The “norm” variable is short for “normalize”—it’s the b0 term by which all the other components are normalized. I decided to rearrange the b1 and b2 calculations to make them more uniform with the rest of the filters—thanks for mentioning it. I suppose that the lack of uniformity in Zölzer’s book was just a minor oversight—either form produces the same results. I’ll list the original interpretation here, and make the change (which is to multiply numerators and denominators for b1 and b2 by V) in the article:

      case "highShelf":
          ...
          else {    // cut
              norm = 1 / (V + Math.sqrt(2*V) * K + K * K);
              a0 = (1 + Math.SQRT2 * K + K * K) * norm;
              a1 = 2 * (K * K - 1) * norm;
              a2 = (1 - Math.SQRT2 * K + K * K) * norm;
              norm = 1 / (1 + Math.sqrt(2/V) * K + K * K / V);
              b1 = 2 * (K * K / V - 1) * norm;
              b2 = (1 - Math.sqrt(2/V) * K + K * K / V) * norm;
          }
      
  7. Malcolm rest says:

    My mistake – just double checked and there is a need for two different “norm” params in the code.

  8. Malcolm rest says:

    Dear Nigel,

    Thanks for your prompt answers, must appreciated.

    I was looking up some other sources for IIR EQ filtering and noticed some of them use trigonometric functions to shift from exp.(-jw) field for the calculation part itself.

    It looks like trigonometric calculations make the code much shorter. Is there a side effect to using trigonometric calculations?
    Are they less stable/precise?

    Thanks again !

    • Nigel Redmon says:

      Hi Malcolm,

      Well, the “j” in exp(-jw) is the square root of minus one, so you have to use trig identities to expand it out then…it’s no longer shorter ;-)

      But I plan to get to some more filters, and look at something more synth-oriented, some time after I get through the last couple of installments on the oscillator…in progress…

  9. Malcolm rest says:

    Dear Nigel,

    I wonder if you can find the time to upload an IIR based formula for “all pass” filter for some “phase correction games”.

    Thanks !

    Malcolm

  10. Kermes says:

    I implemented the code for calculating low pass filter coefficients and got the same coefficients in C code and in GNU Octave (MATLAB-alike software) as in your BiQuad calculator. For retro computing purposes of emulating an ancient sound card, I wish to emulate a fourth order Sallen-Key topology which should implement a Butterworth response to my understanding. Fs=49716, Fc=15392, Q=1.25 for the first BiQuad section, and Fs=49716, Fc=15392, Q=0.5405 for the second.

    But when I tried the coefficients in GNU Octave and tried to plot frequency response of one BiQuad with freqz, it is all wrong. MATLAB (and Octave) seems to want the coefficients in Transposed Direct Form II, while there is no mention what form your calculator gives the coefficients, but I have to assume they are Direct Form I coefficients. In fact based on your Biquad tutorial, the coefficients should be exactly the same in all forms, the arrangement of delays just differ?

    Can you give any pointers how to use those coefficients in MATLAB/Octave to plot a frequency response with freqz command? Do you also know how to plot the frequency response of both BiQuads combined?

    • Nigel Redmon says:

      Yes, the coefficients are the same in all of the direct forms—they are just computational rearrangements.

      The problem you have is that for freqz, the “numerator and denominator coefficients are b and a, respectively” (from the freqz definition). If you look at my derivation here, you’ll see that I have a as the numerator, and b as the denominator (see the fourth equation down). Unfortunately, this is a typical problem—there is no standard “spelling”. I once surveyed all of my DSP books, and the way I spell the coefficients won out—and happened to be the one I prefer—but it just depends on whether you start with the transfer function, where it makes sense to have a on top, or the direct form I implementation, where it makes sense to have a on the left (which is equivalent to a on the bottom of the transfer function).

      • Kermes says:

        Sorry, the problem was that one needs to have also three b coefficients, b1 b2 from your calculator, and b0 is 1. After building a vector b with [1 b1 b2] it works. Of course b and a correctly placed as numerator and denominator. Thank you!

  11. Thanks for this article and the biquad calculator!

    I just discovered this today. I’m trying to add a biquad filter to my own little dsp project. I started out from examples like http://peabody.sapp.org/class/350.838/lab/biquad/ and the pure data source code. I got a basic object ready, based on this algorithm:

    for (Int i = 0; i < length; i++) {
    output = *inPtr++ + fb1 * last + fb2 * previous;
    *outPtr++ = ff1 * output + ff2 * last + ff3 * previous;
    previous = last;
    last = output;
    }

    which works with the example parameters provided. But it would be nice if I could use the object by providing frequency, gain and Q instead of coefs. After looking at this code, that seemed an easy task. Just map a0-2 to ff1-3 and b1-2 to fb1-2, right?

    That didn't work. I haven't written an object yet to plot out what the results are, but at any rate there's no sound and there should be. I've looked at all code I could get my hands on, but I can't figure out what to do different. Do I need to map this differently, or is it not just about mapping?

    I'm sorry if I'm asking a very dumb question here. I'm more an artist and math is not my strongest point.

    Regards,

    yvan

    • Nigel Redmon says:

      Hi Yvan,

      You are correct about the mapping. The code you gave implements the direct form II diagram shown here… [And after further discussion privately, we determined that Yvan needed to negate the feedback coefficients (b1 and b2) as well. This is a common problem—some people use a for the feedforward terms and b for feedback, and others do the opposite; and others subtract the feedback terms and some add negated feedback terms. DSP texts vary on this choice, so you need to pay attention that the coefficient calculations agreed with the filter realization when implementing your biquads.]

      Regards,

      Nigel

  12. Vermeille says:

    Thank you for your amazing blog. It helped me a lot.

    Following your formulas, I don’t see a scale factor, and for Q > 1 values, it may cause overflows. How to calculate a scale factor ?

    Thank you again.

    Vermeille.

    • Nigel Redmon says:

      Hi Vermeille—The calculations are floating point. If you want to implement a filter in fixed-point, you may have to scale and adjust when evaluating the filter. (For instance, coefficients can be in the range of ±2; you’d halve the relevant coefficients, then double products in the filter evaluation.)

      Nigel

  13. Gábor Kakuk says:

    Hy Nigel,

    Could you help me to change the code in HP / LP / BP / etc. to adjust the gain?

    Thank you in advance!

    • Nigel Redmon says:

      That part’s easy—just multiply the result by your gain factor. Or, if you don’t want to add another operation, multiply the coefficients.

      • Gábor Kakuk says:

        Wow!
        Thank youso much!!!

        Another question…or it is an interesting idea I think.

        Is it possible to change the code of the shelf filter (Low & High), to get “asymmetric” slopes?
        Ok, so. The Q could be positive but also negative.
        If it is negative, the only the negative slope will appears on the filters curve, if it is positive, only the positive slope will appears. if it is 0, then no slope.

        • Nigel Redmon says:

          I’m not sure if I follow exactly what you are looking for, but as far as biquads go, there is a limit on possible shapes with only a pair of complex conjugate poles, and zeros. Point me to a response plot if you can.

  14. Niro says:

    Dear Nigel,
    can you tell me, where have you found these formulas (web, book)?
    Thanks

    • Nigel Redmon says:

      All of the formulas were derived from analog (s-domain) prototypes, which describe the frequency responses of the filters, via the bilinear transform, as shown here. Additionally, I got help with the derivation of the shelves from Digital Audio Signal Processing, by Udo Zölzer. Robert Bristow-Johnson’s cookbook (search: rbj cookbook) has some good information on the analog prototypes as well. He uses a equivalent sin/cos substitution instead of the one involving tan that I use, but although the formulas may look different, for the most part they give identical results. However, note that his shelves aren’t symmetrical, so you’ll find the ones here more practical.

Leave a Reply

Your email address will not be published. Required fields are marked *


*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>