Implementation Structures for Recursive Digital Filters
This chapter introduces the four directform filter implementations, and discusses implementation of filters as parallel or series combinations of smaller filter sections. A careful study of filter forms can be important when numerical issues arise, such as when implementing a digital filter in a fixedpoint processor. (The least expensive ``DSP chips'' use fixedpoint numerical processing, sometimes with only a 16bit word length in lowcost audio applications.) For implementations in floatingpoint arithmetic, especially at 32bit wordlengths or greater, the choice of filter implementation structure is usually not critical. In matlab software, for example, one rarely uses anything but the filter function, which is implemented in doubleprecision floating point (typically 64 bits at the time of this writing, and more bits internally in the floatingpoint unit).The Four Direct Forms
DirectForm I
As mentioned in §5.5, the difference equation(10.1) 
specifies the DirectForm I (DFI) implementation of a digital filter [60]. The DFI signal flow graph for the secondorder case is shown in Fig.9.1. The DFI structure has the following properties:
 It can be regarded as a twozero filter section followed in series by a twopole filter section.
 In most fixedpoint arithmetic schemes (such as two's complement, the most commonly used [84]^{10.1}) there is no possibility of internal filter overflow. That is, since there is fundamentally only one summation point in the filter, and since fixedpoint overflow naturally ``wraps around'' from the largest positive to the largest negative number and vice versa, then as long as the final result is ``in range'', overflow is avoided, even when there is overflow of intermediate results in the sum (see below for an example). This is an important, valuable, and unusual property of the DFI filter structure.
 There are twice as many delays as are necessary. As a result, the DFI structure is not canonical with respect to delay. In general, it is always possible to implement an thorder filter using only delay elements.
 As is the case with all directform filter structures (those which have coefficients given by the transferfunction coefficients), the filter poles and zeros can be very sensitive to roundoff errors in the filter coefficients. This is usually not a problem for a simple secondorder section, such as in Fig.9.1, but it can become a problem for higher order directform filters. This is the same numerical sensitivity that polynomial roots have with respect to polynomialcoefficient roundoff. As is well known, the sensitivity tends to be larger when the roots are clustered closely together, as opposed to being well spread out in the complex plane [18, p. 246]. To minimize this sensitivity, it is common to factor filter transfer functions into series and/or parallel secondorder sections, as discussed in §9.2 below.
Two's Complement WrapAround
In this section, we give an example showing how temporary overflow in two's complement fixedpoint causes no ill effects. In 3bit signed fixedpoint arithmetic, the available numbers are as shown in Table 9.1.

Let's perform the sum , which gives a temporary overflow (, which wraps around to ), but a final result () which is in the allowed range :^{10.3}
Direct Form II
The signal flow graph for the DirectFormII (DFII) realization of the secondorder IIR filter section is shown in Fig.9.2. The difference equation for the secondorder DFII structure can be written as It can be regarded as a twopole filter section followed by a twozero filter section.
 It is canonical with respect to delay. This happens because delay elements associated with the twopole and twozero sections are shared.
 In fixedpoint arithmetic, overflow can occur at the delayline input (output of the leftmost summer in Fig.9.2), unlike in the DFI implementation.
 As with all directform filter structures, the poles and zeros are sensitive to roundoff errors in the coefficients and , especially for high transferfunction orders. Lower sensitivity is obtained using series loworder sections (e.g., second order), or by using ladder or lattice filter structures [86].
More about Potential Internal Overflow of DFII
Since the poles come first in the DFII realization of an IIR filter, the signal entering the state delayline (see Fig.9.2) typically requires a larger dynamic range than the output signal . In other words, it is common for the feedback portion of a DFII IIR filter to provide a large signal boost which is then compensated by attenuation in the feedforward portion (the zeros). As a result, if the input dynamic range is to remain unrestricted, the two delay elements may need to be implemented with highorder guard bits to accommodate an extended dynamic range. If the number of bits in the delay elements is doubled (which still does not guarantee impossibility of internal overflow), the benefit of halving the number of delays relative to the DFI structure is approximately canceled. In other words, the DFII structure, which is canonical with respect to delay, may require just as much or more memory as the DFI structure, even though the DFI uses twice as many addressable delay elements for the filter state memory.Transposed DirectForms
The remaining two direct forms are obtained by formally transposing directforms I and II [60, p. 155]. Filter transposition may also be called flow graph reversal, and transposing a SingleInput, SingleOutput (SISO) filter does not alter its transfer function. This fact can be derived as a consequence of Mason's gain formula for signal flow graphs [49,50] or Tellegen's theorem (which implies that an LTI signal flow graph is interreciprocal with its transpose) [60, pp. 176177]. Transposition of filters in statespace form is discussed in §G.5. The transpose of a SISO digital filter is quite straightforward to find: Reverse the direction of all signal paths, and make obviously necessary accommodations. ``Obviously necessary accommodations'' include changing signal branchpoints to summers, and summers to branchpoints. Also, after this operation, the input signal, normally drawn on the left of the signal flow graph, will be on the right, and the output on the left. To renormalize the layout, the whole diagram is usually leftright flipped. Figure 9.3 shows the TransposedDirectFormI (TDFI) structure for the general secondorder IIR digital filter, and Fig.9.4 shows the TransposedDirectFormII (TDFII) structure. To facilitate comparison of the transposed with the original, the inputs and output signals remain ``switched'', so that signals generally flow righttoleft instead of the usual lefttoright. (Exercise: Derive forms TDFI/II by transposing the DFI/II structures shown in Figures 9.1 and 9.2.)Numerical Robustness of TDFII
An advantage of the transposed directform II structure (depicted in Fig.9.4) is that the zeros effectively precede the poles in series order. As mentioned above, in many digital filters design, the poles by themselves give a large gain at some frequencies, and the zeros often provide compensating attenuation. This is especially true of filters with sharp transitions in their frequency response, such as the ellipticfunctionfilter example on page ; in such filters, the sharp transitions are achieved using near polezero cancellations close to the unit circle in the plane.^{10.4}Series and Parallel Filter Sections
In many situations it is better to implement a digital filter in terms of first and/or secondorder elementary sections, either in series or in parallel. In particular, such an implementation may have numerical advantages. In the timevarying case, it is easier to control fundamental parameters of small sections, such as pole frequencies, by means of coefficient variations.Series SecondOrder Sections
For many filter types, such as lowpass, highpass, and bandpass filters, a good choice of implementation structure is often series secondorder sections. In fixedpoint applications, the ordering of the sections can be important. The matlab function tf2sos^{10.5} converts from ``transfer function form'', , to series ``secondordersection'' form. For example, the lineBAMatrix = tf2sos(B,A);converts the real filter specified by polynomial vectors B and A to a series of secondorder sections (biquads) specified by the rows of BAMatrix. Each row of BAMatrix is of the form . The function tf2sos may be implemented simply as a call to tf2zp followed by a call to zp2sos, where the zp form of a digital filter consists of its (possibly complex) zeros, poles, and an overall gain factor:
function [sos,g] = tf2sos(B,A) [z,p,g]=tf2zp(B(:)',A(:)'); % Direct form to (zeros,poles,gain) sos=zp2sos(z,p,g); % (z,p,g) to series secondorder sections
Matlab Example
The following matlab example expands the filterB=[1 0 0 0 0 1]; A=[1 0 0 0 0 .9]; [sos,g] = tf2sos(B,A) sos = 1.00000 0.61803 1.00000 1.00000 0.60515 0.95873 1.00000 1.61803 1.00000 1.00000 1.58430 0.95873 1.00000 1.00000 0.00000 1.00000 0.97915 0.00000 g = 1The g parameter is an input (or output) scale factor; for this filter, it was not needed. Thus, in this example we obtained the following filter factorization:
% Numerically challenging "clustered roots" example: [B,A] = zp2tf(ones(10,1),0.9*ones(10,1),1); [sos,g] = tf2sos(B,A); [Bh,Ah] = sos2tf(sos,g); format long; disp(sprintf('Relative L2 numerator error: %g',... norm(BhB)/norm(B))); % Relative L2 numerator error: 1.26558e15 disp(sprintf('Relative L2 denominator error: %g',... norm(AhA)/norm(A))); % Relative L2 denominator error: 1.65594e15Thus, in this test, the original directform filter is compared with one created from the secondorder sections. Such checking should be done for highorder filters, or filters having many poles and/or zeros close together, because the polynomial factorization used to find the poles and zeros can fail numerically. Moreover, the stability of the factors should be checked individually.
Parallel First and/or SecondOrder Sections
Instead of breaking up a filter into a series of secondorder sections, as discussed in the previous section, we can break the filter up into a parallel sum of first and/or secondorder sections. Parallel sections are based directly on the partial fraction expansion (PFE) of the filter transfer function discussed in §6.8. As discussed in §6.8.3, there is additionally an FIR part when the order of the transferfunction denominator does not exceed that of the numerator (i.e., when the transfer function is not strictly proper). The most general case of a PFE, valid for any finiteorder transfer function, was given by Eq.(6.19), repeated here for convenience:where denotes the number of distinct poles, and denotes the multiplicity of the th pole. The polynomial is the transfer function of the FIR part, as discussed in §6.8.3. The FIR part is typically realized as a tapped delay line, as shown in Fig.5.5.
FirstOrder Complex Resonators
For distinct poles, the recursive terms in the complete partial fraction expansion of Eq.(9.2) can be realized as a parallel sum of complex onepole filter sections, thereby producing a parallel complex resonator filter bank. Complex resonators are efficient for processing complex input signals, and they are especially easy to work with. Note that a complex resonator bank is similarly obtained by implementing a diagonalized statespace model [Eq.(G.22)].Real SecondOrder Sections
In practice, however, signals are typically realvalued functions of time. As a result, for real filters (§5.1), it is typically more efficient computationally to combine complexconjugate onepole sections together to form real secondorder sections (two poles and one zero each, in general). This process was discussed in §6.8.1, and the resulting transfer function of each secondorder section becomeswhere is one of the poles, and is its corresponding residue. This is a special case of the biquad section discussed in §B.1.6. When the two poles of a real secondorder section are complex, they form a complexconjugate pair, i.e., they are located at in the plane, where is the modulus of either pole, and is the angle of either pole. In this case, the ``resonancetuning coefficient'' in Eq.(9.3) can be expressed as
Implementation of Repeated Poles
Fig.9.5 illustrates an efficient implementation of terms due to a repeated pole with multiplicity three, contributing the additive termsFormant Filtering Example
In speech synthesis [27,39], digital filters are often used to simulate formant filtering by the vocal tract. It is well known [23] that the different vowel sounds of speech can be simulated by passing a ``buzz source'' through a only two or three formant filters. As a result, speech is fully intelligible through the telephone bandwidth (nominally only 2003200 Hz). A formant is a resonance in the voice spectrum. A single formant may thus be modeled using one biquad (secondorder filter section). For example, in the vowel as in ``father,'' the first three formant centerfrequencies have been measured near 700, 1220, and 2600 Hz, with halfpower bandwidths^{10.7} 130, 70, and 160 Hz [40]. In principle, the formant filter sections are in series, as can be found by deriving the transfer function of an acoustic tube [48]. As a consequence, the vocaltract transfer function is an allpole filter (provided that the nasal tract is closed off or negligible). As a result, there is no need to specify gains for the formant resonatorsonly centerfrequency and bandwidth are necessary to specify each formant, leaving only an overall scale factor unspecified in a cascade (series) formant filter bank. Numerically, however, it makes more sense to implement disjoint resonances in parallel rather than in series.^{10.8} This is because when one formant filter is resonating, the others will be attenuating, so that to achieve a particular peakgain at resonance, the resonating filter must overcome all combined attenuations as well as applying its own gain. In fixedpoint arithmetic, this can result in large quantizationnoise gains, especially for the last resonator in the chain. As a result of these considerations, our example will implement the formant sections in parallel. This means we must find the appropriate biquad numerators so that when added together, the overall transferfunction numerator is a constant. This will be accomplished using the partial fraction expansion (§6.8).^{10.9} The matlab below illustrates the construction of a parallel formant filter bank for simulating the vowel . For completeness, it is used to filter a bandlimited impulse train, in order to synthesize the vowel sound.F = [700, 1220, 2600]; % Formant frequencies (Hz) BW = [130, 70, 160]; % Formant bandwidths (Hz) fs = 8192; % Sampling rate (Hz) nsecs = length(F); R = exp(pi*BW/fs); % Pole radii theta = 2*pi*F/fs; % Pole angles poles = R .* exp(j*theta); % Complex poles B = 1; A = real(poly([poles,conj(poles)])); % freqz(B,A); % View frequency response: % Convert to parallel complex onepoles (PFE): [r,p,f] = residuez(B,A); As = zeros(nsecs,3); Bs = zeros(nsecs,3); % complexconjugate pairs are adjacent in r and p: for i=1:2:2*nsecs k = 1+(i1)/2; Bs(k,:) = [r(i)+r(i+1), (r(i)*p(i+1)+r(i+1)*p(i)), 0]; As(k,:) = [1, (p(i)+p(i+1)), p(i)*p(i+1)]; end sos = [Bs,As]; % standard secondordersection form iperr = norm(imag(sos))/norm(sos); % make sure sos is ~real disp(sprintf('imag(sos)/sos = %g',iperr)); % 1.6e16 sos = real(sos) % and make it exactly real % Reconstruct original numerator and denominator as a check: [Bh,Ah] = psos2tf(sos); % parallel sos to transfer function % psos2tf appears in the matlabutilities appendix disp(sprintf('AAh = %g',norm(AAh))); % 5.77423e15 % Bh has trailing epsilons, so we'll zeropad B: disp(sprintf('BBh = %g',... norm([B,zeros(1,length(Bh)length(B))]  Bh))); % 1.25116e15 % Plot overlay and sum of all three % resonator amplitude responses: nfft=512; H = zeros(nsecs+1,nfft); for i=1:nsecs [Hiw,w] = freqz(Bs(i,:),As(i,:)); H(1+i,:) = Hiw(:).'; end H(1,:) = sum(H(2:nsecs+1,:)); ttl = 'Amplitude Response'; xlab = 'Frequency (Hz)'; ylab = 'Magnitude (dB)'; sym = ''; lgnd = {'sum','sec 1','sec 2', 'sec 3'}; np=nfft/2; % Only plot for positive frequencies wp = w(1:np); Hp=H(:,1:np); figure(1); clf; myplot(wp,20*log10(abs(Hp)),sym,ttl,xlab,ylab,1,lgnd); disp('PAUSING'); pause; saveplot('../eps/lpcexovl.eps'); % Now synthesize the vowel [a]: nsamps = 256; f0 = 200; % Pitch in Hz w0T = 2*pi*f0/fs; % radians per sample nharm = floor((fs/2)/f0); % number of harmonics sig = zeros(1,nsamps); n = 0:(nsamps1); % Synthesize bandlimited impulse train for i=1:nharm, sig = sig + cos(i*w0T*n); end; sig = sig/max(sig); speech = filter(1,A,sig); soundsc([sig,speech]); % hear buzz, then 'ah'Notes:
 The sampling rate was chosen to be Hz because that is the default Matlab sampling rate, and because that is a typical value used for ``telephone quality'' speech synthesis.
 The psos2tf utility is listed in §J.7.
 The overlay of the amplitude responses are shown in Fig.9.6.
Butterworth Lowpass Filter Example
This example illustrates the design of a 5thorder Butterworth lowpass filter, implementing it using secondorder sections. Since all three sections contribute to the same passband and stopband, it is numerically advisable to choose a series secondordersection implementation, so that their passbands and stopbands will multiply together instead of add.fc = 1000; % Cutoff frequency (Hz) fs = 8192; % Sampling rate (Hz) order = 5; % Filter order [B,A] = butter(order,2*fc/fs); % [0:pi] maps to [0:1] here [sos,g] = tf2sos(B,A) % sos = % 1.00000 2.00080 1.00080 1.00000 0.92223 0.28087 % 1.00000 1.99791 0.99791 1.00000 1.18573 0.64684 % 1.00000 1.00129 0.00000 1.00000 0.42504 0.00000 % % g = 0.0029714 % % Compute and display the amplitude response Bs = sos(:,1:3); % Section numerator polynomials As = sos(:,4:6); % Section denominator polynomials [nsec,temp] = size(sos); nsamps = 256; % Number of impulseresponse samples % Note use of input scalefactor g here: x = g*[1,zeros(1,nsamps1)]; % SCALED impulse signal for i=1:nsec x = filter(Bs(i,:),As(i,:),x); % Series sections end % %plot(x); % Plot impulse response to make sure % it has decayed to zero (numerically) % % Plot amplitude response % (in Octave  Matlab slightly different): figure(2); X=fft(x); % sampled frequency response f = [0:nsamps1]*fs/nsamps; grid('on'); axis([0 fs/2 100 5]); legend('off'); plot(f(1:nsamps/2),20*log10(X(1:nsamps/2)));The final plot appears in Fig.9.7. A Matlab function for frequency response plots is given in §J.4. (Of course, one can also use freqz in either Matlab or Octave, but that function uses subplots which are not easily printable in Octave.) Note that the Matlab Signal Processing Toolbox has a function called sosfilt so that ``y=sosfilt(sos,x)'' will implement an array of secondorder sections without having to unpack them first as in the example above.
Summary of Series/Parallel Filter Sections
In summary, we noted above the following general guidelines regarding series vs. parallel elementarysection implementations: Series sections are preferred when all sections contribute to the same passband, such as in a lowpass, highpass, bandpass, or bandstop filter.
 Parallel sections are usually preferred when the sections have disjoint passbands, such as a formant filter bank used in voice models. Another example would be the phase vocoder filter bank [21].
PoleZero Analysis Problems
See http://ccrma.stanford.edu/~jos/filtersp/Pole_Zero_Analysis_Problems.html.Next Section:
Filters Preserving Phase
Previous Section:
PoleZero Analysis