1
$\begingroup$

I am fairly new DSP, and am using maths that I have not used in 20-30 years. I have been working with Cookbook formulae for audio EQ biquad filter coefficients by Robert Bristow-Johnson and am having a few problems.

I am trying to design a biquad band pass filter with a sampling rate of 48000 Hz, a center frequency of 20000 Hz, and a bandwidth of 500Hz. Q = fc/delta f, so Q should = 20000/500 = 40. I have run these numbers through a few different tools that use the Robert Bristow-Johnson formulas, and calculated them manually, and get the following coefficients:

a1 = 1.72129273 a2 = 0.98757764 b0 = 0.00621118 b1 = 0.00000000e+0 b2 = -0.00621118

I have plotted the frequency response with various tools (https://arachnoid.com/BiQuadDesigner/, source from http://www.earlevel.com/main/2016/12/01/evaluating-filter-frequency-response/ and others). I expected to see a -3db roll off at about 19750Hz/20250Hz, but I see it at about 19950Hz/20050Hz. Could anyone point out what I am doing wrong.

Thank you

$\endgroup$
  • $\begingroup$ BTW, this answer is about the compensation for the scrunched $BW$ that the bilinear transform requires. $\endgroup$ – robert bristow-johnson Mar 2 '18 at 7:00
1
$\begingroup$

i believe that what you're seeing is the compression of bandwidth done by the frequency warping inherent with the bilinear transform. now, if you look at cookbook, the bandwidth $BW$ is expressed in terms of octaves. so, for your spec:

$$ BW = \log_2(20250) - \log_2(19750) = 0.036069 \text{ octave} $$

try using that number for $BW$ or the corresponding $Q$ as shown in the cookbook:

$$ \frac{1}{Q} = 2 \sinh\left( \frac{\ln(2)}{2} BW \frac{\omega_0}{\sin(\omega_0)} \right) $$

i come up with $Q = 7.63359$.

your formula relating bandwidth in linear frequency to Q is accurate for analog filters but is accurate for digital filters (designed using bilinear transform) only for frequencies much lower than Nyquist. the reason why is that frequency warping of the BLT moves frequency points downward and even after compensating the center frequency, it compresses the bandwidth in comparison to the analog counterpart. and 20 kHz is quite close to Nyquist, so we expect this frequency warping and bandwidth compression.

note that for analog, the equation relating $Q$ and $BW$ (in octaves) is slightly different:

$$ \frac{1}{Q} = 2 \sinh\left( \frac{\ln(2)}{2} BW \right) $$

note that the factor $\frac{\omega_0}{\sin(\omega_0)}$ increases the effect of bandwidth, so the digital bandwidth need be smaller than the analog bandwidth especially as $\omega_0 \to \pi = $ Nyquist.

$\endgroup$
  • $\begingroup$ Your answer seems to suggest that the problem in the question cannot be solved exactly (without trial and error), but that's not the case (and I suppose you know that). There is a simple formula where you put in the center frequency and bandwidth (as a difference of frequencies), and you get a discrete-time biquad satisfying all specs exactly. $\endgroup$ – Matt L. Mar 3 '18 at 10:56
  • $\begingroup$ it's not trial and error. but the $\frac{\omega_0}{\sin(\omega_0)}$ compensation for bandwidth scrunching is not exact, but it's pretty good. and there is no iteration. when you say "satisfying all specs exactly", i am not sure there an unambiguous set of "all specs" exists. particularly considering what the gain at Nyquist is. since this is bilinear transform of a BPF, the gain and Nyquist is $-\infty$ dB, but i wouldn't call that a "spec". $\endgroup$ – robert bristow-johnson Mar 3 '18 at 11:57
  • $\begingroup$ The OP's problem was that the cut-off frequencies were not satisfied exactly. What I mean is that it's possible to exactly implement the two band edges and the center frequency, without the need for any compensation. What you'll get is a filter with zeros at DC and Nyquist, a gain of 1 at the center frequency, and a gain of -3dB at the specified edges. I believe that is what the OP wanted and what he wasn't able to get. $\endgroup$ – Matt L. Mar 3 '18 at 12:53
  • $\begingroup$ actually, the OP's problem was that the cut-off frequencies were way off. they were not what he/she expected using the $Q=\frac{f_0}{\Delta f}$ formula and the OP wanted to know how it could be so far off. and the answer i gave is the reason. the Cookbook compensates (to a first-order approximation) for the scrunching of bandwidth (in octaves) and adjusts the bandwidth accordingly, but did not fudge the definition of Q. the reason why is that for LPF, BPF (constant skirt), HPF there are amplitude features, unchanged by BLT mapping, that are related to Q so i didn't change Q. $\endgroup$ – robert bristow-johnson Mar 3 '18 at 22:41
1
$\begingroup$

As explained in this answer there is no analytical solution for designing a discrete-time biquad band pass filter with its bandwidth given as a ratio (i.e., in octaves). However, there is a straightforward analytical solution to the problem if the bandwidth is given as a difference of frequencies, as asked in the above question. So there is no need to use any approximate formulas, you can simply plug your specs in the formula below and you're done.

Let $\Delta\omega$ be the normalized bandwidth:

$$\Delta\omega=2\pi\frac{\Delta f}{f_s}\tag{1}$$

where $\Delta f$ is the given bandwidth in Hertz, and $f_s$ is the sampling frequency. Furthermore, let $\omega_0$ be the given normalized center frequency:

$$\omega_0=2\pi\frac{f_0}{f_s}\tag{2}$$

with $f_0$ the center frequency in Hertz.

The desired transfer function of the discrete-time biquad is given by

$$H(z)=g\,\frac{z^2-1}{z^2+a_1z+a_2}\tag{3}$$

with

$$g=\frac{\beta}{1+\beta},\quad a_1=-\frac{2\alpha}{1+\beta},\quad a_2=\frac{1-\beta}{1+\beta}\tag{4}$$

where $\alpha$ and $\beta$ are defined by

$$\alpha=\cos(\omega_0),\quad\beta=\tan\left(\frac{\Delta\omega}{2}\right)\tag{5}$$

The transfer function $(3)$ with the constants given by $(4)$ and $(5)$ has the specified center frequency and bandwidth.

Note that the transfer function $(3)$ has zeros at DC ($z=1$) and at Nyquist ($z=-1$). The $3$ variables ($g$, $a_1$, and $a_2$) are used to specify the center frequency $\omega_0$, the bandwidth $\Delta\omega$, and the gain at $\omega_0$ (chosen as $0\text{ dB}$).

With your specs I get the following values:

g =  0.0316988960039693
a1 =  1.67714670914616
a2 =  0.936602207992062

The $-3\text{ dB}$ band edges are at $f_1=19743\text{Hz}$ and at $f_2=20243\text{Hz}$.

$\endgroup$
  • $\begingroup$ this answer gives exact coefficients for the case where the -3 dB bandwidth is specified in Hz, but since the result is equivalent to the design one gets from using the bilinear transform, it also does not satisfy: $$Q=\frac{f_0}{\Delta f}=\frac{\omega_0}{\Delta \omega}$$ if the value of $Q$ remains unchanged between the analog and digital realizations. and it does not satisfy that well-established analog filter formula for the same reason as i stated above. $\endgroup$ – robert bristow-johnson Mar 3 '18 at 22:51
  • $\begingroup$ @robertbristow-johnson: what do you mean by "it does not satisfy that well-established analog filter formula". If you look at the problem in the OP (given center frequency and bandwidth as a difference of frequencies) then the formula given in this answer is exactly what is needed, as far as I can tell. $\endgroup$ – Matt L. Mar 4 '18 at 10:34
  • $\begingroup$ just run the numbers, @MattL. why is the Q ca. 7.6 instead of 40? $\endgroup$ – robert bristow-johnson Mar 5 '18 at 17:24
  • $\begingroup$ @robertbristow-johnson: The center frequency and bandwidth of the discrete-time filter exactly satisfy the specs, that's all we should care about. I really don't see what's wrong with that filter, it does what it should do. $\endgroup$ – Matt L. Mar 5 '18 at 17:30
  • $\begingroup$ Matt, i never said that something is wrong with the filter. in a comment following the other answer, we discussed what it was that the OP was asking and how to answer it. i commented after this answer, that your filter is equivalent to one designed with BLT (because it maps s-plane DC to z-plane DC and s-plane $\infty$ to z-plane $-\pi$ (Nyquist). that mapping must necessarily screw up the use of $Q=\frac{f_0}{\Delta f}$ which is what the OP was using in his/her design. he/she had a reason to expect Q=40 to be the correct choice, yet his/her bandwidth was way too small. why? $\endgroup$ – robert bristow-johnson Mar 5 '18 at 20:22

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.