|
| 1 | + |
| 2 | + Cookbook formulae for audio EQ biquad filter coefficients |
| 3 | +---------------------------------------------------------------------------- |
| 4 | + by Robert Bristow-Johnson < [email protected]> |
| 5 | + |
| 6 | + |
| 7 | +All filter transfer functions were derived from analog prototypes (that |
| 8 | +are shown below for each EQ filter type) and had been digitized using the |
| 9 | +Bilinear Transform. BLT frequency warping has been taken into account for |
| 10 | +both significant frequency relocation (this is the normal "prewarping" that |
| 11 | +is necessary when using the BLT) and for bandwidth readjustment (since the |
| 12 | +bandwidth is compressed when mapped from analog to digital using the BLT). |
| 13 | + |
| 14 | +First, given a biquad transfer function defined as: |
| 15 | + |
| 16 | + b0 + b1*z^-1 + b2*z^-2 |
| 17 | + H(z) = ------------------------ (Eq 1) |
| 18 | + a0 + a1*z^-1 + a2*z^-2 |
| 19 | + |
| 20 | +This shows 6 coefficients instead of 5 so, depending on your architechture, |
| 21 | +you will likely normalize a0 to be 1 and perhaps also b0 to 1 (and collect |
| 22 | +that into an overall gain coefficient). Then your transfer function would |
| 23 | +look like: |
| 24 | + |
| 25 | + (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2 |
| 26 | + H(z) = --------------------------------------- (Eq 2) |
| 27 | + 1 + (a1/a0)*z^-1 + (a2/a0)*z^-2 |
| 28 | + |
| 29 | +or |
| 30 | + |
| 31 | + 1 + (b1/b0)*z^-1 + (b2/b0)*z^-2 |
| 32 | + H(z) = (b0/a0) * --------------------------------- (Eq 3) |
| 33 | + 1 + (a1/a0)*z^-1 + (a2/a0)*z^-2 |
| 34 | + |
| 35 | + |
| 36 | +The most straight forward implementation would be the "Direct Form 1" |
| 37 | +(Eq 2): |
| 38 | + |
| 39 | + y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] |
| 40 | + - (a1/a0)*y[n-1] - (a2/a0)*y[n-2] (Eq 4) |
| 41 | + |
| 42 | +This is probably both the best and the easiest method to implement in the |
| 43 | +56K and other fixed-point or floating-point architechtures with a double |
| 44 | +wide accumulator. |
| 45 | + |
| 46 | + |
| 47 | + |
| 48 | +Begin with these user defined parameters: |
| 49 | + |
| 50 | + Fs (the sampling frequency) |
| 51 | + |
| 52 | + f0 ("wherever it's happenin', man." Center Frequency or |
| 53 | + Corner Frequency, or shelf midpoint frequency, depending |
| 54 | + on which filter type. The "significant frequency".) |
| 55 | + |
| 56 | + dBgain (used only for peaking and shelving filters) |
| 57 | + |
| 58 | + Q (the EE kind of definition, except for peakingEQ in which A*Q is |
| 59 | + the classic EE Q. That adjustment in definition was made so that |
| 60 | + a boost of N dB followed by a cut of N dB for identical Q and |
| 61 | + f0/Fs results in a precisely flat unity gain filter or "wire".) |
| 62 | + |
| 63 | + _or_ BW, the bandwidth in octaves (between -3 dB frequencies for BPF |
| 64 | + and notch or between midpoint (dBgain/2) gain frequencies for |
| 65 | + peaking EQ) |
| 66 | + |
| 67 | + _or_ S, a "shelf slope" parameter (for shelving EQ only). When S = 1, |
| 68 | + the shelf slope is as steep as it can be and remain monotonically |
| 69 | + increasing or decreasing gain with frequency. The shelf slope, in |
| 70 | + dB/octave, remains proportional to S for all other values for a |
| 71 | + fixed f0/Fs and dBgain. |
| 72 | + |
| 73 | + |
| 74 | + |
| 75 | +Then compute a few intermediate variables: |
| 76 | + |
| 77 | + A = sqrt( 10^(dBgain/20) ) |
| 78 | + = 10^(dBgain/40) (for peaking and shelving EQ filters only) |
| 79 | + |
| 80 | + w0 = 2*pi*f0/Fs |
| 81 | + |
| 82 | + cos(w0) |
| 83 | + sin(w0) |
| 84 | + |
| 85 | + alpha = sin(w0)/(2*Q) (case: Q) |
| 86 | + = sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) ) (case: BW) |
| 87 | + = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 ) (case: S) |
| 88 | + |
| 89 | + FYI: The relationship between bandwidth and Q is |
| 90 | + 1/Q = 2*sinh(ln(2)/2*BW*w0/sin(w0)) (digital filter w BLT) |
| 91 | + or 1/Q = 2*sinh(ln(2)/2*BW) (analog filter prototype) |
| 92 | + |
| 93 | + The relationship between shelf slope and Q is |
| 94 | + 1/Q = sqrt((A + 1/A)*(1/S - 1) + 2) |
| 95 | + |
| 96 | + 2*sqrt(A)*alpha = sin(w0) * sqrt( (A^2 + 1)*(1/S - 1) + 2*A ) |
| 97 | + is a handy intermediate variable for shelving EQ filters. |
| 98 | + |
| 99 | + |
| 100 | +Finally, compute the coefficients for whichever filter type you want: |
| 101 | + (The analog prototypes, H(s), are shown for each filter |
| 102 | + type for normalized frequency.) |
| 103 | + |
| 104 | + |
| 105 | +LPF: H(s) = 1 / (s^2 + s/Q + 1) |
| 106 | + |
| 107 | + b0 = (1 - cos(w0))/2 |
| 108 | + b1 = 1 - cos(w0) |
| 109 | + b2 = (1 - cos(w0))/2 |
| 110 | + a0 = 1 + alpha |
| 111 | + a1 = -2*cos(w0) |
| 112 | + a2 = 1 - alpha |
| 113 | + |
| 114 | + |
| 115 | + |
| 116 | +HPF: H(s) = s^2 / (s^2 + s/Q + 1) |
| 117 | + |
| 118 | + b0 = (1 + cos(w0))/2 |
| 119 | + b1 = -(1 + cos(w0)) |
| 120 | + b2 = (1 + cos(w0))/2 |
| 121 | + a0 = 1 + alpha |
| 122 | + a1 = -2*cos(w0) |
| 123 | + a2 = 1 - alpha |
| 124 | + |
| 125 | + |
| 126 | + |
| 127 | +BPF: H(s) = s / (s^2 + s/Q + 1) (constant skirt gain, peak gain = Q) |
| 128 | + |
| 129 | + b0 = sin(w0)/2 = Q*alpha |
| 130 | + b1 = 0 |
| 131 | + b2 = -sin(w0)/2 = -Q*alpha |
| 132 | + a0 = 1 + alpha |
| 133 | + a1 = -2*cos(w0) |
| 134 | + a2 = 1 - alpha |
| 135 | + |
| 136 | + |
| 137 | +BPF: H(s) = (s/Q) / (s^2 + s/Q + 1) (constant 0 dB peak gain) |
| 138 | + |
| 139 | + b0 = alpha |
| 140 | + b1 = 0 |
| 141 | + b2 = -alpha |
| 142 | + a0 = 1 + alpha |
| 143 | + a1 = -2*cos(w0) |
| 144 | + a2 = 1 - alpha |
| 145 | + |
| 146 | + |
| 147 | + |
| 148 | +notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1) |
| 149 | + |
| 150 | + b0 = 1 |
| 151 | + b1 = -2*cos(w0) |
| 152 | + b2 = 1 |
| 153 | + a0 = 1 + alpha |
| 154 | + a1 = -2*cos(w0) |
| 155 | + a2 = 1 - alpha |
| 156 | + |
| 157 | + |
| 158 | + |
| 159 | +APF: H(s) = (s^2 - s/Q + 1) / (s^2 + s/Q + 1) |
| 160 | + |
| 161 | + b0 = 1 - alpha |
| 162 | + b1 = -2*cos(w0) |
| 163 | + b2 = 1 + alpha |
| 164 | + a0 = 1 + alpha |
| 165 | + a1 = -2*cos(w0) |
| 166 | + a2 = 1 - alpha |
| 167 | + |
| 168 | + |
| 169 | + |
| 170 | +peakingEQ: H(s) = (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1) |
| 171 | + |
| 172 | + b0 = 1 + alpha*A |
| 173 | + b1 = -2*cos(w0) |
| 174 | + b2 = 1 - alpha*A |
| 175 | + a0 = 1 + alpha/A |
| 176 | + a1 = -2*cos(w0) |
| 177 | + a2 = 1 - alpha/A |
| 178 | + |
| 179 | + |
| 180 | + |
| 181 | +lowShelf: H(s) = A * (s^2 + (sqrt(A)/Q)*s + A)/(A*s^2 + (sqrt(A)/Q)*s + 1) |
| 182 | + |
| 183 | + b0 = A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha ) |
| 184 | + b1 = 2*A*( (A-1) - (A+1)*cos(w0) ) |
| 185 | + b2 = A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha ) |
| 186 | + a0 = (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha |
| 187 | + a1 = -2*( (A-1) + (A+1)*cos(w0) ) |
| 188 | + a2 = (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha |
| 189 | + |
| 190 | + |
| 191 | + |
| 192 | +highShelf: H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A) |
| 193 | + |
| 194 | + b0 = A*( (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha ) |
| 195 | + b1 = -2*A*( (A-1) + (A+1)*cos(w0) ) |
| 196 | + b2 = A*( (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha ) |
| 197 | + a0 = (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha |
| 198 | + a1 = 2*( (A-1) - (A+1)*cos(w0) ) |
| 199 | + a2 = (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha |
| 200 | + |
| 201 | + |
| 202 | + |
| 203 | + |
| 204 | + |
| 205 | +FYI: The bilinear transform (with compensation for frequency warping) |
| 206 | +substitutes: |
| 207 | + |
| 208 | + 1 1 - z^-1 |
| 209 | + (normalized) s <-- ----------- * ---------- |
| 210 | + tan(w0/2) 1 + z^-1 |
| 211 | + |
| 212 | + and makes use of these trig identities: |
| 213 | + |
| 214 | + sin(w0) 1 - cos(w0) |
| 215 | + tan(w0/2) = ------------- (tan(w0/2))^2 = ------------- |
| 216 | + 1 + cos(w0) 1 + cos(w0) |
| 217 | + |
| 218 | + |
| 219 | + resulting in these substitutions: |
| 220 | + |
| 221 | + |
| 222 | + 1 + cos(w0) 1 + 2*z^-1 + z^-2 |
| 223 | + 1 <-- ------------- * ------------------- |
| 224 | + 1 + cos(w0) 1 + 2*z^-1 + z^-2 |
| 225 | + |
| 226 | + |
| 227 | + 1 + cos(w0) 1 - z^-1 |
| 228 | + s <-- ------------- * ---------- |
| 229 | + sin(w0) 1 + z^-1 |
| 230 | + |
| 231 | + 1 + cos(w0) 1 - z^-2 |
| 232 | + = ------------- * ------------------- |
| 233 | + sin(w0) 1 + 2*z^-1 + z^-2 |
| 234 | + |
| 235 | + |
| 236 | + 1 + cos(w0) 1 - 2*z^-1 + z^-2 |
| 237 | + s^2 <-- ------------- * ------------------- |
| 238 | + 1 - cos(w0) 1 + 2*z^-1 + z^-2 |
| 239 | + |
| 240 | + |
| 241 | + The factor: |
| 242 | + |
| 243 | + 1 + cos(w0) |
| 244 | + ------------------- |
| 245 | + 1 + 2*z^-1 + z^-2 |
| 246 | + |
| 247 | + is common to all terms in both numerator and denominator, can be factored |
| 248 | + out, and thus be left out in the substitutions above resulting in: |
| 249 | + |
| 250 | + |
| 251 | + 1 + 2*z^-1 + z^-2 |
| 252 | + 1 <-- ------------------- |
| 253 | + 1 + cos(w0) |
| 254 | + |
| 255 | + |
| 256 | + 1 - z^-2 |
| 257 | + s <-- ------------------- |
| 258 | + sin(w0) |
| 259 | + |
| 260 | + |
| 261 | + 1 - 2*z^-1 + z^-2 |
| 262 | + s^2 <-- ------------------- |
| 263 | + 1 - cos(w0) |
| 264 | + |
| 265 | + |
| 266 | + In addition, all terms, numerator and denominator, can be multiplied by a |
| 267 | + common (sin(w0))^2 factor, finally resulting in these substitutions: |
| 268 | + |
| 269 | + |
| 270 | + 1 <-- (1 + 2*z^-1 + z^-2) * (1 - cos(w0)) |
| 271 | + |
| 272 | + s <-- (1 - z^-2) * sin(w0) |
| 273 | + |
| 274 | + s^2 <-- (1 - 2*z^-1 + z^-2) * (1 + cos(w0)) |
| 275 | + |
| 276 | + 1 + s^2 <-- 2 * (1 - 2*cos(w0)*z^-1 + z^-2) |
| 277 | + |
| 278 | + |
| 279 | + The biquad coefficient formulae above come out after a little |
| 280 | + simplification. |
0 commit comments