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