Free Books
   Free Books    Introduction to Digital Filters  

Elementary Audio Digital Filters

This appendix is devoted to small but useful digital filters that are commonly used in audio applications. Analytical tools from the main chapters are used to analyze these ``audio gems''.

Elementary Filter Sections

This section gives condensed analysis summaries of the four most elementary digital filters: the one-zero, one-pole, two-pole, and two-zero filters. Despite their relative simplicity, they are quite valuable to master in practice. In particular, recall from Chapter 9 that every causal, finite-order, LTI filter (any difference equation of the form Eq.$ \,$(5.1)) may be factored into a series and/or parallel combinationof such sections. Implementing high-order filters as parallel and/or series combinations of low-order sections offers several advantages, such as numerical robustness and easier/safer control in real time.

One-Zero

Figure B.1: Signal flow graph for the general one-zero filter
$ y(n) = b_0x(n) + b_1 x(n - 1)$.
\begin{figure}\input fig/kfig2p17.pstex_t \end{figure}

Figure B.1 gives the signal flow graph for the general one-zero filter. The frequency response for the one-zero filter may be found by the following steps:

\fbox{ \begin{tabular}{rl} Difference equation: & $y(n) = b_0x(n) + b_1x(n - 1)$... ...requency response: & $H(e^{j\omega T}) = b_0 + b_1e^{-j\omega T}$ \end{tabular}}

By factoring out $ e^{-j\omega T/2}$ from the frequency response, to balance the exponents of $ e$, we can get this closer to polar form as follows:

\begin{eqnarray*} H(e^{j\omega T}) &=& b_0 + b_1 e^{-j\omega T}\\ &=& (b_0 - ... ...omega T/2)\\ &=& (b_0 - b_1) + e^{-j\pi f T} 2b_1\cos(\pi f T) \end{eqnarray*}

Figure B.2: Family of frequency responses of the one-zero filter
$ y(n) = x(n) + b_1 x(n - 1)$
for various values of $ b_1$. (a) Amplitude response. (b) Phase response.
\begin{figure}\input fig/kfig2p19.pstex_t \end{figure}

We now apply the general equations given in Chapter 7 for filter gain $ G(\omega)$ and filter phase $ \Theta(\omega)$ as a function of frequency:

\begin{eqnarray*} H(e^{j\omega T}) &=& b_0 + b_1e^{-j\omega T}\\ &=& b_0 + b_1... ...left[\frac{-b_1 \sin(\omega T)}{b_0 + b_1 \cos(\omega T)}\right] \end{eqnarray*}

A plot of $ G(\omega)$ and $ \Theta(\omega)$ for $ b_0 = 1$ and various real values of $ b_1$, is given in Fig.B.2. The filter has a zero at $ z = -b_1/b_0 = -b_1$ in the $ z$ plane, which is always on the real axis. When a point on the unit circle comes close to the zero of the transfer function the filter gain at that frequency is low. Notice that one real zero can basically make either a highpass ( $ b_1/b_0 < 0$) or a lowpass filter ( $ b_1/b_0 > 0$). For the phase response calculation using the graphical method, it is necessary to include the pole at $ z=0$.

One-Pole

Figure B.3: Signal flow graph for the general one-pole filter
$ y(n) = b_0 x(n) - a_1 y(n - 1).$
\begin{figure}\input fig/kfig2p20.pstex_t \end{figure}

Fig.B.3 gives the signal flow graph for the general one-pole filter. The road to the frequency response goes as follows:

Figure B.4: Family of frequency responses of the one-pole filter
$ y(n) = x(n) - a_1 y(n - 1)$
for various real values of $ a_1$. (a) Amplitude response. (b) Phase response.
\begin{figure}\input fig/kfig2p21.pstex_t \end{figure}

\fbox{ \begin{tabular}{rl} Difference equation: & $y(n) = b_0 x(n) - a_1 y(n-1)... ...$H(e^{j\omega T}) = \displaystyle\frac{b_0}{1+a_1e^{-j\omega T}}$ \end{tabular}}

The one-pole filter has a transfer function (hence frequency response) which is the reciprocal of that of a one-zero. The analysis is thus quite analogous. The frequency response in polar form is given by

\begin{eqnarray*} G(\omega) &=& \frac{\vert b_0\vert}{\sqrt{[1 + a_1 \cos(\omega... ... + a_1 \cos(\omega T)}\right], & b_0<0 \\ \end{array} \right.. \end{eqnarray*}

A plot of the frequency response in polar form for $ b_0 = 1$ and various values of $ a_1$ is given in Fig.B.4.

The filter has a pole at $ z = -a_1$, in the $ z$ plane (and a zero at $ z$ = 0). Notice that the one-pole exhibits either a lowpass or a highpass frequency response, like the one-zero. The lowpass character occurs when the pole is near the point $ z = 1$ (dc), which happens when $ a_1$ approaches $ -1$. Conversely, the highpass nature occurs when $ a_1$ is positive.

The one-pole filter section can achieve much more drastic differences between the gain at high frequencies and the gain at low frequencies than can the one-zero filter. This difference is achieved in the one-pole by gain boost in the passband rather than attenuation in the stopband; thus it is usually desirable when using a one-pole filter to set $ b_0$ to a small value, such as $ 1 - \left\vert a_1\right\vert$, so that the peak gain is 1 or so. When the peak gain is 1, the filter is unlikely to overflow.B.1

Finally, note that the one-pole filter is stable if and only if $ \left\vert a_1\right\vert < 1$.

Two-Pole

Figure B.5: Signal flow graph for the general two-pole filter
$ y(n) = b_0 x(n) - a_1 y(n - 1) - a_2 y(n - 2)$.
\begin{figure}\input fig/twopole.pstex_t \end{figure}

The signal flow graph for the general two-pole filter is given in Fig.B.5. We proceed as usual with the general analysis steps to obtain the following:

\fbox{ \begin{tabular}{rl} Difference equation: & $y(n) = b_0 x(n) - a_1 y(n-1)... ... \displaystyle\frac{b_0}{1+a_1e^{-j\omega T}+a_2e^{-j2\omega T}}$ \end{tabular}}

The numerator of $ H(z)$ is a constant, so there are no zeros other than two at the origin of the $ z$ plane.

The coefficients $ a_1$ and $ a_2$ are called the denominator coefficients, and they determine the two poles of $ H(z)$. Using the quadratic formula, the poles are found to be located at

$\displaystyle z = -\frac{a_1}{2} \pm \sqrt{\left(\frac{a_1}{2}\right)^2 -a_2}. $

When the coefficients $ a_1$ and $ a_2$ are real (as we typically assume), the poles must be either real (when $ (a_1/2)^2\geq a_2$) or form a complex conjugate pair (when $ (a_1/2)^2<a_2$).

When both poles are real, the two-pole can be analyzed simply as a cascade of two one-pole sections, as in the previous section. That is, one can multiply pointwise two magnitude plots such as Fig.B.4a, and add pointwise two phase plots such as Fig.B.4b.

When the poles are complex, they can be written as

\begin{eqnarray*} p_1 &=& x_p + j y_p\\ p_2 &=& x_p - j y_p = \overline{p}_1 \end{eqnarray*}

since they must form a complex-conjugate pair when $ a_1$ and $ a_2$ are real. We may express them in polar form as

\begin{eqnarray*} p_1&=&Re^{j\theta_c}\\ p_2&=&Re^{-j\theta_c} \end{eqnarray*}

where

\begin{eqnarray*} R&=&\sqrt{x_p^2 + y_p^2}\,>\,0\\ \theta_c&=&\tan^{-1}\left(\frac{y_p}{x_p}\right). \end{eqnarray*}

$ R$ is the pole radius, or distance from the origin in the $ z$-plane. As discussed in Chapter 8, we must have $ R<1$ for stability of the two-pole filter. The angles $ \pm\theta_c$ are the poles' respective angles in the $ z$ plane. The pole angle $ \theta _c$ corresponds to the pole frequency $ \omega_c$ via the relation

$\displaystyle \theta_c = \omega_c T = 2\pi f_c T $

where $ T$ denotes the sampling interval. See Chapter 8 for a discussion and examples of pole-zero plots in the complex $ z$-plane.

If $ R$ is sufficiently large (but less than 1 for stability), the filter exhibits a resonanceB.2 at radian frequency $ \omega_c = 2\pi f_c = \theta_c/T$. We may call $ \omega_c$ or $ f_c$ the center frequency of the resonator. Note, however, that the resonance frequency is not usually the precise frequency of peak-gain in a two-pole resonator (see Fig.B.9 on page [*]). The peak of the amplitude response is usually a little different because each pole sits on the other's ``skirt,'' which is slanted. (See §B.1.5 and §B.6 for an elaboration of this point.)

Using polar form for the (complex) poles, the two-pole transfer function can be expressed as

$\displaystyle H(z)$ $\displaystyle =$ $\displaystyle \frac{b_0}{(1-Re^{j\theta_c}z^{-1})(1-Re^{-j\theta_c}z^{-1})}$  
  $\displaystyle =$ $\displaystyle \frac{b_0}{1 - 2 R \cos(\theta_c)z^{-1}+ R^2 z^{-2}} \protect$ (B.1)

Comparing this to the transfer function derived from the difference equation, we may identify

\begin{eqnarray*} a_1 &=& -2R \cos(\theta_c)\\ a_2 &=& R^2. \protect \end{eqnarray*}

The difference equation can thus be rewritten as

$\displaystyle y(n) = b_0 x(n) + [2 R \cos(\theta_c)] y(n - 1) - R^2 y(n - 2). \protect$ (B.2)

Note that coefficient $ a_2$ depends only on the pole radius R (which determines damping) and is independent of the resonance frequency, while $ a_1$ is a function of both. As a result, we may retune the resonance frequency of the two-pole filter section by modifying $ a_1$ only.

The gain at the resonant frequency $ \omega=\omega_c$, is found by substituting $ z = e^{j\theta_c}=e^{j\omega_c T}$ into Eq.$ \,$(B.1) to get

$\displaystyle G(\omega_c) \isdef \left\vert H(e^{j\theta_c})\right\vert$ $\displaystyle =$ $\displaystyle \left\vert\frac{b_0}{(1-R)(1-Re^{-j2\theta_c})}\right\vert$  
  $\displaystyle =$ $\displaystyle \frac{\vert b_0\vert}{(1-R)\sqrt{1-2R\cos(2\theta_c)+R^2}} \protect$ (B.3)

See §B.6 for details on how the resonance gain (and peak gain) can be normalized as the tuning of $ \omega_c$ is varied in real time.

Since the radius of both poles is $ R$, we must have $ R<1$ for filter stability (§8.4). The closer $ R$ is to 1, the higher the gain at the resonant frequency $ \omega_c = 2\pi f_c$. If $ R=0$, the filter degenerates to the form $ H(z) = b_0$, which is a nothing but a scale factor. We can say that when the two poles move to the origin of the $ z$ plane, they are canceled by the two zeros there.

Resonator Bandwidth in Terms of Pole Radius

The magnitude $ R$ of a complex pole determines the damping or bandwidth of the resonator. (Damping may be defined as the reciprocal of the bandwidth.)

As derived in §8.5, when $ R$ is close to 1, a reasonable definition of 3dB-bandwidth $ B$ is provided by

$\displaystyle B$ $\displaystyle \isdef$ $\displaystyle - \frac{\ln(R)}{\pi T}$ (B.4)
$\displaystyle R$ $\displaystyle =$ $\displaystyle e^{- \pi B T} \protect$ (B.5)

where $ R$ is the pole radius, $ B$ is the bandwidth in Hertz (cycles per second), and $ T$ is the sampling interval in seconds.

Figure B.6 shows a family of frequency responses for the two-pole resonator obtained by setting $ b_0 = 1$ and varying $ R$. The value of $ \theta _c$ in all cases is $ \pi /4$, corresponding to $ f_c = f_s/8$. The analytic expressions for amplitude and phase response are

\begin{eqnarray*} G(\omega)\! &=& \!\frac{b_0}{\sqrt{[1 + a_1 \cos(\omega T) + a... ... + a_1 \cos(\omega T) + a_2 \cos(2\omega T)}\right]\qquad(b_0>0) \end{eqnarray*}

where $ a_1 = - 2R \cos(\theta_c)$ and $ a_2 = R^2$.

Figure B.6: Frequency response of the two-pole filter
$ y(n) = x(n) + 2R \cos (\theta _c) y(n - 1) - R^2 y(n - 2)$
with $ \theta _c$ fixed at $ \pi /4$ and for various values of $ R$. (a) Amplitude response. (b) Phase response.
\begin{figure}\input fig/kfig2p23.pstex_t \end{figure}

Two-Zero

The signal flow graph for the general two-zero filter is given in Fig.B.7, and the derivation of frequency response is as follows:

\fbox{ \begin{tabular}{rl} Difference equation: & $y(n) = b_0 x(n) + b_1 x(n-1) ... ...+ b_1 \cos(\omega T) + b_2 \cos(2\omega T)}\right]$ \end{tabular}\vspace{10pt} }

Figure B.7: Signal flow graph for the general two-zero filter
$ y(n) = b_0x(n) + b_1x(n - 1) + b_2x(n - 2)$.
\begin{figure}\input fig/twozero.pstex_t \end{figure}

As discussed in §5.1, the parameters $ b_1$ and $ b_2$ are called the numerator coefficients, and they determine the two zeros. Using the quadratic formula for finding the roots of a second-order polynomial, we find that the zeros are located at

$\displaystyle z = \frac{-b_1 \pm \sqrt{b_1^2 - 4 b_0 b_2}}{2b_0} $

If the zeros are real [ $ (b_1/2)^2 \geq b_2$], then the two-zero case reduces to two instances of our earlier analysis for the one-zero. Assuming the zeros to be complex, we may express the zeros in polar form as $ Re^{j\theta_c}$ and $ Re^{-j\theta_c}$, where $ \theta_c = \omega_c T = 2\pi f_c T$.

Forming a general two-zero transfer function in factored form gives

\begin{eqnarray*} H(z) &=& b_0 (1 - Re^{j\theta_c} z^{-1}) (1 - Re^{-j\theta_c} z^{-1})\\ &=& b_0 [1 - 2R\cos(\theta_c) z^{-1}+ R^2 z^{-2}] \end{eqnarray*}

from which we identify $ b_1/b_0 = - 2R \cos(\theta_c)$ and $ b_2/b_0 = R^2$, so that

$\displaystyle y(n) = b_0\{ x(n) - [2R \cos(\theta_c)]x(n - 1) + R^2 x(n - 2)\} $

is again the difference equation of the general two-zero filter with complex zeros. The frequency $ \omega$, is now viewed as a notch frequency, or antiresonance frequency. The closer R is to 1, the narrower the notch centered at $ \omega_c$.

The approximate relation between bandwidth and $ R$ given in Eq.$ \,$(B.5) for the two-pole resonator now applies to the notch width in the two-zero filter.

Figure B.8 gives some two-zero frequency responses obtained by setting $ b_0$ to 1 and varying $ R$. The value of $ \theta _c$, is again $ \pi /4$. Note that the response is exactly analogous to the two-pole resonator with notches replacing the resonant peaks. Since the plots are on a linear magnitude scale, the two-zero amplitude response appears as the reciprocal of a two-pole response. On a dB scale, the two-zero response is an upside-down two-pole response.

Figure B.8: Frequency response of the two-zero filter
$ y(n) = x(n) - 2R\cos (\theta _c) x(n - 1) + R^2 x(n - 2)$
with $ \theta _c$ fixed at $ \pi /4$ and for various values of $ R$. (a) Amplitude response. (b) Phase response.
\begin{figure}\input fig/kfig2p25.pstex_t \end{figure}

Complex Resonator

Normally when we need a resonator, we think immediately of the two-pole resonator. However, there is also a complex one-pole resonator having the transfer function

$\displaystyle H(z) = \frac{g}{1-pz^{-1}} \protect$ (B.6)

where $ p=Re^{j\theta_c}$ is the single complex pole, and $ g$ is a scale factor. In the time domain, the complex one-pole resonator is implemented as

$\displaystyle y(n) = g x(n) + p y(n-1). $

Since $ p$ is complex, the output $ y(n)$ is generally complex even when the input $ x(n)$ is real.

Since the impulse response is the inverse z transform of the transfer function, we can write down the impulse response of the complex one-pole resonator by recognizing Eq.$ \,$(B.6) as the closed-form sum of an infinite geometric series, yielding

$\displaystyle h(n) = u(n) g p^n, $

where, as always, $ u(n)$ denotes the unit step function:

$\displaystyle u(n) \isdef \left\{\begin{array}{ll} 1, & n\geq 0 \\ [5pt] 0, & n<0 \\ \end{array}\right. $

Thus, the impulse response is simply a scale factor $ g$ times the geometric sequence $ p^n$ with the pole $ p$ as its ``term ratio''. In general, $ p^n = R^n e^{j\omega_c nT}$ is a sampled, exponentially decaying sinusoid at radian frequency $ \omega_c=\theta_c/T$. By setting $ p$ somewhere on the unit circle to get

$\displaystyle p \isdef e^{j\omega_c T}, $

we obtain a complex sinusoidal oscillator at radian frequency $ \omega_c$ rad/sec. If we like, we can extract the real and imaginary parts separately to create both a sine-wave and a cosine-wave output:

\begin{eqnarray*} \mbox{re}\left\{h(n)\right\} &=& u(n) g \cos(\omega_c n T)\\ \mbox{im}\left\{h(n)\right\} &=& u(n) g \sin(\omega_c n T) \end{eqnarray*}

These may be called phase-quadrature sinusoids, since their phases differ by 90 degrees. The phase quadrature relationship for two sinusoids means that they can be regarded as the real and imaginary parts of a complex sinusoid.

By allowing $ g$ to be complex,

$\displaystyle g \isdef A e^{j\phi} $

we can arbitrarily set both the amplitude and phase of this phase-quadrature oscillator:

\begin{eqnarray*} \mbox{re}\left\{h(n)\right\} &=& u(n) A \cos(\omega_c n T + \p... ...mbox{im}\left\{h(n)\right\} &=& u(n) A \sin(\omega_c n T + \phi) \end{eqnarray*}

The frequency response of the complex one-pole resonator differs from that of the two-pole real resonator in that the resonance occurs only for one positive or negative frequency $ \omega_c$, but not both. As a result, the resonance frequency $ \omega_c$ is also the frequency where the peak-gain occurs; this is only true in general for the complex one-pole resonator. In particular, the peak gain of a real two-pole filter does not occur exactly at resonance, except when $ \theta_c \isdef \omega_c T = 0$, $ \pi/2$, or $ \pi $. See §B.6 for more on peak-gain versus resonance-gain (and how to normalize them in practice).

Two-Pole Partial Fraction Expansion

Note that every real two-pole resonator can be broken up into a sum of two complex one-pole resonators:

$\displaystyle H(z) = \frac{g}{(1 - p z^{-1}) (1 - \overline{p} z^{-1})} = \frac{g_1}{1-pz^{-1}} + \frac{g_2}{1-\overline{p}z^{-1}} \protect$ (B.7)

where $ g_1$ and $ g_2$ are constants (generally complex). In this ``parallel one-pole'' form, it can be seen that the peak gain is no longer equal to the resonance gain, since each one-pole frequency response is ``tilted'' near resonance by being summed with the ``skirt'' of the other one-pole resonator, as illustrated in Fig.B.9. This interaction between the positive- and negative-frequency poles is minimized by making the resonance sharper ( $ \left\vert p\right\vert\to1$), and by separating the pole frequencies $ 0\ll\angle p \ll \pi$. The greatest separation occurs when the resonance frequency is at one-fourth the sampling rate ( $ \angle p =\pi/2$). However, low-frequency resonances, which are by far the most common in audio work, suffer from significant overlapping of the positive- and negative-frequency poles.

Figure B.9: Frequency response (solid lines) of the two-pole resonator
$ H(z)=1/(1-2R\cos (\theta _c)z^{-1}+ R^2z^{-2})$,
for $ R=0.8$ and $ \theta _c = \pi /8$, overlaid with the frequency responses (dashed lines) of its positive- and negative-frequency complex one-pole components. Also marked (dashed lines) are the two resonance frequencies; the peak frequencies can be seen to lie slightly outside the resonance frequencies.
\includegraphics[width=\twidth ]{eps/tppfe}

To show Eq.$ \,$(B.7) is always true, let's solve in general for $ g_1$ and $ g_2$ given $ g$ and $ p$. Recombining the right-hand side over a common denominator and equating numerators gives

$\displaystyle g = g_1 - g_1 \overline{p}z^{-1}+ g_2 - g_2 pz^{-1} $

which implies

\begin{eqnarray*} g_1+g_2 &=& g\\ g_1 \overline{p} + g_2 p &=& 0. \end{eqnarray*}

The solution is easily found to be

\begin{eqnarray*} g_1 &=& g \frac{p}{2\mbox{im}\left\{p\right\}}\\ g_2 &=& -g \frac{\overline{p}}{2\mbox{im}\left\{p\right\}} \end{eqnarray*}

where we have assumed im$ \left\{p\right\}\neq 0$, as necessary to have a resonator in the first place.

Breaking up the two-pole real resonator into a parallel sum of two complex one-pole resonators is a simple example of a partial fraction expansion (PFE) (discussed more fully in §6.8).

Note that the inverse z transform of a sum of one-pole transfer functions can be easily written down by inspection. In particular, the impulse response of the PFE of the two-pole resonator (see Eq.$ \,$(B.7)) is clearly

$\displaystyle h(n) = g_1 p^n + g_2 \overline{p}^n,\qquad n=0,1,2,\ldots $

Since $ h(n)$ is real, we must have $ g_2=\overline{g}_1$, as we found above without assuming it. If $ \left\vert p\right\vert=1$, then $ h(n)$ is a real sinusoid created by the sum of two complex sinusoids spinning in opposite directions on the unit circle.

The BiQuad Section

The term ``biquad'' is short for ``bi-quadratic'', and is a common name for a two-pole, two-zero digital filter. The transfer function of the biquad can be defined as

$\displaystyle H(z) = g \frac{1 + \beta_1 z^{-1}+ \beta_2 z^{-2}}{1 + a_1 z^{-1}+ a_2 z^{-2}} \protect$ (B.8)

where $ g$ can be called the overall gain of the biquad. Since both the numerator and denominator of this transfer function are quadratic polynomials in $ z^{-1}$ (or $ z$), the transfer function is said to be ``bi-quadratic'' in $ z^{-1}$ (or $ z$).

As derived in §B.1.3, for real second-order polynomials having complex roots, it is often convenient to express the polynomial coefficients in terms of the radius $ R$ and angle $ \theta$ of the positive-frequency pole. For example, denoting the denominator polynomial by $ A(z)=1 + a_1 z^{-1}+ a_2 z^{-2}$, we have

$\displaystyle A(z) = \left(1 - Re^{j\theta}z^{-1}\right)\left(1 - Re^{-j\theta}z^{-1}\right) = 1 - 2R\cos(\theta)z^{-1}+ R^2z^{-2}. $

This representation is most often used for the denominator of the biquad, and we think of $ \theta$ as the resonance frequency (in radians per sample-- $ \theta=2\pi f_c T$, where $ f_c$ is the resonance frequency in Hz), and $ R$ determines the ``Q'' of the resonance (see §B.1.3). The numerator is less often represented in this way, but when it is, we may think of the zero-angle as the antiresonance frequency, and the zero-radius affects the depth and width of the antiresonance (or notch).

As discussed on page [*], a common setting for the zeros when making a resonator is to place one at $ z = 1$ (dc) and the other at $ z = -1$ (half the sampling rate), i.e., $ \beta_1=0$ and $ \beta_2=-1$ in Eq.$ \,$(B.8) above $ \Rightarrow B(z) = 1-z^{-2}= (1-z^{-1})(1+z^{-1})$. This zero placement normalizes the peak gain of the resonator if it is swept using the $ a_1$ parameter.

Using the shift theorem for z transforms, the difference equation for the biquad can be written by inspection of the transfer function as

\begin{eqnarray*} v(n) &=& g\, x(n) \\ y(n) &=& v(n) + \beta_1 v(n-1) + \beta_2 v(n-2) \\ & & \qquad - a_1 y(n-1) - a_2 y(n-2) . \end{eqnarray*}

where $ x(n)$ denotes the input signal sample at time $ n$, and $ y(n)$ is the output signal. This is the form that is typically implemented in software. It is essentially the direct-form I implementation. (To obtain the official direct-form I structure, the overall gain $ g$ must be not be pulled out separately, resulting in feedforward coefficients $ [g,g\beta_1,g\beta_2]$ instead. See Chapter 9 for more about filter implementation forms.)

Biquad Software Implementations

In matlab, an efficient biquad section is implemented by calling 

        outputsignal = filter(B,A,inputsignal);
where

\begin{eqnarray*} \texttt{B} &=& [g, g\beta_1, g\beta_2],\\ \texttt{A} &=& [1, a_1, a_2]. \end{eqnarray*}

A complete C++ class implementing a biquad filter section is included in the free, open-source Synthesis Tool Kit (STK) [15]. (See the BiQuad STK class.)

Figure B.10 lists an example biquad implementation in the C programming language.

Figure B.10: C code implementing a biquad filter section.
  
  typedef double *pp;  // pointer to array of length NTICK
  typedef word double; // signal and coefficient data type

  typedef struct _biquadVars {
      pp output;
      pp input;
      word s2;
      word s1;
      word gain;
      word a2;
      word a1;
      word b2;
      word b1;
  } biquadVars;

  void biquad(biquadVars *a)
  {
      int i;
      dbl A;
      word s0;
      for (i=0; i<NTICK; i++) {
          A = a->gain * a->input[i];
          A -= a->a1 * a->s1;
          A -= a->a2 * a->s2;
          s0 = A;
          A += a->b1 * a->s1;
          a->output[i] = a->b2 * a->s2 + A;
          a->s2 = a->s1;
          a->s1 = s0;
      }
  }

Allpass Filter Sections

The allpass filter passes all frequencies with equal gain. This is in contrast with a lowpass filter, which passes only low frequencies, a highpass which passes high-frequencies, and a bandpass filter which passes an interval of frequencies. An allpass filter may have any phase response. The only requirement is that its amplitude response be constant. Normally, this constant is $ \left\vert H(e^{j\omega T})\right\vert=1$.

From a physical modeling point of view, a unity-gain allpass filter models a lossless system in the sense that it preserves signal energy. Specifically, if $ x(n)$ denotes the input to an allpass filter $ H(z)$, and if $ y(n)$ denotes its output, then we have

$\displaystyle \sum_{n=-\infty}^\infty \left\vert x(n)\right\vert^2 = \sum_{n=-\infty}^\infty \left\vert y(n)\right\vert^2. \protect$ (B.9)

This equation says that the total energy out equals the total energy in. No energy was created or destroyed by the filter. All an allpass filter can do is delay the sinusoidal components of a signal by differing amounts.

Appendix C proves that Eq.$ \,$(B.9) holds if and only if

$\displaystyle \left\vert H(e^{j\omega T})\right\vert = 1, \quad \forall \omega. $

That is, a filter $ H(z)$ is lossless if and only if it is an allpass filter having a gain of $ 1$ at every frequency $ \omega$.

The Biquad Allpass Section

The general biquad transfer function was given in Eq.$ \,$(B.8) to be

$\displaystyle H(z) = g \frac{1 + \beta_1 z^{-1}+ \beta_2 z^{-2}}{1 + a_1 z^{-1}+ a_2 z^{-2}} \isdef \frac{B(z)}{A(z)}. $

To specialize this to a second-order unity-gain allpass filter, we require

$\displaystyle \left\vert H(e^{j\omega T})\right\vert = 1. $

It is easy to show that, given any monic denominator polynomial $ A(z)$, the numerator $ B(z)$ must be, in the real case,B.3

$\displaystyle B(z) = z^{-2}A(z^{-1}) = a_2 + a_1z^{-1}+ z^{-2}. $

Thus, to obtain an allpass biquad section, the numerator polynomial is simply the ``flip'' of the denominator polynomial. To obtain unity gain, we set $ g=a_2$, $ \beta_1 = a_1/a_2$, and $ \beta_2=1/a_2$.

In terms of the poles and zeros of a filter $ H(z)=B(z)/A(z)$, an allpass filter must have a zero at $ z=1/p$ for each pole at $ z=p$. That is if the denominator $ A(z)$ satisfies $ A(p)=0$, then the numerator polynomial $ B(z)$ must satisfy $ B(1/p)=0$. (Show this in the one-pole case.) Therefore, defining $ B(z) = A(1/z)$ takes care of this property for all roots of $ A(z)$ (all poles). However, since we prefer that $ B(z)$ be a polynomial in $ z^{-1}$, we define $ B(z) = z^{-N}A(1/z)$, where $ N$ is the order of $ A(z)$ (the number of poles). $ B(z)$ is then the flip of $ A(z)$.

For further discussion and examples of allpass filters (including muli-input, multi-output allpass filters), see Appendix C. Analog allpass filters are defined and discussed in §E.8.

Allpass Filter Design

There is a fairly large literature thread on the topic of allpass filter design. Generally, they fall into two main categories: parametric and nonparametric methods. Parametric methods can produce allpass filters with optimal group-delay characteristics [42,41]. Nonparametric methods, while suboptimal, can design very large-order allpass filters, and errors can usually be made arbitrarily small by increasing the order [100,70,1], [78, pp. 60,172]. In music applications, it is usually the case that the ``optimality'' criterion is unknown because it depends on aspects of sound perception (see, for example, [35,72]). As a result, perceptually weighted nonparametric methods can often outperform optimal parametric methods in terms of cost/performance. For a nonparametric method that can design very high-order allpass filters according to highly flexible criteria, see [1].

DC Blocker

The dc blocker is an indispensable tool in digital waveguide modeling [86] and other applications.B.4 It is often needed to remove the dc component of the signal circulating in a delay-line loop. It is also often an important tool in multi-track recording, where dc components in the various tracks can add up and overflow the mix.

The dc blocker is a small recursive filter specified by the difference equation

$\displaystyle y(n) = x(n) - x(n-1) + R\, y(n-1) $

where $ R$ is a parameter that is typically somewhere between $ 0.9$ and $ 1$ (for a 44.1 kHz sampling rate, $ R=0.995$ is good). The transfer function is

$\displaystyle H(z) = \frac{1-z^{-1}}{1-Rz^{-1}}. \protect$ (B.10)

Thus, there is a zero at dc ($ z = 1$) and a pole near dc at $ z=R$. Far away from dc, the pole and zero approximately cancel each other. (Recall the graphical method for determining frequency response magnitude described in Chapter 8.)

DC Blocker Frequency Response

Figure B.11 shows the frequency response of the dc blocker for several values of $ R$. The same plots are given over a log-frequency scale in Fig.B.12. The corresponding pole-zero diagrams are shown in Fig.B.13. As $ R$ approaches $ 1$, the notch at dc gets narrower and narrower. While this may seem ideal, there is a drawback, as shown in Fig.B.14 for the case of $ R=0.9$: The impulse response duration increases as $ R\to 1$. While the ``tail'' of the impulse response lengthens as $ R$ approaches 1, its initial magnitude decreases. At the limit, $ R=1$, the pole and zero cancel at all frequencies, the impulse response becomes an impulse, and the notch disappears.

Figure B.11: Frequency response overlays for the dc blocker defined by $ H(z) = (1-z^{-1})/(1-Rz^{-1})$ for various values of pole radius $ R$. (a) Amplitude response. (b) Phase response.
\includegraphics[width=\twidth ]{eps/dcblockerfr}

Figure B.12: Log-frequency response overlays for the dc blocker defined by $ H(z) = (1-z^{-1})/(1-Rz^{-1})$ for various values of pole radius $ R$. (a) Amplitude response. (b) Phase response.
\includegraphics[width=\twidth ]{eps/dcblockerfrlf}

Figure B.13: Pole-zero diagram overlays for the dc blocker defined by $ H(z) = (1-z^{-1})/(1-Rz^{-1})$ for various values of pole radius $ R$.
\includegraphics[width=\twidth]{eps/dcblockerpz}

Figure B.14: Impulse response of the dc blocker defined by $ H(z) = (1-z^{-1})/(1-0.9z^{-1})$.
\includegraphics[width=\twidth ]{eps/dcblockerir}

Note that the amplitude response in Fig.B.11a and Fig.B.12a exceeds 1 at half the sampling rate. This maximum gain is given by $ H(-1)=2/(1+R)$. In applications for which the gain must be bounded by 1 at all frequencies, the dc blocker may be scaled by the inverse of this maximum gain to yield

\begin{eqnarray*} H(z) &=& g\frac{1-z^{-1}}{1-Rz^{-1}}\\ y(n) &=& g[x(n) - x(n-1)] + R\, y(n-1), \quad\hbox{where}\\ g &\isdef & \frac{1+R}{2}. \end{eqnarray*}

DC Blocker Software Implementations

In plain C, the difference equation for the dc blocker could be written as follows: 

  y = x - xm1 + 0.995 * ym1;
  xm1 = x;
  ym1 = y;
Here, x denotes the current input sample, and y denotes the current output sample. The variables xm1 and ym1 hold once-delayed input and output samples, respectively (and are typically initialized to zero). In this implementation, the pole is fixed at $ R=0.995$, which corresponds to an adaptation time-constant of approximately $ 1/(1-R) = 200$ samples. A smaller $ R$ value allows faster tracking of ``wandering dc levels'', but at the cost of greater low-frequency attenuation.

A complete C++ class implementing a dc blocking filter is included in the free, open-source Synthesis Tool Kit (STK) [15]. (See the DCBlock STK class.)

For a discussion of issues and solutions related to fixed-point implementations, see [7].

Low and High Shelving Filters

The analog transfer function for a low shelf is given by [103]

$\displaystyle H(s) \;=\; 1 + \frac{B_0\omega_1}{s+\omega_1} \;=\; \frac{s+\omega_1(B_0+1)}{s+\omega_1} \;\isdef \; \frac{s+\omega_z}{s+\omega_1} $

where $ B_0$ is the dc boost amount (at $ s=0$), and the high-frequency gain ($ s=\infty$) is constrained to be $ 1$. The transition frequency dividing low and high frequency regions is $ \omega_1$. See Appendix E for a development of $ s$-plane analysis of analog (continuous-time) filters.

A high shelf is obtained from a low shelf by the conformal mapping $ s \leftarrow 1/s$, which interchanges high and low frequencies, i.e.,

$\displaystyle H(s) \;=\; 1 + \frac{B_\pi\omega_1}{\frac{1}{s}+\omega_1} \;=\; ... ...{\omega_z}{\omega_1} \cdot \frac{s + \frac{1}{\omega_z}}{s+\frac{1}{\omega_1}} $

In this case, the dc gain is 1 and the high-frequency gain approaches $ 1+B_\pi = \omega_z/\omega_1$.

To convert these analog-filter transfer functions to digital form, we apply the bilinear transform:

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

where $ T$ denotes the sampling interval in seconds.B.5

Low and high shelf filters are typically implemented in series, and are typically used to give a little boost or cut at the extreme low or high end (of the spectrum), respectively. To provide a boost or cut near other frequencies, it is necessary to go to (at least) a second-order section, often called a ``peaking equalizer,'' as described in §B.5 below.

Exercise

Perform the bilinear transform defined above and calculate the coefficients of a first-order digital low shelving filter. Find the pole and zero as a function of $ B_0$, $ \omega_1$, and $ T$. Set $ z = 1$ and verify that you get a gain of $ 1+B_0$. Set $ z = -1$ and verify that you get a gain of 1 there.

Peaking Equalizers

A peaking equalizer filter section provides a boost or cut in the vicinity of some center frequency. It may also be called a parametric equalizer section. The gain far away from the boost or cut is unity, so it is convenient to combine a number of such sections in series. Additionally, a high and/or low shelf (§B.4 above) are nice to include in series with one's peaking eq sections.

The analog transfer function for a peak filter is given by [103,5,6]

$\displaystyle H(s) = 1 + H_R(s) $

where $ H_R(s)$ is a two-pole resonator:

$\displaystyle H_R(s) \isdef R\cdot \frac{Bs}{s^2+Bs + 1} $

The transfer function can be written in the normalized form [103]

$\displaystyle H(s) = \frac{s^2 + gBs + 1}{s^2+Bs + 1}, $

where $ g$ is approximately the desired gain at the boost (or cut), and $ B$ is the desired bandwidth. When $ g>1$, a boost is obtained at frequency $ \omega=1$. For $ g<1$, a cut filter is obtained at that frequency. In particular, when $ g=0$, there are infinitely deep notches at $ \omega=\pm 1$, and when $ g=1$, the transfer function reduces to $ H(s)=1$ (no boost or cut). The parameter $ B$ controls the width of the boost or cut.

It is easy to show that both zeros and both poles are on the unit circle in the left-half $ s$ plane, and when $ g<1$ (a ``cut''), the zeros are closer to the $ j\omega$ axis than the poles.

Again, the bilinear transform can be used to convert the analog peaking equalizer section to digital form.

Figure B.15 gives a matlab listing for a peaking EQ section. Figure B.16 shows the resulting plot for an example call:

boost(2,0.25,0.1);
The frequency-response utility myfreqz, listed in Fig.7.1, can be substituted for freqz.

Figure B.15: Matlab function for designing (and optionally testing) a peaking equalizer section.
  
function [B,A] = boost(gain,fc,bw,fs);
%BOOST - Design a boost filter at given gain, center
%        frequency fc, bandwidth bw, and sampling rate fs
%        (default = 1).
%
% J.O. Smith 11/28/02
% Reference: Zolzer: Digital Audio Signal Processing, p. 124

if nargin<4, fs = 1; end
if nargin<3, bw = fs/10; end

Q = fs/bw;
wcT = 2*pi*fc/fs;

K=tan(wcT/2);
V=gain;

b0 =  1 + V*K/Q + K^2;
b1 =  2*(K^2 - 1);
b2 =  1 - V*K/Q + K^2;
a0 =  1 + K/Q + K^2;
a1 =  2*(K^2 - 1);
a2 =  1 - K/Q + K^2;
A = [a0 a1 a2] / a0;
B = [b0 b1 b2] / a0;

if nargout==0
  figure(1);
  freqz(B,A);
  title('Boost Frequency Response')
end

Figure B.16: Frequency response of a second-order peaking equalizer section tuned to give a 6 dB peak of width $ \approx f_s/10$ at center frequency $ f_s/4$.
\includegraphics[width=\twidth]{eps/tboost}

A Faust implementation of the peaking equalizer is available as the function peak_eq in filter.lib distributed with Faust (Appendix K) starting with version 0.9.9.4k-par.

Time-Varying Two-Pole Filters

It is quite common to want to vary the resonance frequency of a resonator in real time. This is a special case of a tunable filter. In the pre-digital days of analog synthesizers, filter modules were tuned by means of control voltages, and were thus called voltage-controlled filters (VCF). In the digital domain, control voltages are replaced by time-varying filter coefficients. In the time-varying case, the choice of filter structure has a profound effect on how the filter characteristics vary with respect to coefficient variations. In this section, we will take a look at the time-varying two-pole resonator.

Evaluating the transfer function of the two-pole resonator (Eq.$ \,$(B.1)) at the point $ e^{j\theta_c}$ on the unit circle (the filter's resonance frequency $ \omega_c=\theta_c/T$) yields a gain at resonance equal to

$\displaystyle H(e^{j\theta_c})$ $\displaystyle =$ $\displaystyle \frac{b_0}{(1-Re^{j\theta_c}e^{-j\theta_c})(1-Re^{-j\theta_c}e^{-j\theta_c})}$  
  $\displaystyle =$ $\displaystyle \frac{b_0}{1-R}\cdot \frac{1}{1-Re^{-j2\theta_c}} \protect$ (B.11)

For simplicity, let $ b_0 = 1$ in what follows. In the special cases $ \theta_c=0$ (resonance at dc) and $ \theta_c=\pi$ (resonance at $ f=f_s/2$), we have

$\displaystyle H(\pm1) = \frac{1}{(1-R)^2} \protect$ (B.12)

Since $ R$ is real, we have already found the gain (amplitude response) at a dc or $ f_s/2$ resonance:

$\displaystyle G(0) = G(\pi/T) \isdef \left\vert H(\pm1)\right\vert = \left\vert\frac{1}{(1-R)^2}\right\vert = \frac{1}{(1-R)^2} $

In the middle frequency between dc and $ f_s/2$, $ \theta_c = \omega_c T = \pi/2$, Eq.$ \,$(B.11) with $ b_0 = 1$ becomes

$\displaystyle H(j) = \frac{1}{(1-Re^{j2\frac{\pi}{2}})(1-R)} = \frac{1}{(1+R)(1-R)} = \frac{1}{1-R^2} $

and, since $ H(j)$ is real and positive, it coincides with the amplitude response, i.e., $ H(j)=G(\pi/2)=1/(1-R^2)$.

An important fact we can now see is that the gain at resonance depends markedly on the resonance frequency. In particular, the ratio of the two cases just analyzed is

$\displaystyle \left\vert\frac{H(1)}{H(j)}\right\vert = \frac{1-R^2}{(1-R)^2} = ... ...R}{1-R} = \frac{\hbox{maximum resonance gain}}{\hbox{minimum resonance gain}} $

We did not show that resonance gain is maximized at $ e^{j\theta_c}=\pm 1$ and minimized at $ e^{j\theta_c}=\pm j$, but this is straightforward to show, and strongly suggested by Fig.B.17 (and Fig.B.9).

Note that the ratio of the dc resonance gain to the $ f_s/4$ resonance gain is unbounded! The sharper the resonance (the closer $ R$ is to 1), the greater the disparity in the gain.

Figure B.17 illustrates a number of resonator frequency responses for the case $ R=0.99$. (Resonators in practice may use values of $ R$ even closer to 1 than this--even the case $ R=1$ is used for making recursive digital sinusoidal oscillators [90].) For resonator tunings at dc and $ f_s/2$, we predict the resonance gain to be $ 20\log_{10}[1/(1-R)^2] = -20\log_{10}[(1-0.99)^2] = -40\log_{10}(0.01) = 80$ dB, and this is what we see in the plot. When the resonance is tuned to $ f_s/4$, the gain drops well below 40 dB. Clearly, we will need to compensate this gain variation when trying to use the two-pole digital resonator as a tunable filter.

Figure B.17: Frequency response overlays for the two-pole resonator $ H(z)=1/(1-2R\cos (\theta _c)z^{-1}+ R^2z^{-2})$, for $ R=0.99$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines.
\includegraphics[width=\twidth ]{eps/resgain}

Figure B.18 shows the same type of plot for the complex one-pole resonator $ H(z)=1/(1-Re^{j\theta _c}z^{-1})$, for $ R=0.99$ and 10 values of $ \theta _c$. In this case, we expect the frequency response evaluated at the center frequency to be $ H(e^{j\omega_c T}) =1/(1-Re^{j\theta_c}e^{-j\theta_c})=1/(1-R)$. Thus, the gain at resonance for the plotted example is $ 1/(1-0.99)=100=40$ db for all tunings. Furthermore, for the complex resonator, the resonance gain is also exactly equal to the peak gain.

Figure B.18: Frequency response overlays for the one-pole complex resonator $ H(z)=1/(1-Re^{j\theta _c}z^{-1})$, for $ R=0.99$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines.
\includegraphics[width=\twidth ]{eps/cresgain}

Normalizing Two-Pole Filter Gain at Resonance

The question we now pose is how to best compensate the tunable two-pole resonator of §B.1.3 so that its peak gain is the same for all tunings. Looking at Fig.B.17, and remembering the graphical method for determining the amplitude response,B.6 it is intuitively clear that we can help matters by adding two zeros to the filter, one near dc and the other near $ f_s/2$. A zero exactly at dc is provided by the term $ (1-z^{-1})$ in the transfer function numerator. Similarly, a zero at half the sampling rate is provided by the term $ (1+z^{-1})$ in the numerator. The series combination of both zeros gives the numerator $ B(z)=(1-z^{-1})(1+z^{-1})=1-z^{-2}$. The complete second-order transfer function then becomes

$\displaystyle H(z) = \frac{B(z)}{A(z)} = \frac{1 - z^{-2}}{1-2R\cos(\theta_c)z^{-1}+ R^2z^{-2}} $

corresponding to the difference equation

$\displaystyle y(n) = x(n) - x(n-2) + [2R\cos(\theta_c)] y(n-1) - R^2 y(n-2). \protect$ (B.13)

Checking the gain for the case $ \theta_c=\pi/2$, we have

\begin{eqnarray*} H(\pm1) &=& 0\\ \left.H(e^{j\theta_c})\right\vert _{\theta_c=\pi/2} &=& \frac{2}{1-R^2} \end{eqnarray*}

which is better behaved, but now the response falls to zero at dc and $ f_s/2$ rather than being heavily boosted, as we found in Eq.$ \,$(B.12).

Constant Resonance Gain

It turns out it is possible to normalize exactly the resonance gain of the second-order resonator tuned by a single coefficient [89]. This is accomplished by placing the two zeros at $ z=\pm \sqrt{R}$, where $ R$ is the radius of the complex-conjugate pole pair . The transfer function numerator becomes $ B(z)=(1-\sqrt{R}z^{-1})(1+\sqrt{R}z^{-1}) = (1-Rz^{-2})$, yielding the total transfer function

$\displaystyle H(z) = \frac{B(z)}{A(z)} = \frac{1 - Rz^{-2}}{1-2R\cos(\theta_c)z^{-1}+ R^2z^{-2}} $

which corresponds to the difference equation

$\displaystyle y(n) = x(n) - R\, x(n-2) + [2R\cos(\theta_c)] y(n-1) - R^2 y(n-2). $

We see there is one more multiply-add per sample (the term $ -R x(n-2)$) relative to the unnormalized two-pole resonator of Eq.$ \,$(B.13). The resonance gain is now

\begin{eqnarray*} H(e^{j\theta_c}) &=& \frac{1 - R e^{-j2\theta_c}}{1-2R\cos(\th... ...\theta_c}}{(1-R) - (1-R)Re^{-j2\theta_c}}\\ &=& \frac{1}{1 - R} \end{eqnarray*}

Thus, the gain at resonance is $ 1/(1-R)$ for all resonance tunings $ \theta _c$.

Figure B.19 shows a family of amplitude responses for the constant resonance-gain two-pole, for various values of $ \theta _c$ and $ R=0.99$. We see an excellent improvement in the regularity of the amplitude response as a function of tuning.

Figure: Frequency response overlays for the constant resonance-gain two-pole filter $ H(z)=(1-R)(1-Rz^{-2})/(1-2R\cos(\theta_c)z^{-1}+R^2z^{-2})$, for $ R=0.99$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines.
\includegraphics[width=\twidth ]{eps/cgresgain}

Peak Gain Versus Resonance Gain

While the constant resonance-gain filter is very well behaved, it is not ideal, because, while the resonance gain is perfectly normalized, the peak gain is not. The amplitude-response peak does not occur exactly at the resonance frequencies $ \omega T=\pm\theta_c$ except for the special cases $ \theta_c=0$, $ \pm\pi/2$, and $ \pi $. At other resonance frequencies, the peak due to one pole is shifted by the presence of the other pole. When $ R$ is close to 1, the shifting can be negligible, but in more damped resonators, e.g., when $ R<0.9$, there can be a significant difference between the gain at resonance and the true peak gain.

Figure B.20 shows a family of amplitude responses for the constant resonance-gain two-pole, for various values of $ \theta _c$ and $ R=0.9$. We see that while the gain at resonance is exactly the same in all cases, the actual peak gain varies somewhat, especially near dc and $ f_s/2$ when the two poles come closest together. A more pronounced variation in peak gain can be seen in Fig.B.21, for which the pole radii have been reduced to $ R=0.5$.

Figure: Frequency response overlays for the constant resonance-gain two-pole filter $ H(z)=(1-R)(1-Rz^{-2})/(1-2R\cos(\theta_c)z^{-1}+R^2z^{-2})$, for $ R=0.9$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines.
\includegraphics[width=\twidth ]{eps/cgresgaindamped}

Figure: Frequency response overlays for the constant resonance-gain two-pole filter $ H(z)=(1-R)(1-Rz^{-2})/(1-2R\cos(\theta_c)z^{-1}+R^2z^{-2})$, for $ R=0.5$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines. Note the more pronounced variation in peak gain (the resonance gain does not vary).
\includegraphics[width=\twidth ]{eps/cgresgaindampedp5}

Constant Peak-Gain Resonator

It is surprisingly easy to normalize exactly the peak gain in a second-order resonator tuned by a single coefficient [94]. The filter structure that accomplishes this is the one we already considered in §B.6.1:

$\displaystyle H(z) = \frac{B(z)}{A(z)} = \frac{1 - z^{-2}}{1-2R\cos(\theta_c)z^{-1}+ R^2z^{-2}} $ (B.14)

That is, the two-pole resonator normalized by zeros at $ z=\pm 1$ has the constant peak-gain property when it has resonant peaks in its response at all. Note, however, that the peak-gain frequency and the pole-resonance frequency (cf. §B.6.3), are generally two different things, as elaborated below. This structure has the added bonus that its difference equation requires only one more addition relative to the unnormalized two-pole resonator, and no new multiply.

Real-time audio ``plugins'' based on the constant-peak-gain resonator are developed in Appendix K.

The peak gain is $ 2/(1-R^2)$, so multiplying the transfer function by $ (1-R^2)/2$ normalizes the peak gain to one for all tunings. It can also be shown [94] that the peak gain coincides with the variance gain when the resonator is driven by white noise. That is, if the variance of the driving noise is $ \sigma^2$, the variance of the noise at the resonator output is $ 2\sigma^2/(1-R^2)$. Therefore, scaling the resonator input by $ g=\sqrt{(1-R^2)/2}$ will normalize the resonator such that the output signal power equals the input signal power when the input signal is white noise.

Frequency response overlays for the constant-peak-gain resonator are shown in Fig.B.23 ($ R=0.99$), Fig.B.20 ($ R=0.9$), and Fig.B.21 ($ R=0.5$). While the peak frequency may be far from the resonance tuning in the more heavily damped examples, the peak gain is always normalized to unity. The normalized radian frequency $ \psi\in[-\pi,\pi]$ at which the peak gain occurs is related to the pole angle $ \theta_c\in[-\pi,\pi]$ by [94]

$\displaystyle \cos(\theta_c) = \frac{1+R^2}{2R}\cos(\psi). \protect$ (B.15)

When the right-hand side of the above equation exceeds 1 in magnitude, there is no (real) solution for the pole frequency $ \theta _c$. This happens, for example, when $ R$ is less than 1 and $ \psi$ is too close to 0 or $ \pi $. Conversely, given any pole angle $ \theta_c\in(0,\pi)$, there always exists a solution for the peak frequency $ \psi = \arccos[2R\cos(\omega_c)/(1+R^2)]$, since $ \vert 2R/(1+R^2)\vert\leq1$ when $ R\in[0,1]$. However, when $ R$ is small, the peak frequency can be far from the pole resonance frequency, as shown in Fig.B.22.

Figure B.22: Upper and lower peak-gain frequency limits as a function of pole radius $ R$ ( $ \arccos[\pm2R/(1+R^2)]$).
\includegraphics[width=\twidth ]{eps/psivsthetac}

Thus, $ R$ must be close to 1 to obtain a resonant peak near dc (a case commonly needed in audio work) or half the sampling rate (rarely needed in practice). When $ R$ is much less than 1, the peak frequency $ \psi$ cannot leave a small interval near one-fourth the sampling rate, as can be seen at the far left in Fig.B.22.

Figure B.22 predicts that for $ R=0.5$, the lowest peak-gain frequency should be around $ \psi\geq 0.515$ radian per sample. Figure B.21 agrees with this prediction.

As Figures B.23 through B.25 show, the peak gain remains constant even at very low and very high frequencies, to the extent they are reachable for a given $ R$. The zeros at dc and $ f_s/2$ preclude the possibility of peaks at exactly those frequencies, but for $ R$ near 1, we can get very close to having a peak at dc or $ f_s/2$, as shown in Figures B.19 and B.20.

Figure: Frequency response overlays for the constant peak-gain two-pole filter $ H(z)=[(1-R^2)/2](1-z^{-2})/(1-2R\cos(\theta_c)z^{-1}+R^2z^{-2})$, for $ R=0.99$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines.
\includegraphics[width=\twidth ]{eps/cpgresgain}

Figure: Frequency response overlays for the constant peak-gain two-pole filter $ H(z)=[(1-R^2)/2](1-z^{-2})/(1-2R\cos(\theta_c)z^{-1}+R^2z^{-2})$, for $ R=0.9$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines.
\includegraphics[width=\twidth ]{eps/cpgresgaindamped}

Figure: Frequency response overlays for the constant peak-gain two-pole filter $ H(z)=[(1-R^2)/2](1-z^{-2})/(1-2R\cos(\theta_c)z^{-1}+R^2z^{-2})$, for $ R=0.5$ and 10 values of $ \theta _c$ uniformly spaced from 0 to $ \pi $. The 5th case is plotted using thicker lines.
\includegraphics[width=\twidth ]{eps/cpgresgaindampedp5}

Four-Pole Tunable Lowpass/Bandpass Filters

As a practical note, it is worth mentioning that in popular analog synthesizers (both real and virtualB.7), VCFs are typically fourth order rather than second order as we have studied here. Perhaps the best known VCF is the Moog VCF. The four-pole Moog VCF is configured to be a lowpass filter with an optional resonance near the cut-off frequency. When the resonance is strong, it functions more like a resonator than a lowpass filter. Various methods for digitizing the Moog VCF are described in [95]. It turns out to be nontrivial to preserve all desirable properties of the analog filter (such as frequency response, order, and control structure), when translated to digital form by standard means.

Elementary Filter Problems

See http://ccrma.stanford.edu/~jos/filtersp/Elementary_Filter_Problems.html.

Next Section:
Allpass Filters
Previous Section:
Background Fundamentals