Parametric Equalization

Julius O. Smith III
CCRMA, Stanford University
	Nov. 6, 1992

This Mathematica notebook is available online at CCRMA  via ftp to ccrma-ftp.stanford.edu, directory pub/DSP/.

Problem Statement

Design a second-order parametric equalization filter section to be used in cascade with similar sections for the purpose of introducing a narrowband boost or cut in the frequency response of the equalizer.

Theory of Design

A second-order filter can have at most two poles and two zeros.  Consider the case of a boost.  The story for a cut is then obtained by interchanging the role of poles and zeros.  

The pole angles are determined by the center frequency f of the boost.  The angle of one pole must be close to

	thp = 2 Pi f / fs, 
	
where fs is the sampling rate, and the other pole angle must be at -thp in order to obtain real filter coefficients.   The pole radii must be equal to the same value rp which is determined by the bandwidth of the boost.  An approximate formula which is most accurate for narrow bandwidths is given by

	rp = Exp[- Pi b / fs],
	
where b is the desired bandwidth.

Because we want the effect of the boost to be localized to a bandwidth b about the center frequency f, it is common practice to place two zeros “behind” the poles, i.e., along the same radial lines in the complex plane but at slightly smaller radii.  This has the effect of “canceling” the poles far away from the center frequency.  (This can be quickly seen by considering the graphical method for computing the frequency response from a pole-zero diagram.)  Assuming a narrowband boost, the effect of a pole-zero pair at angle thp can be considered in isolation, ignoring the effect of the pole-zero pair at angle -thp which approximately cancel at angle thp).  This transfer function of the isolated pole-zero pair is given by

	H[z] = (1 - rz Exp[I thp] / z) / (1 - rp Exp[I thp] / z)
	
where rz is the pole radius, and z is the complex-plane variable.  At frequency f, we have the frequency response

	H[Exp[I thp]]  = ( 1 - rz ) / (1 - rp )

Given a desired gain g at the center frequency, we can solve for rz in terms of rp to obtain

	rz = 1 - g (1-rp)

Since we have assumed a narrowband boost, 1-rp will always be small.  The cancelation of the pole-zero pair far away from the center frequency requires g(1-rp) be small also, i.e., rz should remain close to 1 while being smaller than rp.

Derivation

fs = 44100;	(* sampling rate *)
f  = 5000;	(* frequency of boost or cut *)
b  = 100;	(* bandwidth of boost or cut *)
gd = 10;	(* gain of boost or cut in dB *)

thp = N[2 Pi f / fs]		(* pole angle in z plane *)

rp = N[Exp[- Pi b / fs]]	(* pole radius *)

g = N[10^(gd/20)]

rz = 1 - g (1-rp)	(* zero radius for desired gain g *)

Amplitude response of two zeros:

fr2[r_,th_,ph_] := N[ 1 - 2 r Cos[th] Exp[- I ph] 
			+ r^2 Exp[-2I ph] ];

Utility for converting complex amplitude to decibels

db[c_] := N[ 20 Log[10,Abs[c]] ];

Plot[db[fr2[rz, thp, ph]/fr2[rp,thp,ph]], {ph,0,Pi},
	PlotRange->All];

Pole-Zero Diagram

Create a plot of the unit circle in the complex plane.

uc = ParametricPlot[{Sin[t],Cos[t]},{t,-Pi,Pi},
	AspectRatio->Automatic];

Display poles as "X"s inside the unit circle in the complex plane:

xp = rp Cos[thp];
yp = rp Sin[thp];
xz = rz Cos[thp];
yz = rz Sin[thp];
Show[uc,Graphics[
	Text[X,{xp,yp}],Text[X,{xp,-yp}],
	Text[O,{xz,yz}],Text[O,{xz,-yz}]]];

Interactive "Spreadsheet" Cell for Trying Cases

f  = 8000;			(* eq frequency *)
b  = 100;			(* filter bandwidth *)
gd = 6;

thp = N[2 Pi f / fs];		(* pole angle *)
g = N[10^(gd/20)];		(* gain at resonance *)
rp = N[Exp[- Pi b / fs]];	(* pole radius *)
rz = 1 - g (1-rp);		(* zero radius *)
Plot[db[fr2[rz, thp, ph]/fr2[rp,thp,ph]], {ph,0,Pi},
	PlotRange->All];

Multiple Parametric Eq Sections

f1  = 4000;
b1  = 100;
gd1 = 3;

f2  = 8000;
b2  = 200;
gd2 = -3;

f3  = 12000;
b3  = 300;
gd3 = 3;

thp1 = N[2 Pi f1 / fs];		(* pole angle *)
g1 = N[10^(gd1/20)];		(* gain at resonance *)
rp1 = N[Exp[- Pi b1 / fs]];	(* pole radius *)
rz1 = 1 - g1 (1-rp1);		(* zero radius *)

thp2 = N[2 Pi f2 / fs];		(* pole angle *)
g2 = N[10^(gd2/20)];		(* gain at resonance *)
rp2 = N[Exp[- Pi b2 / fs]];	(* pole radius *)
rz2 = 1 - g2 (1-rp2);		(* zero radius *)

thp3 = N[2 Pi f3 / fs];		(* pole angle *)
g3 = N[10^(gd3/20)];		(* gain at resonance *)
rp3 = N[Exp[- Pi b3 / fs]];	(* pole radius *)
rz3 = 1 - g3 (1-rp3);		(* zero radius *)

Plot[db[
	(fr2[rz1, thp1, ph]/fr2[rp1,thp1,ph])
	(fr2[rz2, thp2, ph]/fr2[rp2,thp2,ph])
	(fr2[rz3, thp3, ph]/fr2[rp3,thp3,ph])	
	], {ph,0,Pi},
	PlotRange->All];