# Audio EQ Cookbook without frequency warping

The famous http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt offers a set of [real] biquad filter calculation formulas that generally work fine.

However, when filter's frequency approaches Nyquist frequency the Q (bandwidth) specification of a filter becomes greatly distorted - usually it shrinks a lot (even though the author mentioned that he performed a necessary pre-warping).

I'm in quest for filter formulas that do not have such strong bandwidth distortion. I need peaking/bell, band-pass, low-pass, high-pass, high-shelf, notch filters. I know this can be done as previously I've bought a peaking/bell/band-pass filter formulas with less distortion, but they are still not perfect and I need other filter types.

So, I'm also willing to pay for the solution if the price is right.

Alternatively, if one could point me to an optimization algorithm that works with Z-domain filters that would be great, too. Unfortunately, most usual optimization algorithms do not work well in Z-domain - they can't optimize a set of parameters to match a desired frequency response (probably due to periodic functions used to calculate the frequency response).

• In order to effectively compensate the aliasing and warping artefacts for discrete filters you have to use very high orders. FIR would make the design much easier indeed, but there are also other methods. In order to give you a meaningful answer I would like to know what exactly your computational constraints are (including platform), and what your processing requirements are. Why do you need this behaviour, and would oversampling not suffice? – Jazzmaniac Nov 19 '14 at 19:21
• well, there's a limit to what you can get with a 2nd-order IIR filter. there is another paper by Orfanidis that puts the gain at Nyquist to be something other than 0 dB (the Cookbook puts the Nyquist gain at 0 dB because that's how the bilinear transform maps it). for a 44.1 kHz sample rate, and 2nd-order biquad, i don't think you'll do much better than the Orfanidis design. – robert bristow-johnson Nov 19 '14 at 21:07
• The best method I'm aware of is not directly mapping to the discrete domain but rather applying ODE solvers to the analog prototype. Specifically implicit solvers do a great job at preserving the frequency response. – Jazzmaniac Nov 19 '14 at 21:36
• @Jazzmaniac: Why use an analog prototype instead of designing directly in the discrete-time domain? – Matt L. Nov 19 '14 at 22:10
• @aleksv: You could post a concrete design problem with all specs and see what we come up with. – Matt L. Nov 20 '14 at 8:06

here is a quick look at how the 5 degrees of freedom for the parametric EQ can be viewed. it's my take on what Knud Christensen of tc electronic came up with about a decade ago at an AES convention and this patent.

so, forget about the Cookbook (and the issues of Q and bandwidth therein) and consider (in the s-plane) the parametric EQ as the sum of a bandpass filter (with a $Q$ value) in parallel with a wire:

$$H(s) = (G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1$$

$G_\text{boost} = 10^{\frac{dB}{20}}$ is the gain of the peak (or valley, if $dB<0$). the gain at DC and at Nyquist is 0 dB. that's a 2nd-order IIR and there are 3 independent parameters. 2 more to go. so we next add an overall gain parameter:

$$H(s) = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1 \right)$$

that's 4 knobs to twist. one more knob to add (without raising the filter order) and we will be done with adding more independent parameters.

so what Knud does here is replace that "wire" (that trailing "$1$" in the transfer function) with a prototype shelving filter that must have the same poles, the same $Q$ and $\omega_0$ as the BPF, so that the denominator is the same. the transfer function of that shelf is:

$$H_\text{shelf}(s) = \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1}$$

where $R \triangleq 10^{\frac{tilt}{20}}$ and $tilt$ is the gain differential of the shelf in dB. this is what offsets the gain at Nyquist to be different than the gain at DC. After bilinear transformation, Nyquist gets boosted by $tilt$ dB and the gain at DC remains unchanged. Like the $dB$ boost parameter, the $tilt$ parameter can be either positive or negative. $G_0 \cdot R$ is the linear gain at Nyquist.

put that all together and you get:

\begin{align} H(s) & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + H_\text{shelf}(s) \right) \\ & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} \right) \\ & = G_0 \frac{R\left(\frac{s}{\omega_0}\right)^2 + (G_\text{boost}+\sqrt{R}-1)\frac{1}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} \\ \end{align}

no matter how you look at that, this has 5 degrees of freedom and those 5 biquad coefficients are fully defined from these 5 parameters. doesn't matter if you map from $s$ to $z$ using the blinear transform or the trapezoidal rule (effectively the same thing) or any other method that does not change the order of the filter. you may have to fudge the definition of $Q$ or bandwidth, you may have to compensate $\omega_0$ and/or $Q$ for frequency warping effects (like you get with the blinear transform), but if you paid big bucks for something that gets you a 2nd-order IIR filter, it doesn't matter if you implement it with any Direct Form or transposed Direct Form or Lattice or Normalized Ladder or Hal Chamberlin's State-Variable or Andrew Simpson's modeling of linear analog with trapezoidal integration, eventually you get to 5 coefficients and they can be mapped to these 5 independent parameters. it's all the same. whether or not you paid money for a license or not. the math is stronger than any claims made by whomever you are licensing from.

Just FYI, i solved where the true peak or valley frequency is when there is a $tilt$ that is not zero. the frequency where the peak or valley has been nudged over by the tilt is:

$$\omega_\text{peak} \ = \ \omega_0 \ \sqrt{ \frac{Q^2 \left(R - \frac{1}{R}\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \\ \sqrt{\frac{G_\text{boost}^2 - \frac{1}{R} + 2 Q^2 \left(\frac{1}{R} - 1\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \frac{Q^4 \left(R - \frac{1}{R}\right)^2}{\left(G_\text{boost}^2 - R + 2 Q^2 (R - 1)\right)^2} } }$$

you can see that when $tilt = 0$, then $R=1$ and consequently $\omega_\text{peak} = \omega_0$. the peak gain $G_\text{boost}$ might also have to be adjusted a little and that has yet to be worked out. a good first guess would be $G_\text{boost} \leftarrow \frac{G_\text{boost}}{\sqrt{R}}$ or maybe $G_\text{boost} \leftarrow G_\text{boost}-(\sqrt{R}-1)$.

• I cannot quite understand what you are trying to prove? 2nd order IIR filter has 5 degrees of freedom by definition. How your mapping of parameters to biquad coefficients helps to answer the original question? – aleksv Nov 20 '14 at 18:24
• so, once it became clear that "of course it is 2nd order IIR filter I'm talking about", so that increasing the filter order "to effectively compensate the aliasing and warping artefacts for discrete filters [using] very high orders" is not in the cards, then, as best as i can tell from what you said, is that you paid money to license IP from someone to define those 5 coefficients. and all i am trying to say is that this is a fully solved problem. whether you paid for it or not. – robert bristow-johnson Nov 20 '14 at 18:45
• The peaking filter calc code I have also uses 5 parameters: gain at DC, left edge freq and gain, peak freq and gain. I think it's important to specify edge freq to maintain a stable bandwidth. – aleksv Nov 20 '14 at 18:46
• It may be a solved problem, but where I can find a solution? Be it paid or not? – aleksv Nov 20 '14 at 18:47
• okay, now take a look at the Orfanidis paper, but use the math symbols above. this is as good as you can do, because you only have those 5 knobs to twist. you need to figure out what $G_0$ and $A$ are. $\frac{G_0}{A}$ is the DC gain (maybe you want that to be 0 dB or $\frac{G_0}{A} = 1$, i dunno). $A \ G_0$ is the gain at $\omega=\infty$, but, using bilinear transform, that's the gain at Nyquist. so you have to figure out what you want the gain at Nyquist to be. Sophocles has a suggestion. dunno if i agree with it. – robert bristow-johnson Nov 20 '14 at 19:25

@Jazz, one of the things we learned in electrical engineering is that any order of differential equation can be broken down to a set (or "system") of 1st-order diff eqs. so if trapezoidal integration, with the same "time step" $\Delta t$ is being used consistently for all continuous-time integrals, for a $N$th-order linear ODE, you can bust that up into $N$ first-order differential equations. then consider just one of those 1st-order diff eqs:

again, consider emulating a capacitor. let the sampling period be $T=\frac{1}{f_\text{s}}$ be the same as the "$\Delta t$" used in the trapezoid rule.

$$i(t) = C \frac{dv}{dt}$$

or

$$v(t) = \frac{1}{C} \int\limits_{-\infty}^{t} i(u) \ du$$

in the s-domain it's

$$V(s) = \frac{1}{s} \left( \frac{1}{C} I(s) \right)$$

so trapazoidal integration at discrete times is:

\begin{align} v(nT) & = \frac{1}{C} \int\limits_{-\infty}^{nT} i(u) \ du \\ & = \frac{1}{C} \sum\limits_{k=-\infty}^{n} \quad \int\limits_{kT-T}^{kT} i(u) \ du \\ & \approx \frac{1}{C} \sum\limits_{k=-\infty}^{n} \frac{T}{2} \left( i(kT-T) + i(kT) \right) \\ & = v((n-1)T) + \frac{1}{C} \frac{T}{2} \left( i((n-1)T) + i(nT) \right) \\ \end{align}

or as discrete-time sample values

$$v[n] = v[n-1] + \frac{T}{2C} (i[n] + i[n-1])$$

applying the Z transform

$$V(z) = z^{-1} V(z) \ + \ \frac{T}{2C} \left(I(z) + z^{-1} I(z) \right)$$

solving for $V$

$$V(z) = \frac{T}{2} \frac{1 + z^{-1}}{1 - z^{-1}} \left( \frac{1}{C} I(z) \right)$$

looks like we're substituting

$$\frac{1}{s} \leftarrow \frac{T}{2} \frac{1 + z^{-1}}{1 - z^{-1}}$$

or

$$s \leftarrow \frac{2}{T} \frac{z-1}{z+1}$$

which is precisely what the bilinear without compensation for frequency warping does.

• This is what is usually called "proof by example" and hardly general. It also only treats the case that I already said was trivial in the comments above. Again, I'm talking about higher order implicit solvers that don't give a simple discrete recursion. It's also hardly the right place for this discussion as your answer has nothing to do with the original question. If you move this to a new question regarding ODE integration I'm willing to contribute to it. – Jazzmaniac Nov 20 '14 at 18:13
• no, it's completely general for linear diff eqs. you are ignoring the opening statement about we learned long ago in electrical engineering classes. (both in numerical methods and in control theory, an $N$th-order differential equation can be broken down to $N$ first-order differential equations. and then i require consistency; that the same trapezoidal rule for integration is used throughout.) – robert bristow-johnson Nov 20 '14 at 18:29
• I am fully aware of the transformation of higher order ODEs into systems of first order ODEs. And again, it is not what I'm talking about. This statement doesn't make your proof any more general, as you just assert that the bilinear mapping transcends to the entire class while only treating a simple single example. Still, I'm not saying that your result is wrong, I'm talking about different solvers. Higher order solvers are not solvers of higher order DEs. – Jazzmaniac Nov 20 '14 at 18:46
• sorry, Andrew. that's a fail. you truly do not understand if you're continuing to maintain that modeling the continuous-time behavior of the linear analog components using trapezoidal integration, with the step time $\Delta t$ being the same as the sampling period $T=\frac{1}{f_\text{s}}$, then that is the very same thing as using bilinear transform without compensation for frequency warping. same thing and the proof is above. – robert bristow-johnson Nov 20 '14 at 19:18
• Sorry Robert, I'm not Andrew. And your proof above is like the following. I have two maps A and B, and some x. Both A(x)=c and B(x)=c, therefore A = B. That's obviously wrong. Even if A and B are linear, it's still wrong. You're right that a first order trapezoidal solver on a linear system will generate a recursive discrete system, but that's where it ends. You keep ignoring my repeated statements about higher order (adaptive/implicit) integrators. So I see the fail very much on your side. – Jazzmaniac Nov 20 '14 at 19:24

Using optimization methods, we can get a digital filter's frequency response closer to the target analog filter.

In the following experiment, a 6-order bandpass filter is optimized using Adam, an optimization algorithm often used in machine learning. Frequencies above the passband are excluded from the cost function (assigned zero weight). The optimized filter's response becomes higher than target for frequencies very close to Nyquist, but that difference may be offset by the anti-aliasing filter of the signal source (ADC or sample rate converter).

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as clr
from scipy import signal

import tensorflow as tf

# Number of sections
M = 3

# Sample rate
f_s = 24000

# Passband center frequency
f0 = 9000

# Number of frequencies to compute
N = 2048

section_colors = np.zeros([M, 3])
for k in range(M):
section_colors[k] = clr.hsv_to_rgb([(k / (M - 1.0)) / 3.0, 0.5, 0.75])

# Get one of BP poles that maps to LP prototype pole.
def lp_to_bp(s, rbw, w0):
return w0 * (s * rbw / 2 + 1j * np.sqrt(1.0 - np.power(s * rbw / 2, 2)))

# Frequency response

def freq_response(z, b, a):
p = b[0]
q = a[0]
for k in range(1, len(b)):
p += b[k] * np.power(z, -k)
for k in range(1, len(a)):
q += a[k] * np.power(z, -k)
return p / q

# Absolute value in decibel

def abs_db(h):
return 20 * np.log10(np.abs(h))

# Poles of analog low-pass prototype

none, S, none = signal.buttap(M)

# Band limits
c = np.power(2, 1 / 12.0)
f_l = f0 / c
f_u = f0 * c

w0 = 2 * np.pi * f0
w_l = 2 * np.pi * f_l
w_u = 2 * np.pi * f_u

# Relative bandwidth
rbw = (w_u - w_l) / w0

jw0 = 2j * np.pi * f0
z0 = np.exp(jw0 / f_s)

# 1. Analog filter parameters

bc, ac = signal.butter(M, [w_l, w_u], btype='bandpass', analog=True)
ww, H_a = signal.freqs(bc, ac, worN=N)
magnH_a = np.abs(H_a)
f = ww / (2 * np.pi)

omega_d = ww / f_s
z = np.exp(1j * ww / f_s)

# 2. Initial filter design

a = np.zeros([M, 3], dtype=np.double)
b = np.zeros([M, 3], dtype=np.double)
hd = np.zeros([M, N], dtype=np.complex)

# Pre-warp the frequencies

w_l_pw = 2 * f_s * np.tan(np.pi * f_l / f_s)
w_u_pw = 2 * f_s * np.tan(np.pi * f_u / f_s)
w_0_pw = np.sqrt(w_l_pw * w_u_pw)

rbw_pw = (w_u_pw - w_l_pw) / w_0_pw

poles_pw = lp_to_bp(S, rbw_pw, w_0_pw)

# Bilinear transform

T = 1.0 / f_s
poles_d = (1.0 + poles_pw * T / 2) / (1.0 - poles_pw * T / 2)

for k in range(M):
p = poles_d[k]
b[k], a[k] = signal.zpk2tf([-1, 1], [p, np.conj(p)], 1)

g0 = freq_response(z0, b[k], a[k])
g0 = np.abs(g0)
b[k] /= g0
none, hd[k] = signal.freqz(b[k], a[k], worN=omega_d)

plt.figure(2)
plt.title("Initial digital filter (bilinear)")

plt.axis([0, f_s / 2, -90, 10])

plt.plot(f, abs_db(H_a), label='Target response', color='gray', linewidth=0.5)

for k in range(M):
label = "Section %d" % k
plt.plot(f, abs_db(hd[k]), color=section_colors[k], alpha=0.5, label=label)

# Combined frequency response of initial digital filter

Hd = np.prod(hd, axis=0)
plt.legend(loc='upper left')

plt.figure(3)
plt.title("Initial filter - poles and zeros")
plt.axis([-3, 3, -2.25, 2.25])
unitcircle = plt.Circle((0, 0), 1, color='lightgray', fill=False)
ax = plt.gca()

for k in range(M):
zeros, poles, gain = signal.tf2zpk(b[k], a[k])
plt.plot(np.real(poles), np.imag(poles), 'x', color=section_colors[k])
plt.plot(np.real(zeros), np.imag(zeros), 'o', color='none', markeredgecolor=section_colors[k], alpha=0.5)

# Optimizing filter

tH_a = tf.constant(magnH_a, dtype=tf.float32)

# Assign weights

weight = np.zeros(N)
for i in range(N):
# In the passband or below?
if (f[i] <= f_u):
weight[i] = 1.0

tWeight = tf.constant(weight, dtype=tf.float32)
tZ = tf.placeholder(tf.complex64, [1, N])

# Variables to be changed by optimizer
ta = tf.Variable(a)
tb = tf.Variable(b)
ai = a
bi = b

# TF requires matching types for multiplication;
# cast real coefficients to complex
cta = tf.cast(ta, tf.complex64)
ctb = tf.cast(tb, tf.complex64)

xb0 = tf.reshape(ctb[:, 0], [M, 1])
xb1 = tf.reshape(ctb[:, 1], [M, 1])
xb2 = tf.reshape(ctb[:, 2], [M, 1])

xa0 = tf.reshape(cta[:, 0], [M, 1])
xa1 = tf.reshape(cta[:, 1], [M, 1])
xa2 = tf.reshape(cta[:, 2], [M, 1])

# Numerator:   B = b₀z² + b₁z + b₂
tB = tf.matmul(xb0, tf.square(tZ)) + tf.matmul(xb1, tZ) + xb2

# Denominator: A = a₀z² + a₁z + a₂
tA = tf.matmul(xa0, tf.square(tZ)) + tf.matmul(xa1, tZ) + xa2

# Get combined frequency response
tH = tf.reduce_prod(tB / tA, axis=0)

iterations = 2000
learning_rate = 0.0005

# Cost function
cost = tf.reduce_mean(tWeight * tf.squared_difference(tf.abs(tH), tH_a))

zz = np.reshape(z, [1, N])

with tf.Session() as sess:
sess.run(tf.global_variables_initializer())

for epoch in range(iterations):
loss, j = sess.run([optimizer, cost], feed_dict={tZ: zz})
if (epoch % 100 == 0):
print("  Cost: ", j)

b, a = sess.run([tb, ta])

for k in range(M):
none, hd[k] = signal.freqz(b[k], a[k], worN=omega_d)

plt.figure(4)
plt.title("Optimized digital filter")

plt.axis([0, f_s / 2, -90, 10])

# Draw the band limits
plt.axvline(f_l, color='black', linewidth=0.5, linestyle='--')
plt.axvline(f_u, color='black', linewidth=0.5, linestyle='--')

plt.plot(f, abs_db(H_a), label='Target response', color='gray', linewidth=0.5)

Hd = np.prod(hd, axis=0)
for k in range(M):
label = "Section %d" % k
plt.plot(f, abs_db(hd[k]), color=section_colors[k], alpha=0.5, label=label)

magnH_d = np.abs(Hd)
plt.legend(loc='upper left')

plt.figure(5)
plt.title("Optimized digital filter - Poles and Zeros")
plt.axis([-3, 3, -2.25, 2.25])
unitcircle = plt.Circle((0, 0), 1, color='lightgray', fill=False)
ax = plt.gca()

for k in range(M):
zeros, poles, gain = signal.tf2zpk(b[k], a[k])
plt.plot(np.real(poles), np.imag(poles), 'x', color=section_colors[k])
plt.plot(np.real(zeros), np.imag(zeros), 'o', color='none', markeredgecolor=section_colors[k], alpha=0.5)

plt.show()

• "Frequencies above the passband are excluded from the cost function (assigned zero weight)" Isn't there a danger that the filter becomes a sort of high-pass filter? Frequencies above passband could have large gain with no penalty given by such a cost function. – Olli Niemitalo Jan 21 '19 at 8:48
• OK, II thought a symbol regression method would produce a finalized equation...probably for most cases,not all. while real-time derivative-free? optimization is a bit slow – aleksv Jan 22 '19 at 4:15
• @OlliNiemitalo: For some frequencies, it drops to zero after passband, like the initial filter. – igorinov Jan 22 '19 at 6:36
• @aleksv: There is another closed form solution that has not been mentioned here - khabdha.org/wp-content/uploads/2008/03/… – igorinov Jan 22 '19 at 6:42
• With tensorflow-gpu 1.12.0 the optimization is not converging, rather diverging. tensorflow 1.12.0 works better but the optimization does not fully converge, because of the choice of the optimizer (Adam) and its learning rate, and possibly because of the way the cost function is formulated. I will award the bounty to this answer, because it shows a modern approach to optimization of coefficients: automatic differentiation. I would still like to see: more stable optimization and arranging obtained coefficients for real-time use (such as parameter sweeps). – Olli Niemitalo Jan 22 '19 at 8:22

I've come up with a design for the 10dB peaking EQ. I've chosen 20 filters with center frequencies between 500 Hz and 16 kHz (Fs = 48 kHz). The top plot below is the design according to RBJ's Audio-EQ-Cookbook, which is good but which leads to bandwidth distortion when the center frequencies get closer to Nyquist. The bottom plot is the new design where the filters very closely match the analog prototype filters:

And this is how the new notch filters look like compared to the Cookbook (bandwidth = 4 octaves, highest $f_0=23$ kHz):

The following figure shows a lowpass filter design ($Q=2$, $f_0=16$ kHz, $F_s=48$ kHz). Note that the new design approximates the analog prototype and for this reason it does not perform as a conventional lowpass filter (it has no zero at Nyquist):

• Matt, i might suggest plotting the "peakingEQ" for, say, a 10 dB boost with a fixed $BW$, say 1 octave, and for a variety of peak frequencies, going up to close to Nyquist. even though the Cookbook does a first-order compensation on $BW$ by multiplying by $\frac{\omega_0 T}{\sin(\omega_0 T)}$, it still cannot change the fact that the frequency response must be continuous and have continuous derivatives and it must be symmetrical about Nyquist. even if you bump up the gain at Nyquist as does Orfanidis. that, in and of itself, is the source of this warped frequency response. – robert bristow-johnson Nov 21 '14 at 14:20
• i think he's expecting to see a digital filter, with a sample rate as low as 44.1 kHz, match the analog filter, even if the resonant frequency is just below Nyquist. this is something i had been worrying about 2 decades ago and all i could think of doing easily is simply to compensate for the scrunching of bandwidth that bilinear transform does. then Orfanidis changed one assumption (that the gain at Nyquist is 0 dB), and that helped a little more. but the thing that really helps, is simply to bump up the sample rate. – robert bristow-johnson Nov 21 '14 at 15:28
• Thanks for an attempt. The first Cookbook plot is what I want, but as Robert also noted, with the gain at Nyquist being non-zero. It will be also beneficial to draw the responses on the logarithmic scale in the range e.g. 20 Hz to Nyquist. Audio EQs rarely use a linear frequency scale. And bandwidth is usually specified in octaves. Any biquads have no problems with center/corner frequency up to about 300 Hz, but going higher reveals the issue. – aleksv Nov 21 '14 at 16:23
• No, I mean 0 to 300 Hz, lower frequencies. Above 300 Hz up to Nyquist these Cookbook filter formulas produce progressively increasing filter shape distortion. Of course, the sample rate is 44100 or 48000 as usual. – aleksv Nov 22 '14 at 10:57
• You have to plot the graph on the log scale. It is bandwidth distortion on the LOG scale what bothers me. There's no distortion present in the lower frequencies when seen on the log scale. Peaking, low-pass, high-shelf, notch filters always have 0 dB (unity) gain at DC. – aleksv Nov 22 '14 at 12:36