DSPRelated.com
Free Books
   Free Books    Introduction to Digital Filters  

Analysis of a Digital Comb Filter

In Chapter 1, we extensively analyzed the simplest lowpass filter, $ y(n) = x(n) + x(n - 1)$ from a variety of points of view. This served to introduce many important concepts necessary for understanding digital filters. In Chapter 2, we analyzed the same filter using the matlab programming language. This chapter takes the next step by analyzing a more practical example, the digital comb filter, from start to finish using the analytical tools developed in later chapters. Consider this a ``practical motivation'' chapter--i.e., its purpose is to introduce and illustrate the practical utility of tools for filter analysis before diving into a systematic development of those tools.

Suppose you look up the documentation for a ``comb filter'' in a software package you are using, and you find it described as follows: 

    out(n) = input(n) + feedforwardgain * input(n-delay1)
                      - feedbackgain * out(n-delay2)
Does this tell you everything you need to know? Well, it does tell you exactly what is implemented, but to fully understand it, you must be able to predict its audible effects on sounds passing through it. One helpful tool for this purpose is a plot of the frequency response of the filter. Moreover, if delay1 or delay2 correspond to more than a a few milliseconds of time delay, you probably want to see its impulse response as well. In some situations, a pole-zero diagram can give valuable insights.

As a preview of things to come, we will analyze and evaluate the above example comb filter rather thoroughly. Don't worry about understanding the details at this point--just follow how the analysis goes and try to intuit the results. It will also be good to revisit this chapter later, after studying subsequent chapters, as it provides a concise review of the main topics covered. If you already fully understand the analyses illustrated in this chapter, you might consider skipping ahead to the discussion of specific audio filters in Appendix B, §B.4.

Difference Equation

We first write the comb filter specification in a more ``mathematical'' form as a digital filter difference equation:

$\displaystyle y(n) = x(n) + g_1\, x(n-M_1) - g_2\, y(n-M_2), \protect$ (4.1)

where $ x(n)$ denotes the $ n$th sample of the input signal, $ y(n)$ is the output at time $ n$, and $ g_1$ and $ g_2$ denote the feedforward and feedback coefficients, respectively. The signal flow graph is shown in Fig.3.1.

Signal Flow Graph

Figure 3.1 shows the signal flow graph (or system diagram) for the class of digital filters we are considering (digital comb filters). The symbol ``$ z^{-M}$'' means a delay of $ M$ samples (always an integer here). This notation for a pure delay will make more sense after we derive the filter transfer function in §3.7.

Figure 3.1: signal flow graph for digital filters having difference equations of the form $ y(n) = x(n) + g_1\, x(n-M_1) - g_2\, y(n-M_2)$ (digital comb filters).
\begin{figure}\input fig/esd.pstex_t \end{figure}

Software Implementation in Matlab

In Matlab or Octave, this type of filter can be implemented using the filter function. For example, the following matlab4.1 code computes the output signal y given the input signal x for a specific example comb filter: 

g1 = (0.5)^3;        % Some specific coefficients
g2 = (0.9)^5;
B = [1 0 0 g1];      % Feedforward coefficients, M1=3
A = [1 0 0 0 0 g2];  % Feedback coefficients, M2=5
N = 1000;            % Number of signal samples
x = rand(N,1);       % Random test input signal
y = filter(B,A,x);   % Matlab and Octave compatible
The example coefficients, $ g1 = 0.5^3 = 0.125$ and $ g2 = 0.9^5 = 0.59049$, are chosen to place all filter zeros at radius $ 0.5$ and all filter poles at radius $ 0.9$ in the complex $ z$ plane (as we shall see below).

The matlab filter function carries out the following computation for each element of the y array:

$\displaystyle \texttt{y(n)}$ $\displaystyle =$ $\displaystyle \sum_{\texttt{k=0}}^\texttt{NB-1} \texttt{B(k+1) * x(n-k)}$ (4.2)
  $\displaystyle -$ $\displaystyle \sum_{\texttt{k=1}}^\texttt{NA-1} \texttt{A(k+1) * y(n-k)} \protect$ (4.3)

for $ \texttt{n}=1,2,\dots,\texttt{N}$, where NA = length(A) and NB = length(B). Note that the indices of x and y can go negative in this expression. By default, such terms are replaced by zero. However, the filter function has an optional fourth argument for specifying the initial state of the filter, which includes past input and past output samples seen by the filter. This argument is used to forward the filter's state across successive blocks of data:
[y1,state] = filter(B,A,x1);       % filter 1st block x1
[y2,state] = filter(B,A,x2,state); % filter 2nd block x2
...

Sample-Level Implementation in Matlab

For completeness, a direct matlab implementation of the built-in filter function (Eq.$ \,$(3.3)) is given in Fig.3.2. While this code is useful for study, it is far slower than the built-in filter function. As a specific example, filtering $ 10,000$ samples of data using an order 100 filter on a 900MHz Athlon PC required 0.01 seconds for filter and 10.4 seconds for filterslow. Thus, filter was over a thousand times faster than filterslow in this case. The complete test is given in the following matlab listing:

x = rand(10000,1); % random input signal
B = rand(101,1);   % random coefficients
A = [1;0.001*rand(100,1)]; % random but probably stable
tic; yf=filter(B,A,x); ft=toc
tic; yfs=filterslow(B,A,x); fst=toc
The execution times differ greatly for two reasons:
  1. recursive feedback cannot be ``vectorized'' in general, and
  2. built-in functions such as filter are written in C, precompiled, and linked with the main program.

Figure 3.2: Matlab function for implementing a digital filter directly. Do not use this in practice because it is much slower than the built-in filter routine.
  
function [y] = filterslow(B,A,x)
% FILTERSLOW: Filter x to produce y = (B/A) x .
%       Equivalent to 'y = filter(B,A,x)' using
%       a slow (but tutorial) method.

NB = length(B);
NA = length(A);
Nx = length(x);

xv = x(:); % ensure column vector

% do the FIR part using vector processing:
v = B(1)*xv;
if NB>1
  for i=2:min(NB,Nx)
    xdelayed = [zeros(i-1,1); xv(1:Nx-i+1)];
    v = v + B(i)*xdelayed;
  end;
end; % fir part done, sitting in v

% The feedback part is intrinsically scalar,
% so this loop is where we spend a lot of time.
y = zeros(length(x),1); % pre-allocate y
ac = - A(2:NA);
for i=1:Nx, % loop over input samples
  t=v(i);   % initialize accumulator
  if NA>1,
    for j=1:NA-1
      if i>j,
        t=t+ac(j)*y(i-j);
      %else y(i-j) = 0
      end;
    end;
  end;
  y(i)=t;
end;

y = reshape(y,size(x)); % in case x was a row vector

Software Implementation in C++

Figure 3.3 gives a C++ program for implementing the same filter using the Synthesis Tool Kit (STK). The example writes a sound file containing white-noise followed by filtered-noise from the example filter.

Figure: C++ main program for computing filtered noise using the example digital filter. The Synthesis Tool Kit (STK), version 4.2, supplies the Noise and Filter classes used.
  
/* main.cpp - filtered white noise example */

#include "Noise.h"  /* Synthesis Tool Kit (STK) class */
#include "Filter.h" /* STK class */
#include "FileWvOut.h"  /* STK class */
#include <cmath>   /* for pow() */
#include <vector>

int main()
{
  Noise *theNoise = new Noise(); // Noise source

  /* Set up the filter */
  StkFloat bCoefficients[5] = {1,0,0,pow(0.5,3)};
  std::vector<StkFloat> b(bCoefficients, bCoefficients+5);
  StkFloat aCoefficients[7] = {1,0,0,0,0,pow(0.9,5)};
  std::vector<StkFloat> a(aCoefficients, aCoefficients+7);
  Filter *theFilter = new Filter;
  theFilter->setNumerator(b);
  theFilter->setDenominator(a);

  FileWvOut output("main");  /* write to main.wav */

  /* Generate one second of white noise */
  StkFloat amp = 0.1; // noise amplitude
  for (long i=0;i<SRATE;i++)   {
    output.tick(amp*theNoise->tick());
  }

  /* Generate filtered noise for comparison */
  for (long i=0;i<SRATE;i++)   {
    output.tick(theFilter->tick(amp*theNoise->tick()));
  }
}

Software Implementation in Faust

The Faust language for signal processing is introduced in Appendix K. Figure 3.4 shows a Faust program for implementing our example comb filter. As illustrated in Appendix K, such programs can be compiled to produce LADSPA or VST audio plugins, or a Pure Data (PD) plugin, among others.

Figure 3.4: Faust main program implementing the example digital filter. (Tested in Faust version 0.9.9.2a2.)
  
/* GUI Controls */
g1  = hslider("feedforward gain", 0.125, 0, 1, 0.01);
g2  = hslider("feedback gain", 0.59049, 0, 1, 0.01);

/* Signal Processing */
process = firpart : + ~ feedback
with {
  firpart(x) = x + g1 * x''';
  feedback(v) = 0 - g2 * v'''';
};
As discussed in Appendix K, a prime (') denotes delaying a signal by one sample, and a tilde (~) denotes feedback. A colon (:) simply indicates a connection in series. The feedback signal v is delayed only four samples instead of five because there is a free ``pipeline delay'' associated with the feedback loop itself.

Faust's -svg option causes a block-diagram to be written to disk for each Faust expression (as further discussed in Appendix K). The block diagram for our example comb filter is shown in Fig.3.5.

Figure 3.5: Block diagram generated by the Faust -svg option.
\includegraphics[width=\twidth]{eps/faustcfbd}

Compiling the Faust code in Fig.3.4 for LADSPA plugin format produces a plugin that can be loaded into ``JACK Rack'' as depicted in Fig.3.6.

Figure 3.6: JACK Rack screenshot for the example comb filter.
\includegraphics[width=4in]{eps/faustcfjr}

At the risk of belaboring this mini-tour of filter embodiments in common use, Fig.3.7 shows a screenshot of a PD test patch for the PD plugin generated from the Faust code in Fig.3.4.

Figure: Pure Data (PD) screenshot for a test patch exercising a PD plugin named cf.pd that was generated automatically from the Faust code in Fig.3.4 using the faust program and faust2pd script (see Appendix K for details).
\includegraphics[width=3.5in]{eps/faustcfpd}

By the way, to change the Faust example of Fig.3.4 to include its own driving noise, as in the STK example of Fig.3.3, we need only add the line 

  import("music.lib");
at the top to define the noise signal (itself only two lines of Faust code), and change the process definition as follows: 
  process = noise : firpart : + ~ feedback

In summary, the Faust language provides a compact representation for many digital filters, as well as more general digital signal processing, and it is especially useful for quickly generating real-time implementations in various alternative formats. Appendix K gives a number of examples.

Impulse Response

Figure 3.8 plots the impulse response of the example filter, as computed by the matlab script shown in Fig.3.9. (The plot-related commands are also included for completeness.4.2) The impulse signal consists of a single sample at time 0 having amplitude 1, preceded and followed by zeros (an ideal ``click'' at time 0). Linear, time-invariant filters are fully characterized by their response to this simple signal, as we will show in Chapter 4.

Figure 3.8: Impulse response of the example filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$.
\includegraphics[width=4in]{eps/eir}

Figure 3.9: Matlab script for computing and plotting the impulse response of the example digital filter.
  
% Example comb filter:
g1 = 0.5^3;  B = [1 0 0 g1];      % Feedforward coeffs
g2 = 0.9^5;  A = [1 0 0 0 0 g2];  % Feedback coefficients

h = filter(B,A,[1,zeros(1,50)]);  % Impulse response
% h = impz(B,A,50); % alternative in octave-forge or MSPTB

% Matlab-compatible plot:
clf; figure(1); stem([0:50],h,'-k'); axis([0 50 -0.8 1.1]);
ylabel('Amplitude'); xlabel('Time (samples)'); grid;
cmd = 'print -deps ../eps/eir.eps'; disp(cmd); eval(cmd);

The impulse response of this filter is also easy to compute by hand. Referring to the difference equation

$\displaystyle y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5) $

and setting the input $ x(n)$ to $ \delta(n) = [1,0,0,\dots]$, we compute, for $ n=0,1,2,\dots$,

\begin{eqnarray*} y(0) &=& 1\\ y(1) &=& 0\\ y(2) &=& 0\\ y(3) &=& 0.5^3 = 0.1... ...0.9^5)^3 = -0.9^{15} = -0.205891132094649\\ \vdots && \vdots\\ \end{eqnarray*}


Transfer Function

As we cover in Chapter 6, the transfer function of a digital filter is defined as $ H(z)=Y(z)/X(z)$ where $ X(z)$ is the z transform of the input signal $ x(n)$, and $ Y(z)$ is the z transform of the output signal $ y(n)$. We may find $ H(z)$ from Eq.$ \,$(3.1) by taking the z transform of both sides and solving for $ Y(z)/X(z)$:

\begin{eqnarray*} \boldmath {{\cal Z}}\{y(n)\} &=& \boldmath {{\cal Z}}\{x(n) + ... ..._1}}{1 + g_2 z^{-M_2}} \qquad\hbox{(Transfer Function)} \protect \end{eqnarray*}

Some principles of this analysis are as follows:

  1. The z transform $ {\cal Z}\{\cdot\}$ is a linear operator which means, by definition,

    $\displaystyle \boldmath {{\cal Z}}\{\alpha x_1(n) + \beta x_2(n)\} = \alpha \bo... ...\} + \beta \boldmath {{\cal Z}}\{x_2(n)\} \isdef \alpha X_1(z) + \beta X_2(z). $

  2. $ {\cal Z}\{x(n-M)\} = z^{-M}X(z)$. That is, the z transform of a signal $ x(n)$ delayed by $ M$ samples, $ x(n-M),\,n=0,1,2,\ldots$, is $ z^{-M}X(z)$. This is the shift theorem for z transforms, which can be immediately derived from the definition of the z transform, as shown in §6.3.
Note that these two properties of the z transform are all we really need to find the transfer function of any linear, time-invariant digital filter from its difference equation (its implementation formula in the time domain).

In matlab, difference-equation coefficients are specified as transfer-function coefficients (vectors B and A in Fig.3.9). This is why a minus sign is needed in Eq.$ \,$(3.3).

Ok, so finding the transfer function is not too much work. Now, what can we do with it? There are two main avenues of analysis from here: (1) finding the frequency response by setting $ z=e^{j\omega T}$, and (2) factoring the transfer function to find the poles and zeros of the filter. One also uses the transfer function to generate different implementation forms such as cascade or parallel combinations of smaller filters to achieve the same overall filter. The following sections will illustrate these uses of the transfer function on the example filter of this chapter.

Frequency Response

Given the transfer function $ H(z)$, the frequency response is obtained by evaluating it on the unit circle in the complex plane, i.e., by setting $ z=e^{j\omega T}$, where $ T$ is the sampling interval in seconds, and $ \omega$ is radian frequency:4.3

$\displaystyle H(e^{j\omega T}) \isdef \frac{1 + g_1 e^{-jM_1\omega T}}{1 + g_2 ... ... T}}, \qquad -\pi \le \omega T < \pi \qquad\hbox{(Frequency Response)} \protect$ (4.4)

In the special case $ g_1=g_2=1$, we obtain

\begin{eqnarray*} H(e^{j\omega T}) &=& \frac{1 + e^{-jM_1\omega T}}{1 + e^{-jM_2... ...\cos\left(M_1\omega T/2\right)}{\cos\left(M_2\omega T/2\right)}. \end{eqnarray*}

When $ M_1\neq M_2$, the frequency response is a ratio of cosines in $ \omega$ times a linear phase term $ e^{j(M_2-M_1)\omega T/2}$ (which corresponds to a pure delay of $ M_1-M_2$ samples). This special case gives insight into the behavior of the filter as its coefficients $ g_1$ and $ g_2$ approach 1.

When $ M_1=M_2\isdef M$, the filter degenerates to $ H(z)=1$ which corresponds to $ y(n)=x(n)$; in this case, the delayed input and output signals cancel each other out. As a check, let's verify this in the time domain:

\begin{eqnarray*} y(n) &=& x(n) + x(n-M) - y(n-M)\\ &=& x(n) + x(n-M) - [x(n-M... ...) - y(n-3M)]\\ &=& x(n) + y(n-3M)\\ &=& \cdots\\ &=& x(n). \end{eqnarray*}

Amplitude Response

Since the frequency response is a complex-valued function, it has a magnitude and phase angle for each frequency. The magnitude of the frequency response is called the amplitude response (or magnitude frequency response), and it gives the filter gain at each frequency $ \omega$.

In this example, the amplitude response is

$\displaystyle \left\vert H(e^{j\omega T})\right\vert = \left\vert\frac{1 + g_1 e^{-jM_1\omega T}}{1 + g_2 e^{-jM_2\omega T}}\right\vert \protect$ (4.5)

which, for $ g_1=g_2=1$, reduces to

$\displaystyle \left\vert H(e^{j\omega T})\right\vert = \frac{\left\vert\cos\lef... ...a T/2\right)\right\vert}{\left\vert\cos\left(M_2\omega T/2\right)\right\vert}. $

Figure 3.10a shows a graph of the amplitude response of one case of this filter, obtained by plotting Eq.$ \,$(3.5) for $ \omega T \in[-\pi,\pi]$, and using the example settings $ g_1 = 0.5^3$, $ g_2 = 0.9^5$, $ M_1 = 3$, and $ M_2=5$.

Figure 3.10: Frequency response of the example filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$. (a) Amplitude response. (b) Phase response.
\includegraphics[width=\twidth]{eps/efr}

Phase Response

The phase of the frequency response is called the phase response. Like the phase of any complex number, it is given by the arctangent of the imaginary part of $ H(e^{j\omega T})$ divided by its real part, and it specifies the delay of the filter at each frequency. The phase response is a good way to look at short filter delays which are not directly perceivable as causing an ``echo''.4.4 For longer delays in audio, it is usually best to study the filter impulse response, which is output of the filter when its input is $ [1,0,0,0,\ldots]$ (an ``impulse''). We will show later that the impulse response is also given by the inverse z transform of the filter transfer function (or, equivalently, the inverse Fourier transform of the filter frequency response).

In this example, the phase response is

$\displaystyle \angle{H(e^{j\omega T})} = \angle \left(1 + g_1 e^{-jM_1\omega T}\right) - \angle \left ( 1 + g_2 e^{-jM_2\omega T} \right). $

A specific case is plotted in Fig.3.10b, corresponding to the amplitude response in Fig.3.10a. The impulse response is plotted in Fig.3.8. The matlab code for producing these figures is shown in Fig.3.11. (The plotting utility plotfr is given in §J.4.) In Octave or the Matlab Signal Processing Toolbox, a figure similar to Fig.3.10 can be produced by typing simply freqz(B,A,Nspec).

Figure 3.11: Matlab script for computing and displaying the frequency response of the example filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$.
  
% efr.m - frequency response computation in Matlab/Octave

% Example filter:
g1 = 0.5^3; B = [1 0 0 g1];      % Feedforward coeffs
g2 = 0.9^5; A = [1 0 0 0 0 g2];  % Feedback coefficients

Nfft = 1024;         % FFT size
Nspec = 1+Nfft/2;    % Show only positive frequencies
f=[0:Nspec-1]/Nfft;  % Frequency axis
Xnum = fft(B,Nfft);  % Frequency response of FIR part
Xden = fft(A,Nfft);  % Frequency response, feedback part
X = Xnum ./ Xden;    % Should check for divide by zero!

clf; figure(1);      % Matlab-compatible plot
plotfr(X(1:Nspec),f);% Plot frequency response
cmd = 'print -deps ../eps/efr.eps'; disp(cmd); eval(cmd);

Pole-Zero Analysis

Since our example transfer function

$\displaystyle H(z) = \frac{1 + g_1 z^{-M_1}}{1 + g_2 z^{-M_2}} $

(from Eq.$ \,$(3.4)) is a ratio of polynomials in $ z$, and since every polynomial can be characterized by its roots plus a scale factor, we may characterize any transfer function by its numerator roots (called the zeros of the filter), its denominator roots (filter poles), and a constant gain factor:

$\displaystyle H(z) = g \frac{(1-q_1z^{-1})(1-q_2z^{-1})\cdots(1-q_{M_1}z^{-1})}{(1-p_1z^{-1})(1-p_2z^{-1})\cdots(1-p_{M_2}z^{-1})} $

The poles and zeros for this simple example are easy to work out by hand. The zeros are located in the $ z$ plane at

$\displaystyle q_k = - g_1^{\frac{1}{M_1}} e^{j2\pi\frac{k}{M_1}}, \quad k=0,2,\dots,M_1-1 $

where we assume $ g_1>0$, and the poles are similarly given by

$\displaystyle p_k = - g_2^{\frac{1}{M_2}} e^{j2\pi\frac{k}{M_2}}, \quad k=0,2,\dots,M_2-1. $

Figure 3.12 gives the pole-zero diagram of the specific example filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$. There are three zeros, marked by `O' in the figure, and five poles, marked by `X'. Because of the simple form of digital comb filters, the zeros (roots of $ z^3+0.5^3$) are located at 0.5 times the three cube roots of -1 ( $ -e^{2k\pi/3}, k=0,1,2$), and similarly the poles (roots of $ z^5+0.9^5$) are located at 0.9 times the five 5th roots of -1 ( $ -e^{k2\pi/5}, k=0,\dots,4$). (Technically, there are also two more zeros at $ z=0$.) The matlab code for producing this figure is simply

[zeros, poles, gain] = tf2zp(B,A); % Matlab or Octave
zplane(zeros,poles); % Matlab Signal Processing Toolbox
                     % or Octave Forge
where B and A are as given in Fig.3.11. The pole-zero plot utility zplane is contained in the Matlab Signal Processing Toolbox, and in the Octave Forge collection. A similar plot is produced by
sys = tf2sys(B,A,1);
pzmap(sys);
where these functions are both in the Matlab Control Toolbox and in Octave. (Octave includes its own control-systems tool-box functions in the base Octave distribution.)

Figure 3.12: Pole-Zero diagram of the example filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$.
\includegraphics[width=\twidth]{eps/epz}

Alternative Realizations

For actually implementing the example digital filter, we have only seen the difference equation

$\displaystyle y(n) = x(n) + g_1\, x(n-M_1) - g_2\, y(n-M_2) $

(from Eq.$ \,$(3.1), diagrammed in Fig.3.1). While this structure, formally known as ``direct form I'', happens to be a good choice for digital comb filters, there are many other structures to consider in other situations. For example, it is often desirable, for numerical reasons, to implement low-pass, high-pass, and band-pass filters as series second-order sections. On the other hand, digital filters for simulating the vocal tract (for synthesized voice applications) are typically implemented as parallel second-order sections. (When the order is odd, there is one first-order section as well.) The coefficients of the first- and second-order filter sections may be calculated from the poles and zeros of the filter.

We will now illustrate the computation of a parallel second-order realization of our example filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$. As discussed above in §3.11, this filter has five poles and three zeros. We can use the partial fraction expansion (PFE), described in §6.8, to expand the transfer function into a sum of five first-order terms:

\begin{eqnarray*} H(z) &=& \frac{1 + 0.5^3 z^{-3}}{1 + 0.9^5 z^{-5}} \mathrel{=... ...+} \frac{0.4555 + 0.0922z^{-1}}{1 + 0.5562z^{-1}+ 0.8100z^{-2}}, \end{eqnarray*}

where, in the last step, complex-conjugate one-pole sections are combined into real second-order sections. Also, numerical values are given to four decimal places (so `$ =$' is replaced by `$ \approx$' in the second line). In the following subsections, we will plot the impulse responses and frequency responses of the first- and second-order filter sections above.

First-Order Parallel Sections

Figure 3.13 shows the impulse response of the real one-pole section

$\displaystyle H_1(z) = \frac{0.1657}{1 + 0.9z^{-1}}, $

and Fig.3.14 shows its frequency response, computed using the matlab utility myfreqz listed in §7.5.1. (Both Matlab and Octave have compatible utilities freqz, which serve the same purpose.) Note that the sampling rate is set to 1, and the frequency axis goes from 0 Hz all the way to the sampling rate, which is appropriate for complex filters (as we will soon see). Since real filters have Hermitian frequency responses (i.e., an even amplitude response and odd phase response), they may be plotted from 0 Hz to half the sampling rate without loss of information.

Figure 3.13: Impulse response of section 1 of the example filter.
\includegraphics[width=\twidth]{eps/arir1}

Figure 3.14: Frequency response of section 1 of the example filter.
\includegraphics[width=\twidth]{eps/arfr1}

Figure 3.15 shows the impulse response of the complex one-pole section

$\displaystyle H_2(z) = \frac{0.1894 - j 0.0326}{1 - (0.7281 + j 0.5290)z^{-1}}, $

and Fig.3.16 shows the corresponding frequency response.

Figure 3.15: Impulse response of complex one-pole section 2 of the full partial-fraction-expansion of the example filter.
\includegraphics[width=\twidth]{eps/arcir2}

Figure 3.16: Frequency response of complex one-pole section 2.
\includegraphics[width=\twidth]{eps/arcfr2}

The complex-conjugate section,

$\displaystyle H_3(z) = \frac{0.1894 + j 0.0326}{1 - (0.7281 - j 0.5290)z^{-1}}, $

is of course quite similar, and is shown in Figures 3.17 and 3.18.

Figure 3.17: Impulse response of complex one-pole section 3 of the full partial-fraction-expansion of the example filter.
\includegraphics[width=\twidth]{eps/arcir3}

Figure 3.18: Frequency response of complex one-pole section 3.
\includegraphics[width=\twidth]{eps/arcfr3}

Figure 3.19 shows the impulse response of the complex one-pole section

$\displaystyle H_4(z) = \frac{0.2277 + j 0.0202}{1 + (0.2781 + j 0.8560)z^{-1}}, $

and Fig.3.20 shows its frequency response. Its complex-conjugate counterpart, $ H_5(z)$, is not shown since it is analogous to $ H_4(z)$ in relation to $ H_3(z)$.

Figure 3.19: Impulse response of complex one-pole section 4 of the full partial-fraction-expansion of the example filter.
\includegraphics[width=\twidth]{eps/arcir4}

Figure 3.20: Frequency response of complex one-pole section 4.
\includegraphics[width=\twidth]{eps/arcfr4}

Parallel, Real, Second-Order Sections

Figure 3.21 shows the impulse response of the real two-pole section

$\displaystyle H^r_2(z) \isdef H_2(z) + H_3(z) = \frac{0.3788 -0.2413z^{-1}}{1 - 1.4562z^{-1}+ 0.8100z^{-2}}, $

and Fig.3.22 shows its frequency response. The frequency axis unnecessarily extends all the way to the sampling rate ($ f_s=1$).

Figure 3.21: Impulse response of real two-pole section $ H^r_2(z)$ of the real partial-fraction-expansion of the example filter.
\includegraphics[width=\twidth]{eps/arir2}

Figure 3.22: Frequency response of real two-pole section $ H^r_2(z)$.
\includegraphics[width=\twidth]{eps/arfr2}

Finally, Fig.3.23 gives the impulse response of the real two-pole section

$\displaystyle H^r_3(z) \isdef H_4(z) + H_5(z) = \frac{0.4555 + 0.0922z^{-1}}{1 + 0.5562z^{-1}+ 0.8100z^{-2}}, $

and Fig.3.24 its frequency response.

Figure 3.23: Impulse response of real two-pole section $ H^r_3(z)$ of the real partial-fraction-expansion of the example filter.
\includegraphics[width=\twidth]{eps/arir3}

Figure 3.24: Frequency response of real two-pole section $ H^r_3(z)$.
\includegraphics[width=\twidth]{eps/arfr3}

Parallel Second-Order Signal Flow Graph

Figure 3.25 shows the signal flow graph for the implementation of our example filter using parallel second-order sections (with one first-order section since the number of poles is odd). This is the same filter as that shown in Fig.3.1 with $ g_1 = 0.5^3$, $ g_2 = 0.9^5$, $ M_1 = 3$, and $ M_2=5$. The second-order sections are special cases of the ``biquad'' filter section, which is often implemented in software (and chip) libraries. Any digital filter can be implemented as a sum of parallel biquads by finding its transfer function and computing the partial fraction expansion.

Figure 3.25: Signal flow graph for the re-implementation of the example filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$ as a parallel bank of real first- and second-order digital filter sections.
\begin{figure}\input fig/epfe.pstex_t \end{figure}

Figure: Same as Fig.3.25 except using transposed direct-form-II biquad sections in place of direct-form-II biquad sections.
\begin{figure}\input fig/epfet.pstex_t \end{figure}

The two second-order biquad sections in Fig.3.25 are in so-called ``Direct-Form II'' (DF-II) form. In Chapter 9, a total of four direct-form filter implementations will be discussed, along with some other commonly used implementation structures. In particular, it is explained there why Transposed Direct-Form II (TDF-II) is usually a better choice of implementation structure for IIR filters when numerical dynamic range is limited (as it is in fixed-point ``DSP chips''). Figure 3.26 shows how our example looks using TDF-II biquads in place of the DF-II biquads of Fig.3.25.

Series, Real, Second-Order Sections

Converting the difference equation $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$ to a series bank of real first- and second-order sections is comparatively easy. In this case, we do not need a full blown partial fraction expansion. Instead, we need only factor the numerator and denominator of the transfer function into first- and/or second-order terms. Since a second-order section can accommodate up to two poles and two zeros, we must decide how to group pairs of poles with pairs of zeros. Furthermore, since the series sections can be implemented in any order, we must choose the section ordering. Both of these choices are normally driven in practice by numerical considerations. In fixed-point implementations, the poles and zeros are grouped such that dynamic range requirements are minimized. Similarly, the section order is chosen so that the intermediate signals are well scaled. For example, internal overflow is more likely if all of the large-gain sections appear before the low-gain sections. On the other hand, the signal-to-quantization-noise ratio will deteriorate if all of the low-gain sections are placed before the higher-gain sections. For further reading on numerical considerations for digital filter sections, see, e.g., [103].

Summary

This chapter has provided an overview of practical digital filter implementation, analysis, transformation, and display. The principal time-domain views were the difference equation and the impulse response. The signal flow graph applies to either time- or frequency-domain signals. The frequency-domain analyses were derived from the z transform of the difference equation, yielding the transfer function, frequency response, amplitude response, phase response, pole-zero analysis, alternate realizations, and related topics. We will take up these topics and more in the following chapters.

In the next chapter, we begin a more systematic presentation of the main concepts of linear systems theory. After that, we will return to practical digital filter analysis, implementation, and (elementary) design.

Next Section:
Linear Time-Invariant Digital Filters
Previous Section:
Matlab Analysis of the Simplest Lowpass Filter