Chapter 10

Analysis and Design of Filters

Chapter Objectives

• Develop the concept of analog and discrete-time filters.
• Learn ideal filter behavior and the conditions for distortionless transmission.
• Learn design specifications for realizable filters.
• Explore methods for designing realizable analog filters using Butterworth and Chebyshev approximation formulas.
• Study design methods for obtaining IIR discrete-time filters from analog prototype filters.
• Study design methods for obtaining FIR discrete-time filters.

10.1 Introduction

In many signal processing applications the need arises to change the strength, or the relative significance, of various frequency components in a given signal. Sometimes we may need to eliminate certain frequency components in a signal; at other times we may need to boost the strength of a range of frequencies over others. This act of changing the relative amplitudes of frequency components in a signal is referred to as filtering, and the system that facilitates this is referred to as a filter.

In Chapters 4 and 7 we have developed techniques for transform-domain analysis of CTLTI systems. Specifically, the relationship between the input and the output signals of a CTLTI system is expressed in the transform domain as

$$\begin{array}{rr}\hfill Y\left(\omega \right)=H\left(\omega \right)X\left(\omega \right)& \hfill \left(10.1\right)\end{array}$$

in terms of the Fourier transforms of the signals involved, or as

$$\begin{array}{ll}Y\left(s\right)=H\left(s\right)X\left(s\right)\hfill & \left(10.2\right)\hfill \end{array}$$

using the Laplace transforms of the signals. In either case, the system function H (ω) or H (s) serves as a multiplier function that shapes the spectrum of the input signal to create an output signal.

In a general sense any CTLTI system can be seen as a filter. Consider, for example, a communication channel the purpose of which is to transmit an information bearing signal from one point to another. Assuming that the particular channel we are considering can be reasonably modeled as a CTLTI system, we have the model shown in Fig. 10.1.

Ideally we would like the signal at the output of the channel, $\overline{x}\left(t\right)$, to be identical to the signal transmitted, x (t). This is generally not the case, however, and the signal $\overline{x}\left(t\right)$ is a distorted version of x (t). Therefore the channel acts as a filter the action of which is an undesired one. In order to compensate for the distortion effect of the channel, a second system may be connected in cascade with the channel as shown in Fig. 10.2. This second system is referred to as an equalizer. Its purpose is to reduce the distortion effect caused by the channel (or filter) so that the output y (t) is a better representation of the original signal x (t).

If we want the output y (t) to be identical to the original input x (t), it follows from Fig. 10.2 that the cascade combination of the channel and the equalizer must form an identity system, that is,

$$\begin{array}{ll}{H}_e\left(s\right)H_e\left(s\right)=1\hfill & \left(10.3\right)\hfill \end{array}$$

This requires the system function of the equalizer to be

$$\begin{array}{ll}{H}_e\left(s\right)=\frac{1}{H_e\left(s\right)}\hfill & \left(10.4\right)\hfill \end{array}$$

In this case Hc (s) represents a filter (channel) with undesired effects on the signal transmitted through it. In contrast, He (s) represents a filter that we may design for the purpose of eliminating, or at least reducing, the undesired effects of the channel.

10.2 Distortionless Transmission

A system driven by an input signal x (t) produces an output signal y (t) that is, in general, different in shape from the input signal. Thus the system modifies the signal as it transmits it from its input to its output. This modification may be a desired effect of the system or an undesired side effect. Regardless, it is referred to as distortion.

Consider a CTLTI system with system function H (s) driven by an input signal x (t). The input and the output signals are related to each other through the convolution relationship

$$\begin{array}{ll}y\left(t\right)=h\left(t\right)*x\left(t\right)={\int }_{-\infty }^{\infty }h\left(\lambda \right)x\left(t-\lambda \right)d\lambda \hfill & \left(10.5\right)\hfill \end{array}$$

The frequency spectra of input and output signals are related to each other by

$$\begin{array}{ll}Y\left(\omega \right)=H\left(\omega \right)X\left(\omega \right)\hfill & \left(10.6\right)\hfill \end{array}$$

The type of distortion caused by a CTLTI system is linear distortion , and it exhibits itself in two distinct forms, namely amplitude distortion and phase distortion. Let us write the magnitude spectrum of the output signal from Eqn. (10.6):

$$\begin{array}{ll}\left|Y\left(\omega \right)\right|=\left|H\left(\omega \right)\right|\left|X\left(\omega \right)\right|\hfill & \left(10.7\right)\hfill \end{array}$$

Considering that the input signal x (t) is a specific mixture of sinusoidal components at various frequencies, Eqn. (10.7) describes how the system modifies the amplitude of each of these sinusoidal components. The scale factor applied by the system to a sinusoid at a particular frequency may be different from the scale factor applied to a sinusoid at another frequency. This non-uniform scaling of sinusoidal components of the input signal leads to amplitude distortion.

The phase spectrum of the output signal is found from Eqn. (10.6) as

$$\begin{array}{ll}\measuredangle Y\left(\omega \right)=\measuredangle H\omega +\measuredangle X\left(\omega \right)\hfill & \left(10.8\right)\hfill \end{array}$$

The system adds an offset to the phase of each sinusoidal component of the input signal. Phase-shifting a sinusoidal component of the input signal causes a corresponding time shift (delay). If signal components are delayed by differing amounts, their relative alignment in time is disturbed, causing the output signal to be different from the input signal. This effect is known as phase distortion.

What type of a system function would be needed in order to obtain a distortionless version of the input signal at the output of the system? As a starting point for discussion let us assume that what we mean by distortionless transmission is the output signal being identical to the input signal, that is,

$$\begin{array}{ll}y\left(t\right)=x\left(t\right)\hfill & \left(10.9\right)\hfill \end{array}$$

This clearly requires that we have

$$\begin{array}{ll}Y\left(\omega \right)=X\left(\omega \right)\hfill & \left(10.10\right)\hfill \end{array}$$

which in turn requires

$$\begin{array}{lll}H\left(\omega \right)=1,\hfill & \text{all }\omega \hfill & \left(10.11\right)\hfill \end{array}$$

and

$$\begin{array}{ll}h\left(t\right)=\delta \left(t\right)\hfill & \left(10.12\right)\hfill \end{array}$$

Thus, we need the impulse response of the system to be an impulse. While Eqns. (10.11) and (10.12) represent a mathematically satisfying result, the corresponding system is physically unrealizable. No physical system is capable of responding to an input signal instantaneously, without a time delay between the input and the output. In order to come up with a realizable result, we will relax our definition of distortionless transmission, and accept an output signal in the form

$$\begin{array}{ll}y\left(t\right)=Kx\left(t-t_d\right)\hfill & \left(10.13\right)\hfill \end{array}$$

where K is a constant scale factor and td is a constant time-delay. Thus, y (t) in Eqn. (10.13) is a scaled and delayed version of the input signal x (t). Taking the Fourier transform of both sides of Eqn. (10.13) yields

$$\begin{array}{ll}Y\left(\omega \right)=K{e}^{-j\omega t_d}X\left(\omega \right)\hfill & \left(10.14\right)\hfill \end{array}$$

and the system function needed for distortionless transmission is

$$\begin{array}{ll}H\left(\omega \right)=K{e}^{-j\omega t_d}\hfill & \left(10.15\right)\hfill \end{array}$$

Magnitude and phase characteristics of the system function are

$$\begin{array}{ll}\left|H\left(\omega \right)\right|=K\hfill & \left(10.16\right)\hfill \end{array}$$

and

$$\begin{array}{cc}\Theta \left(\omega \right)=\measuredangle H\left(\omega \right)=-\omega t_d& \left(10.17\right)\end{array}$$

A distortion-free system has a magnitude that is constant, and a phase that is a linear function of frequency as shown in Fig. 10.3.

Example 10.1: Distortion caused by the RC circuit

Consider again the simple RC circuit that was used in several examples in Chapters 2, 4 and 7.

In some cases the RC circuit can be used as a simple model for a wired electrical connection between two points. Such a model would provide a reasonably good approximation to the characteristics of the wired connection provided that the distance between the two points being connected is not too long. (For long distance connections the resistance and the capacitance would have to be modeled as distributed parameters rather than being concentrated at one point each.) The system function for the RC circuit is

$$H\left(\omega \right)=\frac{1}{1+j\omega RC}$$

Clearly this is not a distortionless system as defined by Eqn. (10.15), however, under certain circumstances it can be considered near distortion-free. The magnitude of the system function is

$$\begin{array}{cc}\left|H\left(\omega \right)\right|=\frac{1}{\sqrt{1+{\left(\omega /{\omega }_c\right)}^2}}& \left(10.18\right)\end{array}$$

with the 3-dB cutoff frequency defined as

$$\begin{array}{cc}{\omega }_c\triangleq \frac{1}{RC}& \left(10.19\right)\end{array}$$

This magnitude characteristic is shown in Fig. 10.5(a). The phase of the system function is

$$\begin{array}{cc}\Theta \left(\omega \right)=\measuredangle H\left(\omega \right)=-{tan}^{-1}\left(\frac{\omega }{{\omega }_c}\right)& \left(10.20\right)\end{array}$$

and the corresponding time delay is

$$\begin{array}{cc}{t}_d\left(\omega \right)=-\frac{\Theta \left(\omega \right)}{\omega }=\frac{1}{\omega }{tan}^{-1}\left(\frac{\omega }{{\omega }_c}\right)& \left(10.21\right)\end{array}$$

The phase characteristic is shown in Fig. 10.5(b).

The magnitude characteristic of the system seems to be relatively flat for frequencies in the vicinity of ω = 0. Similarly, the phase characteristic is close to a straight line, and the time delay is relatively constant for low frequencies. Thus, if we limit the use of the system to input signals that contain only the frequencies ω < < ωc, we could consider the system as near distortion-free.

What is the frequency range (in terms of ωc) in which one can consider the RC circuit approximately distortionless? Of course, the answer depends on the amount of variation in magnitude and/or time delay that we consider negligible. Since both the magnitude and the time-delay characteristics exhibit their peak values at ω = 0, and both are symmetric around that point, we could simply look for a frequency at which the characteristic drops down by a certain percentage. Let ω0 be the frequency at which the magnitude of the system function drops below the peak value by p percent. Using Eqn. (10.18) and realizing that the maximum magnitude is H (0) = 1, we can write

$$\begin{array}{cc}\frac{1}{\sqrt{1+{\left({\omega }_0/{\omega }_c\right)}^2}}\le \frac{100-p}{100}& \left(10.22\right)\end{array}$$

Solving this inequality for ω0 yields

$$\begin{array}{cc}{\omega }_0\le {\omega }_c\sqrt{{\left(\frac{100}{100-p}\right)}^2-1}& \left(10.23\right)\end{array}$$

For the frequencies |ω| ≤ 0.1425 ωc, the variation in the magnitude of the system function will be 1 percent or less. If we are willing to tolerate up to 5 percent magnitude variation, then the usable frequency range is |ω| ≤ 0.3287 ωc.

Next consider the phase characteristic given by Eqn. (10.20). Differentiating it with respect to frequency, we obtain

$$\begin{array}{cc}\frac{d\Theta }{d\omega }=-\frac{1}{{\omega }_c}\left[\frac{1}{1+{\left(\omega /{\omega }_c\right)}^2}\right]& \left(10.24\right)\end{array}$$

and the slope of the phase characteristic at frequency ω = 0 is

$$\begin{array}{cc}{\frac{d\Theta }{d\omega }|}_{\omega =0}=-\frac{1}{{\omega }_c}& \left(10.25\right)\end{array}$$

The line that is tangent to the phase characteristic at the frequency ω = 0 is also shown in Fig. 10.5(b) for comparison. The phase remains fairly close to the tangent line for low frequencies. Similar analysis can be performed for the time-delay characteristic of the system; however, because of the form of Eqn. (10.21), arriving at an analytical result is tedious. Instead of trying to duplicate the methodology used in Eqns. (10.22) and (10.23), let us just take the frequency ranges obtained for specified magnitude variations above, and determine the corresponding variation in time delay. It can easily be shown that

$$\begin{array}{cc}{t}_d\left(0\right)=0.1592/f_c& \left(10.26\right)\end{array}$$

$$\begin{array}{cc}{t}_d\left(0.1425f_c\right)=0.1581/f_c& \left(10.27\right)\end{array}$$

$$\begin{array}{cc}{t}_d\left(0.3287f_c\right)=0.1538/f_c& \left(10.28\right)\end{array}$$

If we accept 1 percent variation in magnitude, the corresponding variation in time delay is 0.67 percent. For a 5 percent magnitude variation, the time delay varies by 3.38 percent.

Software resources:

ex_10_1.m

10.3 Ideal Filters

The conditions for a CTLTI system to be distortionless were identified in Eqns. (10.16) and (10.17) of the previous section, and were graphically illustrated in Fig. 10.3. We observe that, in order to be truly distortionless, a system must have infinite bandwidth. In other words, it must be able to respond to signals containing an infinite range of frequencies. No physically realizable system is capable of doing this.

Furthermore, signals encountered in engineering applications have limited bandwidth, that is, they contain a limited range of frequencies. Consider such a signal x (t). A CTLTI system that satisfies Eqns. (10.16) and (10.17) not necessarily at every frequency, but at all frequencies contained in the signal x (t), would still facilitate distortionless transmission of this signal. Ideal frequency-selective filters are systems that allow distortionless transmission of signal components in a particular frequency range while eliminating all other signal components.

One type is an ideal lowpass filter that is characterized by the system function

$$\begin{array}{ccc}{H}_{LP}\left(\omega \right)\hfill & =\prod \left(\frac{\omega }{2{\omega }_0}\right)e^{-j\omega t_d}\hfill & \hfill \\ \hfill & =\{\begin{array}{ll}{e}^{-j\omega t_d},\hfill & \left|\omega \right|\le {\omega }_0\hfill \\ 0,\hfill & \left|\omega \right|>{\omega }_0\hfill \end{array}\hfill & \left(10.29\right)\hfill \end{array}$$

Any signal that contains only the frequencies in the range |ω| ≤ ω0 is transmitted through this system without any distortion of its magnitude or phase characteristics. The magnitude and the phase of the system function HLP (ω) of an ideal lowpass filter are graphed in Fig. 10.6.

The magnitude of the system function is constant for frequencies in the range −ω0 ≤ ω ≤ ω0. The phase of the system function is a straight line that passes through the origin, with the slope −ωtd. Signals that do not contain any frequencies outside the range |ω| ≤ ω0 would be transmitted through this system without experiencing any distortion. The frequency ω0 in Eqn. (10.29) is the cutoff frequency of the ideal lowpass filter. The parameter td indicates the time delay caused by the filter.

Some signals are bandpass in nature, that is, they contain frequencies in a bandpass frequency range ω1 ≤ |ω| < ω2. Examples of these signals are encountered often in the study of communication systems. The system function of an ideal bandpass filter that would allow distortionless transmission of the appropriate bandpass signal is given by

$$\begin{array}{ll}{H}_{BP}\left(\omega \right)=\{\begin{array}{ll}{e}^{-j\omega t_d},\hfill & {\omega }_1\le \left|\omega \right|\le {\omega }_2\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\hfill & \left(10.30\right)\hfill \end{array}$$

and is illustrated in Fig. 10.7.

Similarly, ideal highpass and band-reject filters can also be defined. An ideal highpass filter transmits, without distortion, signals that contain only frequencies greater than a minimum value, that is, |ω| ≥ ω0. An ideal band-reject filter, on the other hand, facilitates distortionless transmission of signals that contain only those frequencies that are outside a specified range. Frequencies that satisfy either |ω| ≤ ω1 or |ω| ≥ ω2 are transmitted; all other frequencies are suppressed. Figs. 10.8 and 10.9 illustrate the system functions of ideal highpass and band-reject filters respectively.

The ideal filters discussed above are not practically realizable. In order to see why this is the case, let us determine the impulse response of an ideal filter. For the ideal lowpass filter of Eqn. (10.29), the impulse response can be easily computed from the system function through inverse Fourier transform. Using the inverse Fourier transform relationship given by Eqn. (4.126) on the ideal lowpass filter system function in Eqn. (10.29) we obtain

$$\begin{array}{ll}{h}_{LP}\left(t\right)={ℱ}^{-1}\left\{H_{LP}\left(\omega \right)\right\}=\frac{{\omega }_0}{\pi }sin\text{c}\left(\frac{{\omega }_0}{\pi }\left(t-t_d\right)\right)\hfill & \left(10.31\right)\hfill \end{array}$$

For convenience, let us express hLP (t) using f0 instead of ω0:

$$\begin{array}{ll}{h}_{LP}\left(t\right)=2f_{0}sin\text{c}\left(2f_0\left(t-t_d\right)\right)\hfill & \left(10.32\right)\hfill \end{array}$$

The resulting impulse response is shown in Fig. 10.10.

A few observations are in order regarding the impulse response of the ideal lowpass filter:

1. The impulse response hLP (t) given by Eqn. (10.32) is in the form of a sinc function with a peak amplitude of 2f0. This peak occurs at the time instant t = td.
2. The zero crossings of the impulse response are uniformly spaced. Consecutive zero crossings are 1/(2f0) apart from each other.
3. The impulse response exists for all time instants including t < 0. Thus, the ideal lowpass filter is a non-causal system, and is therefore physically unrealizable.
4. Furthermore, no amount of additional time delay would make the impulse response in Eqn. (10.32) causal.

It is also interesting and instructive to derive the impulse response of an ideal bandpass filter the spectrum of which is shown in Fig. 10.7. We could compute the impulse response as the inverse Fourier transform of the system function in Eqn. (10.30), using the approach that was employed in finding the impulse response of the ideal lowpass filter. Instead we will take a simpler approach and express the ideal bandpass filter spectrum in terms of frequency shifts applied to the ideal lowpass filter spectrum. In order to accomplish this, we will initially consider an ideal lowpass filter with zero phase and no time delay, i.e., td = 0. A comparison of Figs. 10.6 and 10.7 reveals that the ideal bandpass filter spectrum with zero phase can be written as

$$\begin{array}{lll}{H}_{BP}\left(\omega \right)\hfill & =H_{LP}\left(\omega -{\omega }_b\right)+H_{LP}\left(\omega +{\omega }_b\right)\hfill & \hfill \\ \hfill & =\prod \left(\frac{\omega -{\omega }_b}{2{\omega }_0}\right)+\prod \left(\frac{\omega +{\omega }_b}{2{\omega }_0}\right)\hfill & \left(10.33\right)\hfill \end{array}$$

with parameter ωb set to the midpoint frequency of the passband

$$\begin{array}{ll}{\omega }_b=\frac{{\omega }_2+{\omega }_1}{2}\hfill & \left(10.34\right)\hfill \end{array}$$

and the cutoff frequency of the lowpass filter chosen as

$$\begin{array}{ll}{\omega }_0=\frac{{\omega }_2-{\omega }_1}{2}\hfill & \left(10.35\right)\hfill \end{array}$$

The construction of the bandpass filter spectrum from the lowpass filter spectrum is illustrated in Fig. 10.11.

Using the frequency-shifting property of the Fourier transform, we are now ready to write the impulse response of the ideal bandpass filter as

$$\begin{array}{lll}{h}_{BP}\left(t\right)\hfill & =h_{LP}\left(t\right)e^{j2\pi f_{b}t}+h_{LP}\left(t\right)e^{-j2\pi f_{b}t}\hfill & \hfill \\ \hfill & =2h_{LP}\left(t\right)cos\left(2\pi f_{b}t\right)\hfill & \hfill \\ \hfill & =4f_{0}sin\text{c}\left(2f_{0}t\right)cos\left(2\pi f_{b}t\right)\hfill & \left(10.36\right)\hfill \end{array}$$

For convenience we have used frequencies in Hertz rather than the radian frequencies in Eqn. (10.36). This impulse response is for an ideal bandpass filter with no time delay. It may be generalized by adding a linear phase term and the corresponding constant time delay to obtain

$$\begin{array}{ll}{h}_{BP}\left(t\right)=4f_{0}sin\text{c}\left[2f_0\left(t-t_d\right)\right]cos\left[2\pi f_b\left(t-t_d\right)\right]\hfill & \left(10.37\right)\hfill \end{array}$$

Finally, by substituting Eqns. (10.34) and (10.35) into Eqn. (10.37) we get

$$\begin{array}{ll}{h}_{BP}\left(t\right)=2\left(f_2-f_1\right)sin\text{c}\left[\left(f_2-f_1\right)\left(t-t_d\right)\right]cos\left[\pi \left(f_2+f_1\right)\left(t-t_d\right)\right]\hfill & \left(10.38\right)\hfill \end{array}$$

which is shown graphically in Fig. 10.12.

10.4 Design of Analog Filters

In this section we will focus on designing analog filters to meet a set of user specifications. We have seen in Section 10.3 that filters with ideal characteristics are not realizable since they are non-causal. In practical applications of filtering we seek realizable approximations to ideal filter behavior. The requirements on the frequency spectrum of the filter must be relaxed in order to obtain realizable filters.

In the design of analog filters, specifications for the desired filter are usually given in terms of a set of tolerance limits for the magnitude response |H (ω)|. Specification diagrams for the four frequency-selective filter types, namely lowpass, highpass, bandpass and band-reject filters, are shown in Fig. 10.13.

Consider the lowpass filter specifications illustrated by Fig. 10.13(a). The frequency interval 0 < ω < ω1 is referred to as the passband of the filter. Within the passband the magnitude of the filter spectrum is required to remain in the range 1 − Δ1 ≤ |H (ω) ≤ 1. For frequencies greater than ω2, that is, in the stopband of the filter, the magnitude is required to stay below Δ2. The parameters Δ1 and Δ2 are called the passband tolerance and the stopband tolerance respectively.

The frequency range ω1 < ω < ω2 between the passband and the stopband is referred to as the transition band where there are no constraints on the magnitude behavior of the filter. The inclusion of a transition band along with the tolerance limits on passband and stopband behavior should allow us to come up with practically realizable designs. It should be noted that the impulse response h (t) is required to be real-valued. Recall from the discussion in Section 4.3.5 of Chapter 4 that if the impulse response is real then the system function is conjugate symmetric, and therefore the magnitude of the system function is an even function of ω. This is the reason why the specification diagrams in Fig. 10.13 display only the positive frequencies.

Specifications for a highpass filter, shown in Fig. 10.13(b) mirror those of the lowpass filter. Bandpass and band-reject filter specifications, shown in Fig. 10.13(c),(d) have two transition bands each.

Tolerance limits for a desired filter can also be given on a decibel (dB) scale as shown in Fig. 10.14.

The maximum value of the magnitude response is set as 0 dB. The maximum allowed dB passband ripple Rp and the minimum required dB stopband attenuation As are related to the parameters Δ1 and Δ2 as

$$\begin{array}{ll}{R}_p=20{\text{ log}}_{10}\left(\frac{1}{1-{\Delta }_1}\right)\hfill & \left(10.39\right)\hfill \end{array}$$

and

$$\begin{array}{ll}{A}_s=20{\text{ log}}_{10}\left(\frac{1}{{\Delta }_2}\right)\hfill & \left(10.40\right)\hfill \end{array}$$

Given the specifications of the desired filter in the form of one of the specification diagrams in Fig. 10.13, the design problem can be stated as follows:

In general, the analog filter design problem begins with the design of a lowpass prototype filter regardless of the type of filter that is ultimately needed. If a highpass, bandpass or band-reject filter is desired, it is obtained subsequently by means of a frequency transformation to be applied to the lowpass prototype filter. This approach is motivated by the availability of well established methods for the design of analog lowpass filters.

A number of approximation formulas exist in the literature for specifying the squared-magnitude of a realizable lowpass filter spectrum. Some of the better known formulas will be mentioned here.

Butterworth lowpass filters are characterized by the squared-magnitude function

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=\frac{1}{1+{\left(\omega /{\omega }_c\right)}^{2N}}\hfill & \left(10.41\right)\hfill \end{array}$$

where the parameters N and ωc are the filter-order and the cutoff frequency respectively.

As an alternative to the Butterworth approximation formula, the Chebyshev type-I approximation formula for the squared-magnitude function of a lowpass filter is

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=\frac{1}{1+{\varepsilon }^2{C}_N^2\left(\omega /{\omega }_1\right)}\hfill & \left(10.42\right)\hfill \end{array}$$

where ε is a positive constant. The frequency ω1 is the passband edge frequency. The function CN (ν) in the denominator represents the Chebyshev polynomial of order N.

Chebyshev type-II approximation formula for the squared-magnitude function is similar to the Chebyshev type-I formula.

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=\frac{{\varepsilon }^2{C}_N^2\left({\omega }_2/\omega \right)}{1+{\varepsilon }^2{C}_N^2\left({\omega }_2/\omega \right)}\hfill & \left(10.43\right)\hfill \end{array}$$

The parameter ω2 is the stopband edge frequency.

Finally, an elliptic lowpass filter is characterized by the squared-magnitude function

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=\frac{1}{1+{\varepsilon }^2{\psi }_N^2\left(\omega /{\omega }_1\right)}\hfill & \left(10.44\right)\hfill \end{array}$$

The parameter ε is a positive constant. The function φN (ν) is called a Chebyshev rational function, and is defined in terms of Jacobi elliptic functions.

In the rest of this section we will provide a brief introduction to analog filter design methods. A thorough treatment of analog filter design would be well beyond the scope of this text. There are a number of excellent texts on the subject.

10.4.1 Butterworth lowpass filters

Consider again the Butterworth squared-magnitude function given by Eqn. (10.41). At the frequency ω = ωc the magnitude is equal to $1/\sqrt{2}$. Since this approximately corresponds to the −3 dB point on a logarithmic scale, the parameter ωc is also referred to as the 3-dB cutoff frequency of the Butterworth filter. Furthermore, it can easily be shown that the first 2N − 1 derivatives of the magnitude spectrum |H (ω)| are equal to zero at ω = 0, that is,

$$\begin{array}{lll}{\frac{d^n}{d{\omega }^n}\left\{\left|H\left(\omega \right)\right|\right\}|}_{\omega =0}=0,\hfill & n=1,...,2N-1\hfill & \left(10.45\right)\hfill \end{array}$$

The magnitude characteristic |H (ω)| for a Butterworth filter is said to be maximally flat. Magnitude spectra for Butterworth filters of various orders are shown in Fig. 10.15.

We know that the squared-magnitude function can be expressed as the product of the system function and its complex conjugate, that is

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=H\left(\omega \right)H^{*}\left(\omega \right)\hfill & \left(10.46\right)\hfill \end{array}$$

In addition, since the impulse response h (t) of the filter is real-valued, its system function exhibits conjugate symmetry (see Section 4.3.5 of Chapter 4):

$$\begin{array}{ll}{H}^{*}\left(\omega \right)=H\left(-\omega \right)\hfill & \left(10.47\right)\hfill \end{array}$$

Combining Eqns. (10.46) and (10.47) yields

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=H\left(\omega \right)H\left(-\omega \right)\hfill & \left(10.48\right)\hfill \end{array}$$

Using the s-domain system function H (s), the problem of designing a Butterworth filter can be stated as follows: Given the values of the two filter parameters ωc and N, find the s-domain system function H (s) for which the squared-magnitude of the function H (ω) matches the right side of Eqn. (10.41).

For a system function H (s) with real coefficients, it can be shown that

$$\begin{array}{ll}{{\left|H\left(\omega \right)\right|}^2=H\left(s\right)H\left(-s\right)|}_{s=j\omega }\hfill & \left(10.49\right)\hfill \end{array}$$

The procedure in Eqn. (10.49) needs to be reversed to find the system function H (s). Substituting s for jω, or equivalently, −s2 for ω2 in Eqn. (10.41) we obtain

$$\begin{array}{ll}H\left(s\right)H\left(-s\right)=\frac{1}{1+{\left(-s^2/{\omega }_c^2\right)}^N}\hfill & \left(10.50\right)\hfill \end{array}$$

If H (s) has a pole in the left half s-plane at p1 = σ1 + jω1, then H (−s) has a corresponding pole at ${\overline{p}}_1=-{\sigma }_1-j{\omega }_1$. Thus the poles of H (s) are mirror images of the poles of H (−s) with respect to the origin. In order to extract H (s) from the product H (s) H (−s), we need to find the poles of this product, and separate the two sets of poles that belong to H (s) and H (−s) respectively. The cases of even and odd filter-orders N need to be treated separately.

A few observations are in order:

1. All the poles of the product H (s) H (−s) are located on a circle with radius equal to ωc.
2. Furthermore, the poles are equally spaced on the circle, and the angle between any two adjacent poles is π/N radians.
3. Since H (s) and H (−s) have only real coefficients, all complex poles appear in conjugate pairs.

Fig. 10.16 depicts the situation for ωc = 2 rad/s, and N = 5.

We are interested in obtaining a filter that is both causal and stable. It is therefore necessary to use the poles in the left half of the s-plane for H (s). The remaining poles, the ones in the right half of the s-plane, must be associated with H (−s). The system function H (s) for the Butterworth lowpass filter is constructed in the form

$$\begin{array}{llll}H\left(s\right)=\frac{{\omega }_c^N}{\prod _k\left(s-p_k\right)},\hfill & \text{for all }k\text{ that satisfy}\hfill & \frac{\pi }{2}<\frac{k\pi }{N}<\frac{3\pi }{2}\hfill & \left(10.54\right)\hfill \end{array}$$

For even filter orders, poles of the magnitude-squared system function are still equally spaced on a circle. Compared to odd filter orders, the only difference in this case is that no poles appear on the real axis. The system function H (s) can be constructed in the same way as it was done for odd values of N.

Example 10.2: Butterworth filter design

Determine the system function H (s) of a third-order Butterworth lowpass filter with a 3-dB cutoff frequency of ωc = 20π rad/s.

Solution: Using the Butterworth approximation formula in Eqn. (10.41), the squared-magnitude function of the desired filter can be written as

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=\frac{1}{1+{\left(\omega /20\pi \right)}^6}\hfill & \left(10.58\right)\hfill \end{array}$$

Substituting ω2 → −s2 in Eqn. (10.58) we obtain

$$H\left(s\right)H\left(-s\right)=\frac{1}{1+{\left(-s^2/400{\pi }^2\right)}^3}$$

or equivalently

$$H\left(s\right)H\left(-s\right)=\frac{{\left(20\pi \right)}^6}{{\left(20\pi \right)}^6-s^6}$$

The poles of the function H (s) H (−s) are those values of s that satisfy

$$s^6={\left(20\pi \right)}^6$$

Solutions of Eqn. (10.59) are

$$\begin{array}{ll}{p}_k=\left(20\pi \right)e^{j2\pi k/6},\hfill & k=0,...,5\hfill \end{array}$$

as shown in Fig. 10.17.

We are particularly interested in the poles that lie in the left half of the s-plane since those are the poles that are associated with H (s). By inspection of Fig. 10.17 it is clear that the poles of H (s) are

$$\begin{array}{l}{p}_2=\left(20\pi \right)e^{j2\pi /3}\hfill \\ p_3=\left(20\pi \right)e^{j\pi }=-20\pi \hfill \\ p_4=\left(20\pi \right)e^{j4\pi /3}=\left(20\pi \right)e^{-j2\pi /3}\hfill \end{array}$$

We are now ready to construct the transfer function by using these poles:

$$\begin{array}{lll}H\left(s\right)\hfill & =\frac{{\left(20\pi \right)}^3}{\left(s-p_2\right)\left(s-p_3\right)\left(s-p_4\right)}\hfill & \hfill \\ \hfill & =\frac{{\left(20\pi \right)}^3}{\left(s+20\pi \right)\left(s^2+20\pi s+400{\pi }^2\right)}\hfill & \left(10.59\right)\hfill \end{array}$$

In order to compute the magnitude and the phase characteristics of the system we need to substitute s = jω into Eqn. (10.59) to obtain

$$H\left(\omega \right)=\frac{{\left(20\pi \right)}^3}{\left(j\omega +20\pi \right)\left(-{\omega }^2+j20\pi \omega +400{\pi }^2\right)}$$

The magnitude of the system function is

$$\left|H\left(\omega \right)\right|=\frac{{\left(20\pi \right)}^3}{\sqrt{{\omega }^2+{\left(20\pi \right)}^2}\sqrt{{\left(-{\omega }^2+400{\pi }^2\right)}^2+{\left(20\pi \omega \right)}^2}}$$

and the phase of the system function is

$$\measuredangle H\left(\omega \right)=-{tan}^{-1}\left(\frac{\omega }{20\pi }\right)-{tan}^{-1}\left(\frac{20\pi \omega }{-{\omega }^2+400{\pi }^2}\right)$$

Magnitude and phase responses of the designed filter are shown in Fig. 10.18.

Software resources:

ex_10_2.m

Obtaining N and ωc for the Butterworth lowpass filter

The desired behavior of a lowpass filter to be designed is usually specified in terms of the two critical frequencies ω1 and ω2 as well as the decibel tolerance values Rp and As as shown in Fig. 10.14. The corresponding parameter values for ωc and N need to be determined from the provided set of specifications so that the filter can be designed. If Rp = 3 dB, then we have ωc = ω1. Otherwise the requirements at the two critical frequencies can be expressed as inequalities.

In summary, the parameters of the Butterworth lowpass filter are computed as follows:

10.4.2 Chebyshev lowpass filters

The squared-magnitude function for a Chebyshev lowpass filter is given by

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=\frac{1}{1+{\varepsilon }^2{C}_N^2\left(\omega /{\omega }_1\right)}\hfill & \left(10.66\right)\hfill \end{array}$$

Parameter ε is a positive constant, N is the order of the filter, and ω1 is the passband edge frequency in rad/s (see Fig. 10.13(a)). The term CN (ν) represents the Chebyshev polynomial of order N. Chebyshev filters of this type are sometimes called Chebyshev type-I filters. A variant of the squared-magnitude response in Eqn. (10.66) will be introduced in the next section as the inverse Chebyshev filter or the Chebyshev type-II filter.

Once the values of the parameters ε, N and ω1 are specified, the design procedure proceeds as follows:

1. Find CN (ν), the Chebyshev polynomial of order N.
2. Substitute ν = ω/ω1 in the polynomial CN (ν).
3. Form the Chebyshev type I squared-magnitude function as specified in Eqn. (10.66).
4. Find the poles of the system function H (s) from the squared-magnitude function. As in the case of the Butterworth lowpass filter, this involves substituting jω = s into the squared-magnitude function, finding the roots of the denominator polynomial, and separating those poles into two sets, one to be associated with H (s) and the other to be associated with H (s).

Before we carry out the steps outlined above, it will be useful to review the properties of Chebyshev polynomials.

Chebyshev polynomials

The Chebyshev polynomial of order N is defined as

$$\begin{array}{ll}{C}_N\left(\nu \right)=cos\left(N{cos}^{-1}\left(\nu \right)\right)\hfill & \left(10.67\right)\hfill \end{array}$$

A better approach for understanding the definition in Eqn. (10.67) would be to split it into two equations as

$$\begin{array}{ll}\nu =cos\left(\theta \right)\hfill & \left(10.68\right)\hfill \end{array}$$

and

$$\begin{array}{ll}{C}_N\left(\nu \right)=cos\left(N\theta \right)\hfill & \left(10.69\right)\hfill \end{array}$$

The Chebyshev polynomial of a specified order N can be obtained by

1. Using trigonometric identities to express cos (Nθ) as a function of cos (θ).
2. Replacing each cos (θ) term with ν.

The first two polynomials are easy to obtain from this definition:

$$\begin{array}{l}{C}_0\left(\nu \right)=1\\ C_1\left(\nu \right)=\nu \end{array}$$

To obtain C2 (ν) let us write cos (2θ) in terms of cos (θ) using a trigonometric identity:

$$cos\left(2\theta \right)=2{cos}^2\left(\theta \right)-1$$

Thus, the second-order Chebyshev polynomial is

$$C_2\left(\nu \right)=2{\nu }^2-1$$

Higher-order Chebyshev polynomials may be obtained by continuing in this fashion. As the order increases, however, the procedure outlined above becomes increasingly tedious. Fortunately, it is possible to derive a recursive formula to facilitate the derivation of higher-order Chebyshev polynomals.

Thus, the recursive formula for obtaining Chebyshev polynomials is

$$\begin{array}{ll}{C}_{N+1}\left(\nu \right)=2\nu C_N\left(\nu \right)-C_{N-1}\left(\nu \right)\hfill & \left(10.70\right)\hfill \end{array}$$

The recursive formula in Eqn. (10.70) allows any order Chebyshev polynomial to be found provided that the polynomials for the previous two orders are known. With the knowledge of C0 (ν) = 1 and C1 (ν) = ν, the polynomial C2 (ν) can be found as

$$\begin{array}{ll}{C}_2\left(\nu \right)\hfill & =2\nu C_1\left(\nu \right)-C_0\left(\nu \right)\hfill \\ \hfill & =2{\nu }^2-1\hfill \end{array}$$

Similarly C3 (ν) is found as

$$\begin{array}{ll}{C}_3\left(\nu \right)\hfill & =2\nu C_2\left(\nu \right)-C_1\left(\nu \right)\hfill \\ \hfill & =2\nu \left(2{\nu }^2-1\right)-\nu \hfill \\ \hfill & =4{\nu }^3-3\nu \hfill \end{array}$$

The behaviors of the first few orders of Chebyshev polynomials are illustrated in Fig. 10.19.

A few interesting characteristics of Chebyshev polynomials can be readily observed:

1. All Chebyshev polynomials pass through the point (1, 1), that is, CN (1) = 1 for all N.
2. For all Chebyshev polynomials, if |ν| ≤ 1 then |CN (ν)| ≤ 1.
3. For |ν| > 1 all Chebyshev polynomials grow monotonically without bound.
4. At ν = 0 the behavior of Chebyshev polynomials depends on the order N · CN (0) = 0 if N is odd, and CN (0) = ±1 if N is even.

Fig. 10.20 shows the squared Chebyshev polynomials for the first few orders N.

Squared Chebyshev polynomials $C_N^2\left(\nu \right)$ oscillate in the amplitude range [0, 1] for |ν| ≤ 1 and N ≥ 2. For |ν| > 1 they grow monotonically without bound.

The squared magnitude function given by Eqn. (10.66) is shown in Fig. 10.21 for several values of N using ε = 0.4 and ω1 = 1 rad/s.

Poles for the Chebyshev lowpass filter

In constructing the system function for a Chebyshev lowpass filter we will use an approach similar to that taken with a Butterworth lowpass filter in Section 10.4.1. First we will determine the poles of the product H (s) H (−s). Recall that the product H (s) H (−s) is obtained by starting with the squared magnitude function |H (ω)|2 and replacing jω with s:

$$\begin{array}{ll}H\left(s\right)H\left(-s\right)=\frac{1}{1+{\varepsilon }^2{C}_N^2\left(\frac{s}{j{\omega }_1}\right)}\hfill & \left(10.71\right)\hfill \end{array}$$

The poles pk of H (s) H (−s) are the solutions of the equation

$$\begin{array}{ll}1+{\varepsilon }^2{C}_N^2\left(\frac{p_k}{j{\omega }_1}\right)=0\hfill & \left(10.72\right)\hfill \end{array}$$

for k = 0,..., 2N − 1.

Using the values of αk and βk found, the poles of H (s) H (−s) are

$$\begin{array}{ll}{p}_k=j{\omega }_1\left[cos\left({\alpha }_k\right)cosh\left({\beta }_k\right)-jsin\left({\alpha }_k\right)sinh\left({\beta }_k\right)\right]\hfill & \left(10.86\right)\hfill \end{array}$$

for k = 0,..., 2N − 1. It can be shown that the poles given by Eqn. (10.86) are on an elliptical trajectory. The ones in the left half s-plane are associated with H (s) in order to obtain a causal and stable filter.

Example 10.3: Design of a Chebyshev lowpass filter

Design a third-order Chebyshev type-I analog lowpass filter with a passband edge frequency of ω1 = 1 rad/s and ε = 0.4. Afterwards compute and graph the magnitude response of the designed filter.

Solution: Using Eqn. (10.82) the parameter αk is found as

$$\begin{array}{cc}{\alpha }_k=\frac{\left(2k+1\right)\pi }{6},& k=0,...,5\end{array}$$

and the parameter βk is obtained from Eqn. (10.84) as

$${\beta }_k=\frac{{sinh}^{-1}\left(1/0.4\right)}{3}=0.5491$$

Poles of H (s) H (−s) are found through the use of Eqn. (10.86):

$$\begin{array}{cc}{p}_k=j\left[cos\left(\frac{\left(2k+1\right)\pi }{6}\right)cosh\left(0.5491\right)-jsin\left(\frac{\left(2k+1\right)\pi }{6}\right)sinh\left(0.5491\right)\right],k=0,...,5& \left(10.87\right)\end{array}$$

Evaluating Eqn. (10.87) for k = 0,..., 5 yields the following poles:

$$\begin{array}{ll}{p}_0\hfill & =0.2885+j1\hfill \\ p_1\hfill & =0.5771\hfill \\ p_2\hfill & =0.2885-j1\hfill \\ p_3\hfill & =-0.2885-j1\hfill \\ p_4\hfill & =-0.5771\hfill \\ p_5\hfill & =-0.2885+j1\hfill \end{array}$$

The first three poles, namely p0, p1 and p2 are in the right half s-plane; therefore they belong to H (−s). The system function H (s) should be constructed using the remaining three poles.

$$\begin{array}{ll}H\left(s\right)\hfill & =\frac{0.6250}{\left(s+0.2885+j1\right)\left(s+0.5771\right)\left(s+0.2885-j1\right)}\hfill \\ \hfill & =\frac{0.6250}{s^3+1.1542s^2+1.4161s+0.6250}\hfill \end{array}$$

The numerator of H (s) was adjusted to achieve |H (0)| = 0. The magnitude and the phase of H (s) are shown in Fig. 10.22.

Software resources:

ex_10_3.m

Obtaining N and ε for the Chebyshev lowpass filter

In some design problems it is desired to find the lowest-order filter that satisfies a set of design specifications. The desired behavior of the lowpass filter is typically specified in terms of the two critical frequencies ω1 and ω2 as well as the dB tolerance values Rp and As for the passband and the stopband respectively.

In summary, the parameters of a Chebyshev type-I lowpass filter are computed as follows:

10.4.3 Inverse Chebyshev lowpass filters

The squared-magnitude function for an inverse Chebyshev lowpass filter, also referred to as a Chebyshev type-II lowpass filter, is

$$\begin{array}{ll}{\left|H\left(\omega \right)\right|}^2=\frac{{\varepsilon }^2{C}_N\left({\omega }_2/\omega \right)}{1+{\varepsilon }^2{C}_N^2\left({\omega }_2/\omega \right)}\hfill & \left(10.102\right)\hfill \end{array}$$

As in the case of Chebyshev type-I filters, ε is a positive constant and CN (ν) represents the Chebyshev polynomial of order N. The parameter ω2 is the stopband edge frequency (see Fig. 10.13(a)). Magnitude responses of several Chebyshev type-II filters are shown in Fig. 10.23. The magnitude response of the Chebyshev type-II filter is smooth in the passband and has equiripple behavior in the stopband.

Poles and zeros of the inverse Chebyshev filter

The denominator of the squared magnitude function |H (ω)|2 in Eqn. (10.102) is very similar to that of the type-I Chebyshev filter squared magnitude response given by Eqn. (10.66). Consequently, most of the results obtained in the previous section through the derivation of the poles of Chebyshev type-I lowpass filter will be usable. For an inverse Chebyshev filter the poles pk of the product H (s) H (−s) are the solutions of the equation

$$\begin{array}{ll}1+{\varepsilon }^2{C}_N^2\left(\frac{j{\omega }_2}{p_k}\right)=0\hfill & \left(10.103\right)\hfill \end{array}$$

for k = 0,..., 2N − 1.

A summary of the results is given below:

Obtaining N and ε for the inverse Chebyshev lowpass filter

The inverse Chebyshev lowpass filter is uniquely described by the choice of parameters ω2, ε and N. If the desired filter is described by the specification diagram of Fig. 10.14 then we have the set of parameters Rp, As, ω1 and ω2 from which ε and N must be determined.

The first step is to choose the parameter ε so that the magnitude response of the filter fills the allowed area completely in the stopband.

In summary, the parameters of a Chebyshev type-II lowpass filter are computed as follows:

10.4.4 Analog filter transformations

In preceding sections we have discussed the design of analog lowpass filters. Design formulas for Butterworth, Chebyshev and elliptic squared-magnitude functions were presented only for lowpass filters. If a different filter type such as a highpass, bandpass or band-reject filter is needed, it is obtained from a lowpass filter by means of an analog filter transformation.

Let G (s) be the system function of an analog lowpass filter, and let H (λ) represent the new filter to be obtained from it. For the new filter we use λ as the Laplace transform variable. H (λ) is obtained from G (s) through a transformation such that

$$\begin{array}{cc}{H\left(\lambda \right)=G\left(s\right)|}_{s=f\left(\lambda \right)}& \left(10.119\right)\end{array}$$

The function s = f (λ) is the transformation that converts the lowpass filter into the type of filter desired.

Lowpass to highpass transformation

Consider the specification diagrams shown in Fig. 10.24 for a lowpass filter with magnitude response |G (ω)| and a highpass filter with magnitude response |H (ω)|.

It is desired to obtain the highpass filter system function H (λ) from the lowpass filter system function G (s). The transformation

$$\begin{array}{cccc}\frac{s}{{\omega }_{L1}}=\frac{{\omega }_{H2}}{\lambda }& \Rightarrow & s=\frac{{\omega }_{L1}{\omega }_{H2}}{\lambda }=\frac{{\omega }_0}{\lambda }& \left(10.120\right)\end{array}$$

can be used for this purpose The magnitudes of the two filters are identical at their respective passband edge frequencies, that is,

$$\begin{array}{cc}\left|H\left({\omega }_{H2}\right)\right|=\left|G\left({\omega }_{L1}\right)\right|& \left(10.121\right)\end{array}$$

The stopband edges of the two filters are related by

$$\begin{array}{cc}{\omega }_{L2}{\omega }_{H1}={\omega }_{L1}{\omega }_{H2}& \left(10.122\right)\end{array}$$

so that we have

$$\begin{array}{cc}\left|H\left({\omega }_{H1}\right)\right|=\left|G\left({\omega }_{L2}\right)\right|& \left(10.123\right)\end{array}$$

Example 10.4: Design of a Chebyshev highpass filter

Recall that a third-order Chebyshev type-I analog lowpass filter was designed in Example 10.3 with a cutoff frequency of ω1 = 1 rad/s and ε = 0.4. The system function of the filter was found to be

$$G\left(s\right)=\frac{0.6250}{s^3+1.1542s^2+1.4161s+0.6250}$$

Convert this filter to a highpass filter with a critical frequency of ωH2 = 3 rad/s. Afterwards compute and graph the magnitude response of the designed filter.

Solution: The critical frequency of the lowpass filter is ωL1 = 1 rad/s; therefore, we will use the transformation

$$s=\frac{3}{\lambda }$$

which leads to the system function

$$\begin{array}{cc}\hfill H\left(\lambda \right)& =\frac{0.6250}{{\left(\frac{3}{\lambda }\right)}^3+1.1542{\left(\frac{3}{\lambda }\right)}^2+1.4161\left(\frac{3}{\lambda }\right)+0.6250}\hfill \\ \hfill & =\frac{0.6250{\lambda }^3}{0.6250{\lambda }^3+4.2482{\lambda }^2+10.3875\lambda +27}\hfill \end{array}$$

The magnitude and the phase of H (λ) are shown in Fig. 10.25.

Software resources:

ex_10_4.m

Lowpass to bandpass transformation

Specification diagrams in Fig. 10.26 show a lowpass filter with magnitude response |G (ω)| and a bandpass filter with magnitude response |H (ω)|.

It is desired to obtain the bandpass filter system function H (λ) from the lowpass filter system function G (s). The transformation

$$\begin{array}{cc}s=\frac{{\lambda }^2+{\omega }_0^2}{B\lambda }& \left(10.124\right)\end{array}$$

can be used for this purpose. We require

$$\begin{array}{cc}\left|H\left({\omega }_{B2}\right)\right|=\left|G\left(-{\omega }_{L1}\right)\right|& \left(10.125\right)\end{array}$$

and

$$\begin{array}{cc}\left|H\left({\omega }_{B3}\right)\right|=\left|G\left({\omega }_{L1}\right)\right|& \left(10.126\right)\end{array}$$

It can be shown that (see Problem 10.10) these requirements are met if the parameters ω 0 and B are chosen as

$$\begin{array}{cccc}{\omega }_0=\sqrt{{\omega }_{B2}{\omega }_{B3}}& \text{and}& B=\frac{{\omega }_{B3}-{\omega }_{B2}}{{\omega }_{L1}}& \left(10.127\right)\end{array}$$

The parameter ω0 is the geometric mean of the passband edge frequencies of the bandpass filter. The parameter B is the ratio of the bandwidth of the bandpass filter to the bandwidth of the lowpass filter.

Lowpass to band-reject transformation

Specification diagrams in Fig. 10.27 show a lowpass filter with magnitude response |G (ω)| and a band-reject filter with magnitude response |H (ω)|.

In order to obtain the band-reject filter system function H (λ) from the lowpass filter system function G (s), the transformation to be used is in the form

$$\begin{array}{cc}s=\frac{B\lambda }{{\lambda }^2+{\omega }_0^2}& \left(10.128\right)\end{array}$$

For a lowpass filter to be converted into a highpass filter we require

$$\begin{array}{cc}\left|H\left({\omega }_{S4}\right)\right|=\left|G\left(-{\omega }_{L1}\right)\right|& \left(10.129\right)\end{array}$$

and

$$\begin{array}{cc}\left|H\left({\omega }_{S1}\right)\right|=\left|G\left(-{\omega }_{L1}\right)\right|& \left(10.130\right)\end{array}$$

It can be shown that (see Problem 10.12) these requirements are met if the parameters ω0 and B are chosen as

$$\begin{array}{cccc}{\omega }_0=\sqrt{{\omega }_{S1}{\omega }_{S4}}& \text{and}& B=\left({\omega }_{S4}-{\omega }_{S1}\right){\omega }_{L1}& \left(10.131\right)\end{array}$$

10.5 Design of Digital Filters

In this section we will briefly discuss the design of discrete-time filters. Discrete-time filters are viewed under two broad categories: infinite impulse response (IIR) filters, and finite impulse response (FIR) filters. The system function of an IIR filter has both poles and zeros. Consequently the impulse response of the filter is of infinite length. It should be noted, however, that the impulse response of a stable IIR filter must also be absolute summable, and must therefore decay over time. In contrast, the behavior of an FIR filter is controlled only by the placement of the zeros of its system function. For causal FIR filters all of the poles are at the origin, and they do not contribute to the magnitude characteristics of the filter.

For a given set of specifications, an IIR filter is generally more efficient than a comparable FIR filter. On the other hand, FIR filters are always stable. Additionally, a linear phase characteristic is possible with FIR filters whereas causal and stable IIR filters cannot have linear phase. The significance of a linear phase characteristic is that the time delay is constant for all frequencies. This is desirable in some applications, and requires the use of an FIR filter.

10.5.1 Design of IIR filters

Consider again the specification diagram for an analog lowpass filter shown in Fig. 10.13(a). It can be adapted for the system function of a discrete-time filter as shown in Fig. 10.28.

The magnitude specifications of the discrete-time filter are given in the angular frequency range 0 ≤ Ω < π. Since the impulse response of the filter to be designed is real, the magnitude is an even function of Ω, and therefore the part of the magnitude response in the range −π ≤ Ω < 0 is obtained as the mirror image of what is shown in Fig. 10.28. Furthermore, the magnitude response is 2π-periodic.

Specification diagrams for other filter types (highpass, bandpass or band-reject) can be adapted to discrete-time filters in a similar manner.

The most common method of designing IIR filters is to start with an appropriate analog filter, and to convert its system function to the system function of a discrete-time filter by means of some transformation. Designing a discrete-time filter by this approach involves a three step procedure:

Step 2 of the design procedure, namely the design of an analog prototype, was discussed in Section 10.4. The problem of converting an analog system function to a discrete-time system function (step 3) will be discussed next, followed by the conversion of specifications (step 1). The reason for this particular order is because step 1 is dependent on which method is to be used in step 3.

Analog to discrete-time conversion

Step 3 of the design procedure discussed above requires a discrete-time system function to be obtained from the analog system function. The goal is to come up with a system function H (z) the characteristics of which resemble those of the analog system function G (s) in some way. The goal is somewhat loosely defined; there is no clear definition of when two filters resemble each other. There is also the question of which characteristics of the two filters (magnitude, phase, time-delay, impulse response, etc.) exhibit similarity. In most cases analog prototype filters are designed based on their magnitude responses. Therefore, our goal will be to obtain a reasonable degree of similarity between the magnitude characteristics of the two filters.

The main challenge is the fact that the system function H (Ω) is 2π-periodic whereas the analog prototype filter system function G (ω) is not. Consequently, an exact match between the system functions is not possible. Multiple methods exist for obtaining an approximate match. Two of these methods, namely impulse invariance and bilinear transformation, will be discussed in the remainder of this section.

Impulse invariance

A possible method of obtaining a discrete-time filter from an analog prototype is to ensure that the impulse response is preserved in the conversion process. Specifically, the discrete-time filter is chosen such that its impulse response is equal to a scaled and sampled version of the impulse response of the analog prototype:

$$\begin{array}{cc}h\left[n\right]=Tg\left(nT\right)& \left(10.132\right)\end{array}$$

The relationship in Eqn. (10.132) is the sampling relationship discussed in Chapter 6. The reason for the scale factor T will become obvious when we consider the magnitude spectrum of the resulting discrete-time filter.

Let the system function for a causal analog prototype filter be G (s) which can be expanded into partial fractions in the form

$$\begin{array}{cc}G\left(s\right)=\sum _{i=1}^N\frac{k_i}{s-p_i}& \left(10.133\right)\end{array}$$

where pi are the poles, and ki are the corresponding residues. The impulse response is

$$\begin{array}{cc}g\left(t\right)=\sum _{i=1}^N{k}_i{e}^{p_{i}t}u\left(t\right)& \left(10.134\right)\end{array}$$

Under impulse-invariant design the impulse response of the discrete-time system is required to be

$$\begin{array}{cc}h\left[n\right]=Tg\left(nT\right)=\sum _{i=1}^{N}T{k}_i{e}^{p_{i}nT}u\left[n\right]& \left(10.135\right)\end{array}$$

The system function of the discrete-time filter is found by taking the z transform of h [n] as

$$\begin{array}{cc}H\left(z\right)=\sum _{i=1}^N\frac{T{k}_{i}z}{z-e^{p_{i}T}}& \left(10.136\right)\end{array}$$

Comparing Eqns. (10.133) and (10.136) we conclude that H (z) can be obtained directly from the partial fraction expansion of G (z) without the intermediate steps of computing g (t) and sampling it.

An important issue in the conversion of an analog filter to a discrete-time filter is the requirement that a stable discrete-time filter be obtained from a stable analog filter. The poles of the analog prototype filter are at

$$\begin{array}{cc}\begin{array}{cc}{p}_i={\sigma }_i+j{\omega }_i,& i=1,...,N\end{array}& \left(10.137\right)\end{array}$$

The corresponding poles of the discrete-time filter obtained through impulse invariance are at

$$\begin{array}{ccc}{\overline{p}}_i=e^{p_{i}T}=e^{{\sigma }_{i}T}{e}^{j{\omega }_{i}T},& i=1,...,N& \left(10.138\right)\end{array}$$

For a casual and stable analog prototype filter the poles are in the left half s-plane, and therefore σi < 0 for i = 1,..., N. This implies that

$$\begin{array}{ccc}\left|{\overline{p}}_i\right|=e^{{\sigma }_{i}T}<1,& i=1,...,N& \left(10.139\right)\end{array}$$

We are assured that a casual and stable analog prototype filter leads to a casual and stable discrete-time filter when the impulse invariance technique is used.

Finally we need to understand the relationship between the magnitude responses of the analog prototype filter and the resulting discrete-time filter. This step is easy owing to the sampling relationship between g (t) and h [n] given by Eqn. (10.132). Using Eqn. (6.25) from Chapter 6 we obtain

$$\begin{array}{cc}H\left(\Omega \right)=\sum _{k=-\infty }^{\infty }G\left(\frac{\Omega -2\pi k}{T}\right)& \left(10.140\right)\end{array}$$

In order to avoid aliasing, the analog prototype filter must be bandlimited such that

$$\begin{array}{ccc}G\left(\omega \right)=0,& \left|\omega \right|\ge \frac{\pi }{T}& \left(10.141\right)\end{array}$$

so that

$$\begin{array}{ccc}H\left(\Omega \right)=G\left(\frac{\Omega }{T}\right)& -\pi \le \omega <\pi & \left(10.142\right)\end{array}$$

and the spectrum of the resulting discrete-time filter is a frequency-scaled version of the analog prototype filter spectrum in the range −π ≤ < ω < π. Nevertheless, the bandlimiting condition given by Eqn. (10.141) is often not satisfied by the analog prototype filter (none of the analog designs considered in Section 10.4 are perfectly bandlimited). Therefore, the spectrum of the resulting discrete-time filter is an aliased version of the spectrum of the analog prototype. This is illustrated in Fig. 10.29.

Based on Eqn. (10.140) the value of the system function H (Ω) at a specific frequency Ω = Ω0 is determined by contributions from the analog prototype system function G (ω) at frequencies

$$\begin{array}{ccc}{\omega }_k=\frac{{\Omega }_0-2\pi k}{T},& \text{all }k& \left(10.143\right)\end{array}$$

Due to the aliasing of the spectrum, impulse-invariant designs are only useful for lowpass filters for which aliasing can be kept at acceptable levels.

Example 10.5: Impulse invariant design of a lowpass filter

Consider the third-order Chebyshev type-II analog lowpass filter designed in Example 10.3. Convert this filter to a discrete-time filter using the impulse invariance technique with T = 0.2 s. Afterwards compute and graph the magnitude response of the discrete-time filter.

Solution: The system function of the analog prototype filter designed in Example 10.3 was

$$G\left(s\right)=\frac{0.6250}{s^3+1.1542s^2+1.4161s+0.6250}$$

which can be written in partial fraction form as

$$G\left(s\right)=\frac{-0.2885-j0.0833}{s+0.2886-j1}+\frac{-0.2885+j0.0833}{s+0.2886+j1}+\frac{0.5771}{s+0.5771}$$

Using Eqn. (10.136) the discrete-time filter system function is written in partial fraction form as

$$H\left(z\right)=\frac{\left(-0.0577-j0.0167\right)z}{z-0.9251-j0.1875}+\frac{\left(-0.0577+j0.016\right)z}{z-0.9251+j0.1875}+\frac{0.1154z}{z-0.8910}$$

The closed form expression for H (z) is

$$H\left(z\right)=\frac{0.0023z^2+0.0021z}{z^3-2.7412z^2+2.5395z-0.7939}$$

The magnitude and the phase of H (z) are shown in Fig. 10.30.

Note that aliasing can be kept at a negligible level with the choice of T, and the analog frequency ω1 = 1 rad/s corresponds to the discrete-time frequency Ω = 0.2 radians.

Software resources:

ex_10_5a.m

ex_10_5b.m

Example 10.5 demonstrates that, if the analog prototype filter system function G (s) is fixed, then the sampling interval should be chosen as small as possible to minimize the effect of aliasing. However, a typical IIR filter design problem begins with the critical frequencies Ω1 and Ω2 of the discrete-time filter. Let the corresponding critical frequencies of the analog prototype filter to be designed be ω1 and ω2 respectively. Recall from the sampling relationship that the analog frequencies and the discrete-time frequencies are related by ω = Ω/T. If the sampling interval is T, the sampling rate is ωs = 2π/T, and we obtain

$$\begin{array}{cc}\frac{{\omega }_1}{{\omega }_s}=\frac{{\Omega }_1}{2\pi }=\text{constant}& \left(10.144\right)\end{array}$$

and

$$\begin{array}{cc}\frac{{\omega }_2}{{\omega }_s}=\frac{{\Omega }_2}{2\pi }=\text{constant}& \left(10.145\right)\end{array}$$

For specified discrete-time critical frequencies, increasing the sampling rate ωs does not improve the aliasing condition since the bandwidth of the signal g (t) to be sampled also increases proportionally. Therefore the choice of T is irrelevant and we often use T = 1 s for simplicity. Example 10.6 will illustrate this.

Example 10.6: Impulse invariant design specifications

A lowpass IIR filter is to be designed with the following specifications:

$$\begin{array}{cccc}{\Omega }_1=0.2\pi ,& {\Omega }_2=0.25\pi ,& R_p=1\text{ dB},& A_s=30\text{ dB}\end{array}$$

Impulse invariance technique is to be used for converting an analog prototype filter to a discrete-time filter. Determine the specifications of the analog prototype if the sampling interval is to be

1. T = 1 s
2. T = 2 s

Solution:

1. Using T = 1 s, the critical frequencies of the analog prototype filter are

   $$\begin{array}{cc}{\omega }_1=\frac{{\Omega }_1}{T}=\frac{{\Omega }_1}{1}=0.2\pi ,& {\omega }_2=\frac{{\Omega }_2}{T}=\frac{{\Omega }_2}{1}=0.25\pi \end{array}$$

   The dB tolerance limits are unchanged:

   $$\begin{array}{cc}{R}_p=1\text{ dB},& A_s=30\text{ dB}\end{array}$$

   Let the system function for the analog prototype filter be G1 (ω) yielding an impulse response g1 (t). The impulse response of the discrete-time filter is

   $$h_1\left[n\right]=g_1\left(nT\right)=g_1\left(n\right)$$

2. With T = 2 s, the critical frequencies of the analog prototype filter are

   $$\begin{array}{cc}{\omega }_1=\frac{{\Omega }_1}{T}=\frac{{\Omega }_1}{2}=0.1\pi ,& {\omega }_2=\frac{{\Omega }_2}{T}=\frac{{\Omega }_2}{2}=0.125\pi \end{array}$$

   The dB tolerance limits are unchanged as before. Let the system function for the analog prototype filter be G2 (ω) so that its impulse response is g2 (t). The impulse response of the discrete-time filter is obtained by sampling g2 (t) every 2 seconds, that is,

   $$h_2\left[n\right]=g_2\left(nT\right)=g_2\left(2n\right)$$

We have thus obtained two discrete-time filters with the two choices of the sampling interval T. What is the relationship between these two filters? Let us realize that G1 (0.2π) = G2 (0.1π) and G1 (0.25π) = G2 (0.125π). Generalizing these relationships we have

$$G_1\left(\omega \right)=G_2\left(\frac{\omega }{2}\right)$$

Based on the scaling property of the Fourier transform (see Section 4.3.5 of Chapter 4) this implies that

$$g_1\left(n\right)=g_2\left(2n\right)$$

and therefore

$$h_1\left[n\right]=h_2\left[n\right]$$

The two IIR filters designed are identical and independent from the choice of T.

Bilinear transformation

The impulse invariance technique has a severe shortcoming that prevents its adoption as a general method applicable to the design of all filter types. The aliasing effect that results from sampling the impulse response limits its usefulness to lowpass filters only. An alternative technique known as bilinear transformation avoids aliasing and overcomes this limitation.

In order to develop the idea we will work with a first-order system; however, the technique developed will be applicable to higher-order systems as well. Consider a first-order analog prototype filter described by the system function

$$\begin{array}{cc}G\left(s\right)=\frac{Y_a\left(s\right)}{X_a\left(s\right)}=\frac{a}{s-b}& \left(10.146\right)\end{array}$$

The differential equation for the filter is

$$\begin{array}{cc}\frac{d{y}_a\left(t\right)}{dt}=b{y}_a\left(t\right)+a{x}_a\left(t\right)& \left(10.147\right)\end{array}$$

Let us define an intermediate variable ωa (t) as

$${\omega }_a\left(t\right)=\frac{d{y}_a\left(t\right)}{dt}$$

so that Eqn. (10.147) becomes

$$\begin{array}{cc}{\omega }_a\left(t\right)=b{y}_a\left(t\right)+a{x}_a\left(t\right)& \left(10.148\right)\end{array}$$

Since the goal is to obtain a discrete-time filter, we will write the discrete-time version of Eqn. (10.148) as

$$\begin{array}{cc}\omega \left[n\right]=by\left[n\right]+ax\left[n\right]& \left(10.149\right)\end{array}$$

where we have defined ω[n] = ωa (nT), y[n] = ya (nT) and x[n] = xa (nT). The equivalent relationship between ya(t) and ωa(t) is

$$\begin{array}{cc}{y}_a\left(t\right)={\int }_{t_0}^t{\omega }_a\left(t\right)dt+y_a\left(t_0\right)& \left(10.150\right)\end{array}$$

Letting t0 = (n − 1)T and t = nT, a discrete-time approximation to ya (t) can be obtained using the trapezoidal approximation method (see Problem 3.35 at the end of Chapter 3):

$$\begin{array}{cc}y\left[n\right]=\frac{T}{2}\left(\omega \left[n\right]+\omega \left[n-1\right]+y\left[n-1\right]\right)& \left(10.151\right)\end{array}$$

Substituting Eqn. (10.149) into Eqn. (10.151) the difference equation between x[n] and y[n] is found as

$$\begin{array}{cc}y\left[n\right]=\frac{T}{2}\left(by\left[n\right]+ax\left[n\right]+by\left[n-1\right]+ax\left[n-1\right]+y\left[n-1\right]\right)& \left(10.152\right)\end{array}$$

which can be simplified to

$$\begin{array}{cc}\left(1-b\frac{T}{2}\right)y\left[n\right]=\left(1+b\frac{T}{2}\right)y\left[n-1\right]+a\frac{T}{2}x\left[n\right]+a\frac{T}{2}x\left[n-1\right]& \left(10.153\right)\end{array}$$

The system function of the corresponding discrete-time filter is obtained from Eqn. (10.153) as

$$\begin{array}{cc}H\left(z\right)=\frac{a}{\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}-b}& \left(10.154\right)\end{array}$$

By comparing Eqns. (10.146) and (10.154) we conclude that the transformation between the s-plane and the z-plane is given by

$$\begin{array}{cc}s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}& \left(10.155\right)\end{array}$$

This is referred to as the bilinear transformation. Solving Eqn. (10.155) for z yields

$$\begin{array}{cc}z=\frac{2+sT}{2-sT}& \left(10.156\right)\end{array}$$

as the inverse of bilinear transformation. Since the relationship between s and z is invertible, bilinear transformation represents a one-to-one mapping between the s-plane and the z-plane. For each point in the s-plane, there is one and only one point in the z-plane, and vice versa. Contrast this with the impulse invariance technique that maps multiple frequencies on the jω-axis of the s-plane to a single point on the unit circle of the z-plane, resulting in aliasing.

In order to understand how the jω-axis of the s-plane is mapped to the z-plane we will evaluate Eqn. (10.156) for s = jω:

$$\begin{array}{cc}z=\frac{2+j\omega T}{2-j\omega T}& \left(10.157\right)\end{array}$$

The magnitude of the expression is

$$\begin{array}{cc}\left|z\right|=\frac{\sqrt{4+{\omega }^2{T}^2}}{\sqrt{4+{\omega }^2{T}^2}}=1& \left(10.158\right)\end{array}$$

Since the magnitude of z is equal to unity for all ω, the jω-axis of the s-plane is mapped onto the unit circle of the z-plane. Furthermore, the entire jω-axis of the s-plane is mapped onto the unit circle of the z-plane only once, and there is no aliasing. Bilinear transform does not suffer from the limitations of the impulse invariance technique, and can be used for designing all four frequency-selective filter types. Setting z = ejΩ in Eqn. (10.158) we have

$$\begin{array}{cc}{e}^{j\Omega }=\frac{\sqrt{4+{\omega }^2{T}^2}}{\sqrt{4+{\omega }^2{T}^2}}=1& \left(10.159\right)\end{array}$$

and

$$\begin{array}{cc}\Omega =\measuredangle e^{j\Omega }=2{tan}^{-1}\left(\frac{\omega T}{2}\right)& \left(10.160\right)\end{array}$$

Eqn. (10.160) gives the relationship between analog frequencies ω and discrete-time frequencies Ω. It is shown graphically in Fig. 10.31.

It is evident from Fig. 10.31 that the mapping of frequencies is highly nonlinear. For low frequencies the relationship between ω and Ω is almost linear, however, for high frequencies we see that significant changes in ω trigger minuscule changes in Ω. This behavior is necessary if we want to map an infinite range of analog frequencies to a 2π-radian range of discrete-time frequencies. Fig. 10.32 illustrates the mapping from the jω-axis of the s-plane and the unit circle of the z-plane.

The jω-axis of the s-plane is mapped onto the unit circle of the z-plane. Each point on the jω-axis of the s-plane maps to one and only one point on the unit circle of the z-plane. Furthermore, each point in the left half s-plane maps to one and only one point inside the unit circle of the z-plane. This ensures that a stable analog prototype filter leads to a stable discrete-time filter under bilinear transformation.

Example 10.7: Applying bilinear transformation

A second-order analog filter has the system function

$$G\left(s\right)=\frac{2}{\left(s+1\right)\left(s+2\right)}$$

Convert this filter to a discrete-time filter using bilinear transformation with T = 2 s. Afterwards compute and graph the magnitude response of the discrete-time filter.

Solution: Substituting

$$\begin{array}{cc}s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}=\frac{1-z^{-1}}{1+z^{-1}}& \left(10.161\right)\end{array}$$

we obtain

$$H\left(z\right)=\frac{2}{\left(\frac{1-z^{-1}}{1+z^{-1}}+1\right)\left(\frac{1-z^{-1}}{1+z^{-1}}+2\right)}$$

which can be simplified to

$$H\left(z\right)=\frac{z^2+2z+1}{z\left(3z+1\right)}$$

The magnitude and the phase of H (z) are shown in Fig. 10.33.

Software resources:

ex_10_7.m

Obtaining analog prototype specifications

As discussed earlier, the first step in designing IIR filters is to convert the specifications of the desired discrete-time filter to the corresponding specifications of an appropriate analog prototype filter. For lowpass and highpass filters two edge frequencies Ω1 and Ω2 are given along with the tolerance values Δ1 and Δ2, or their decibel equivalents Rp and As. The critical frequencies for the discrete-time IIR filter must be translated to the critical frequencies of the analog prototype based on which method is to be used in step 3 of the design process. If impulse invariance will be used in converting the analog prototype to a discrete-time filter, critical frequencies of the analog prototype are computed as

$$\begin{array}{cccc}{\omega }_1=\frac{{\Omega }_1}{T}& \text{and}& {\omega }_2=\frac{{\Omega }_2}{T}& \left(10.162\right)\end{array}$$

since the conversion process is based on sampling the impulse response of the analog prototype filter. If, on the other hand, bilinear transformation is to be used, then the distortion of the frequency axis must be taken into account. In this case the critical frequencies for the analog prototype filter are computed as

$$\begin{array}{cccc}{\omega }_1=\frac{2}{T}tan\left(\frac{{\Omega }_1}{2}\right)& \text{and}& {\omega }_2=\frac{2}{T}tan\left(\frac{{\Omega }_2}{2}\right)& \left(10.163\right)\end{array}$$

This is called prewarping. Recall that, when bilinear transformation is applied to the analog prototype filter in step 3 of the design process, the frequency axis is distorted as described by Eqn. (10.160). Prewarping counteracts this distortion at the critical frequencies.

Example 10.8: IIR filter design using bilinear transformation

Using bilinear transformation design a Butterworth lowpass filter with the following specifications:

$$\begin{array}{cccc}{\Omega }_1=0.2\pi ,& {\Omega }_2=0.36\pi ,& R_p=2\text{ dB},& A_s=20\text{ dB}\end{array}$$

Solution: We will use T = 1 s. The critical frequencies of the analog prototype filter are found using the prewarping equation applied to Ω1 and Ω2:

$$\begin{array}{c}{\omega }_1=2tan\left(0.2\pi /2\right)=0.6498\text{ rad}\\ {\omega }_2=2tan\left(0.36\pi /2\right)=1.2692\text{ rad}\end{array}$$

Next the filter order N needs to be determined using Eqn. (10.65).

$$N\ge \frac{{log}_{10}\sqrt{\left({10}^{2/10}-1\right)\left({10}^{20/10}-1\right)}}{{log}_{10}\left(0.6498/1.2692\right)}=3.8326$$

The filter order needs to be chosen as N = 4. The 3-dB cutoff frequency is found by setting the magnitude at the stopband edge to −As dB by solving

$${\left(\frac{{\omega }_2}{{\omega }_c}\right)}^{2N}={10}^{A_s/10}-1$$

so that

$${\left(\frac{1.2692}{{\omega }_c}\right)}^8={10}^{20/10}-1=99$$

which yields ωc = 0.7146 rad/s. The analog prototype filter can now be designed using Butterworth lowpass filter design technique described in Section 10.4.1 with the values of N and ωc found. The system function is

$$G\left(s\right)=\frac{0.2608}{s^4+1.8675s^3+1.7437s^2+0.9537s+0.2608}$$

The system function for the discrete-time filter is found through bilinear transformation using the replacement

$$s\to 2\frac{1-z^{-1}}{1+z^{-1}}$$

which results in

$$H\left(z\right)=\frac{0.0065z^4+0.0260z^3+0.0390z^2+0.0260z+0.0065}{z^4-2.2209z^3+2.0861z^2-0.9204z+0.1594}$$

The dB magnitude of the resulting filter is shown in Fig. 10.34 with the tolerance limits.

Software resources:

ex_10_8a.m

ex_10_8b.m

10.5.2 Design of FIR filters

A length-N FIR filter is completely characterized by its finite-length impulse response h[n] for n =0,...,N − 1. The system function for such a filter is computed as

$$\begin{array}{lll}H\left(z\right)\hfill & =\sum _{n=0}^{N-1}h\left[n\right]z^{-n}\hfill & \hfill \\ \hfill & =h\left[0\right]+h\left[1\right]z^{-1}+...+h\left[N-1\right]z^{-\left(N-1\right)}\hfill & \left(10.164\right)\hfill \end{array}$$

The system function can also be written using non-negative powers of z in the form

$$\begin{array}{ll}H\left(z\right)=\frac{h\left[0\right]z^{N-1}+h\left[1\right]z^{N-2}+...+h\left[N-2\right]z+h\left[N-1\right]}{z^{N-1}}\hfill & \left(10.165\right)\hfill \end{array}$$

which leads us to the following conclusions:

1. Length-N FIR filter has a system function that is of order N − 1.
2. The system function H (z) has N − 1 zeros and as many poles.
3. The placement of zeros of the filter is determined by the sample amplitudes of the impulse response h[n].
4. In contrast, all poles of the length-N FIR filter are at the origin, independent of the impulse response. Therefore, FIR filters are inherently stable.
5. The magnitude response of the filter is determined only by the locations of the zeros since the poles at the origin do not contribute to the magnitude response (see Section 8.5.4).

The sample amplitudes of the impulse response h[n] for n = 0,...,N − 1 are often referred to as the filter coefficients. In this section we will present two approaches to the problem of designing FIR filters: one very simplistic approach that is nevertheless instructive, and one elegant approach that is based on computer optimization. In a general sense, the design procedure for FIR filters consists of the following steps:

1. Start with the desired frequency response. Select the appropriate length N of the filter. This choice may be an educated guess, or may rely on empirical formulas.
2. Choose a design method that attempts to minimize, in some way, the difference between the desired frequency response and the actual frequency response that results.
3. Determine the filter coefficients h[n] using the design method chosen.
4. Compute the system function and decide if it is satisfactory. If not, repeat the process with a different value of N and/or a different design method.

It was discussed earlier in this chapter that one of the main reasons for preferring FIR filters over IIR filters is the possibility of a linear phase characteristic. Linear phase is desirable since it leads to a time-delay characteristic that is constant independent of frequency. As far as real-time implementations of IIR and FIR filters are concerned, we have seen that IIR filters are mathematically more efficient and computationally less demanding of hardware resources compared to FIR filters. If linear phase is not a significant concern in a particular application, then an IIR filter may be preferred. If linear phase is a requirement, on the other hand, an FIR filter must be chosen even though its implementation may be more costly. Therefore, in the discussion of FIR design we will focus on linear-phase FIR filters.

Conditions for linear phase

It can be shown that a length-N FIR filter with impulse response h[n] has a linear phase characteristic if

$$\begin{array}{llll}h\left[n\right]=\mp h\left[N-1-n\right]\hfill & \text{for}\hfill & n=0,...,K\hfill & \left(10.166\right)\hfill \end{array}$$

where the upper limit K is set as follows:

$$\begin{array}{ll}K=\{\begin{array}{cc}N/2-1,& N\text{ even}\\ \left(N-1\right)/2,& N\text{ odd}\end{array}\hfill & \left(10.167\right)\hfill \end{array}$$

The condition given by Eqns. (10.166) and (10.167) leads to four possible scenarios:

1. N : even, and h[n] = h[N − 1 − n] for N = 0,...,N/2 − 1
2. N : even, and h[n] = −h[N − 1 − n] for N = 0,...,N/2 − 1
3. N : odd, and h[n] = h[N − 1 − n] for N = 0,..., (N − 1)/2
4. N : odd, and h[n] = −h[N − 1 − n] for N = 0,..., (N − 1)/2

In addition, it can be shown that the constant time delay for a length-N linear-phase FIR filter is

$$\begin{array}{cc}{n}_d=\frac{N-1}{2}& \left(10.168\right)\end{array}$$

for all four possibilities. It is interesting to see that the time-delay nd is not an integer when N is even.

A complete proof of linear-phase conditions must address each of the four possibilities listed above (see Problems 10.17 and 10.19).

Example 10.9: Phase of FIR filter with symmetric impulse response

A length-5 FIR filter has the impulse response

$$h\left[n\right]=\{\underset{\stackrel{\uparrow }{n=0}}{3},2,1,2,3\}$$

Show that the phase characteristic of H (Ω) is a linear function of Ω.

Solution: The DTFT of the impulse response is

Let us factor out e−j2Ω and write the result as

$$\begin{array}{ll}H(\Omega )=\left[3e^{j2\Omega }+2e^{j\Omega }+1+2e^{-j\Omega }+3e^{-j2\Omega }\right]e^{-j2\Omega }\hfill & \left(10.169\right)\hfill \end{array}$$

The expression in square brackets contains symmetric exponentials. Using Euler’s formula, it becomes

$$\begin{array}{lll}A(\Omega )\hfill & =3e^{j2\Omega }+2e^{j\Omega }+1+2e^{-j\Omega }+3e^{-j2\Omega }\hfill & \hfill \\ \hfill & =1+cos(\Omega )+\frac{3}{2}cos\left(2\Omega \right)\hfill & \left(10.170\right)\hfill \end{array}$$

which is purely real. Using Eqn. (10.170) in Eqn. (10.169) we obtain

$$H(\Omega )=A(\Omega )e^{-j2\Omega }$$

The phase response of the filter is

$$\Theta (\Omega )=-j2\Omega$$

corresponding to a time delay of nd = 2 samples. The magnitude and the phase of the system function are shown in Fig. 10.35.

It should be noted that the phase characteristic is linear in spite of the vertical jumps in Fig. 10.35(b). These vertical sections are due to the fact that A (Ω) is negative for some frequencies. Therefore, the magnitude |H (Ω)| is

$$\left|H(\Omega )\right|=\left|A(\Omega )\right|$$

and the phase is

$$\measuredangle H(\Omega )=\{\begin{array}{cc}\Theta (\Omega ),& \text{if }A(\Omega )\ge 0\\ \Theta (\Omega )\mp \pi ,& \text{if }A(\Omega )<0\end{array}$$

What matters for linear phase is that all non-vertical sections of the phase characteristic have the same slope value.

Software resources:

ex_10_9.m

We are now ready to discuss FIR filter design. First the Fourier series design method will be discussed. The use of window functions to overcome the fundamental problems encountered will be explored. Afterwards we will briefly discuss the Parks McClellan technique.

Fourier series design of FIR filters

An ideal lowpass discrete-time filter with a cutoff frequency Ωc has the system function

$$\begin{array}{ll}{H}_d(\Omega )=\{\begin{array}{ll}1,\hfill & |\Omega |<{\Omega }_c\hfill \\ 0,\hfill & {\Omega }_c<|\Omega |<\pi \hfill \end{array}\hfill & \left(10.171\right)\hfill \end{array}$$

as shown in Fig. 10.36. The filter has zero phase and zero time delay.

The impulse response of the ideal discrete-time lowpass filter is found by taking the inverse DTFT of Hd (Ω) given by Eqn. (10.171).

$$\begin{array}{lll}{h}_d\left[n\right]\hfill & =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{H}_d(\Omega )e^{j\Omega n}d\Omega \hfill & \hfill \\ \hfill & =\frac{1}{2\pi }{\int }_{-{\Omega }_c}^{{\Omega }_c}\left(1\right)e^{j\Omega n}d\Omega =\frac{{\Omega }_c}{\pi }\text{sinc}\left(\frac{{\Omega }_{c}n}{\pi }\right)\hfill & \left(10.172\right)\hfill \end{array}$$

The result in Eqn. (10.172) is shown in Fig. 10.37 for Ωc = 0.4π.

The result obtained in Eqn. (10.172) for hd[n] is valid for all N. Therefore, hd[n] infinitely long, and cannot be the impulse response of an FIR filter. On the other hand, sample amplitudes of hd[n] get smaller as the sample index is increased in both directions. A finite-length impulse response can be obtained by truncating hd[n] as follows:

$$\begin{array}{ll}{h}_T\left[n\right]=\{\begin{array}{ll}{h}_d\left[n\right],\hfill & -M\ge n\ge M\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\hfill & \left(10.173\right)\hfill \end{array}$$

The truncated impulse response has 2M + 1 samples for n = −M,...,M. Truncation of the ideal impulse response causes the spectrum of the filter to deviate from the ideal spectrum. The system function for the resulting filter is

$$\begin{array}{ll}{H}_T(\Omega )={\int }_{n=-\infty }^{\infty }{h}_T\left[n\right]e^{-j\Omega n}={\int }_{n=-M}^M{h}_d\left[n\right]e^{-j\Omega n}\hfill & \left(10.174\right)\hfill \end{array}$$

An example is shown in Fig. 10.38. The truncated impulse response with M = 11 is shown in Fig. 10.38(a), and the magnitude spectrum of the resulting length-23 filter is shown in Fig. 10.38(b).

The truncated impulse response hT[n] is still non-causal owing to the fact that the filter has zero phase. In order to obtain a causal filter, a delay of M samples must be incorporated into the impulse response, resulting in

$$\begin{array}{lll}h\left[n\right]=h_T\left[n-M\right],\hfill & n=0,...,2M\hfill & \left(10.175\right)\hfill \end{array}$$

Thus, h[n] corresponds to a causal FIR filter. Using the time shifting property of the DTFT (see Section 5.3.5 of Chapter 5), the system function H (Ω) of the causal filter is related to HT (Ω) by

$$\begin{array}{ll}H(\Omega )=H_T(\Omega )e^{-j\Omega M}\hfill & \left(10.176\right)\hfill \end{array}$$

The addition of the M-sample delay only affects the phase of the system function and not its magnitude. Since h[n] is both causal and finite-length, it is the impulse response of a valid FIR filter.

It is worth observing that the impulse response hd[n] given by Eqn. (10.172) is an even function of n. Even symmetry is preserved in hT[n] obtained in Eqn. (10.173) by truncating hd[n]. Finally, when h[n] is obtained by time shifting hT[n] by M samples, the symmetry necessary for linear phase is preserved. Therefore, any filter designed using this technique will have linear phase.

Example 10.10: Fourier series design

Using the Fourier series method, design a length-15 FIR lowpass filter to approximate an ideal lowpass filter with Ωc = 0.3π rad.

Solution: The impulse response of the ideal lowpass filter is

$$h_d\left[n\right]=\frac{0.3\pi }{\pi }\text{sinc}\left(\frac{0.3\pi n}{\pi }\right)=0.3\text{ sinc}\left(0.3n\right)$$

Since N = 2M + 1 = 15 we have M = 7. The truncated impulse response is

$$h_T\left[n\right]=\{\begin{array}{ll}0.3\text{ sinc}\left(0.3n\right),\hfill & -7\le n\le 7\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}$$

The impulse response of the FIR filter is

$$\begin{array}{ll}h\left[n\right]=h_T\left[n-7\right]=0.3\text{ sinc}\left(0.3\left(n-7\right)\right),\hfill & n=0,...,14\hfill \end{array}$$

which can be evaluated as

$$h\left[n\right]=\{\underset{\stackrel{\uparrow }{n=0}}{0.014},-0.031,-0.064,-0.047,0.033,0.151,0.258,0.3,0.258,0.151,0.032,-0.047,-0.064,-0.031,-0.014\}$$

The magnitude response of the designed filter is shown in Fig. 10.39.

Software resources:

ex_10_10a.m

ex_10_10b.m

A by-product of truncating the impulse response in Eqn. (10.173) is the oscillatory behavior of the frequency response that is particularly evident around the cutoff frequency Ωc. This effect was observed as the Gibbs phenomenon in earlier chapters. An alternative way of expressing the truncation of hd[n] given by Eqn. (10.173) is

$$\begin{array}{ll}{h}_T\left[n\right]=h_d\left[n\right]\omega \left[n\right]\hfill & \left(10.177\right)\hfill \end{array}$$

where ω[n] is a rectangular window function defined as

$$\begin{array}{ll}\omega \left[n\right]=\{\begin{array}{ll}1,\hfill & -M\le n\le M\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\hfill & \left(10.178\right)\hfill \end{array}$$

Based on the multiplication property of the DTFT (see Section 5.3.5 of Chapter 5), the spectrum of hT[n] is the convolution of the spectra of hd[n] and ω[n].

$$\begin{array}{ll}{H}_T(\Omega )=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{H}_d\left(\lambda \right)W\left(\Omega -\lambda \right)d\lambda \hfill & \left(10.179\right)\hfill \end{array}$$

which is illustrated in Fig. 10.40.

As the spectra in Fig. 10.40 reveal, the reason behind the oscillatory behavior of HT (Ω) especially around Ω = Ωc is the shape of the spectrum W (Ω). The high-frequency content of W (Ω), on the other hand, is mainly due to the abrupt transition of the rectangular window from unit amplitude to zero amplitude at its two edges. The solution is to use an alternative window function in Eqn. (10.177), one that smoothly tapers down to zero at its edges, in place of the rectangular window. The chosen window function must also be an even function of n in order to keep the symmetry of hT[n] for linear phase. A large number of window functions exist in the literature. A few of them are listed below:

Triangular (Bartlett) window:

$$\begin{array}{lll}\omega \left[n\right]=1-\frac{\left|n\right|}{M},\hfill & -M\le n\le M\hfill & \left(10.180\right)\hfill \end{array}$$

Hamming window:

$$\begin{array}{lll}\omega \left[n\right]=0.54-0.46cos\left(\frac{\pi \left(n+M\right)}{M}\right),\hfill & -M\le n\le M\hfill & \left(10.181\right)\hfill \end{array}$$

Hanning window:

$$\begin{array}{lll}\omega \left[n\right]=0.5-0.5cos\left(\frac{\pi \left(n+M\right)}{M}\right),\hfill & -M\le n\le M\hfill & \left(10.182\right)\hfill \end{array}$$

Blackman window:

$$\begin{array}{lll}\omega \left[n\right]=0.42-0.5cos\left(\frac{\pi \left(n+M\right)}{M}\right)+0.08cos\left(\frac{2\pi \left(n+M\right)}{M}\right),\hfill & -M\le n<M\hfill & \left(10.183\right)\hfill \end{array}$$

The four window functions listed are shown in Fig. 10.41 for M = 12.

Fig 10.42 shows the decibel magnitude spectra for rectangular, Hamming and Blackman windows for comparison.

It is seen that the rectangular window has stronger side lobes than the other two window functions. Hamming window provides side lobes that are at least 40 dB below the main lobe, and the side lobes for the Blackman window are at least 60 dB below the main lobe. The downside to side lobe suppression is the widening of the main lobe.

We are now ready to summarize the Fourier series design method for linear-phase FIR filters:

Example 10.11: Fourier series design using window functions

Redesign the filter of Example 10.10 using Hamming and Blackman windows.

Solution: Using a window function the truncated impulse response is

$$\begin{array}{ll}{h}_T\left[n\right]=0.3\text{ sinc}\left(0.3n\right)\omega \left[n\right],\hfill & -7\le n\le 7\hfill \end{array}$$

and the impulse response of the FIR filter is

$$\begin{array}{lll}h\left[n\right]=h_T\left[n-7\right]=0.3\text{ sinc}\left(0.3\left(n-7\right)\right)\omega \left[n-7\right],\hfill & n=0,...,14\hfill & \left(10.184\right)\hfill \end{array}$$

For the Hamming window, Eqn. (10.184) yields

$$h\left[n\right]=\{\underset{\stackrel{\uparrow }{n=0}}{0.006},-0.017-0.042,-0.036,0.028,0.142,0.254,0.3,0.254,0.142,0.028,-0.036,-0.042,-0.017,0.006\}$$

When the Blackman window is used, we get

$$h\left[n\right]=\{\underset{\stackrel{\uparrow }{n=0}}{0.003},-0.011-0.031,-0.030,0.025,0.135,0.250,0.3,0.250,0.135,0.025,-0.030,-0.031,-0.011,0.003\}$$

Magnitude responses of the two filters are shown in Fig. 10.43.

Software resources:

ex_10_11a.m

ex_10_11b.m

Filter types other than lowpass

In the discussion above we have concentrated on the design of lowpass FIR filters. Other types of filters can also be designed; the only modification that is needed is in determining the ideal impulse response hd[n]. Expressions for the ideal impulse responses of highpass, bandpass and band-reject filters can be derived from Eqn. (10.172) as follows:

Highpass:

$$\begin{array}{ll}{H}_d(\Omega )=\{\begin{array}{ll}1,\hfill & {\Omega }_c<|\Omega |<\pi \hfill \\ 0,\hfill & |\Omega |<{\Omega }_c\hfill \end{array}\hfill & h_d\left[n\right]=\delta \left[n\right]-\frac{{\Omega }_c}{\pi }\text{sinc}\left(\frac{{\Omega }_{c}n}{\pi }\right)\hfill \end{array}$$

Bandpass:

$$\begin{array}{ll}{H}_d(\Omega )=\{\begin{array}{ll}1,\hfill & {\Omega }_1<|\Omega |<{\Omega }_2\hfill \\ 0,\hfill & |\Omega |<{\Omega }_1\hfill \\ \hfill & \text{or }{\Omega }_2<|\Omega |<\pi \hfill \end{array}\hfill & h_d\left[n\right]=\frac{{\Omega }_2}{\pi }\text{sinc}\left(\frac{{\Omega }_{2}n}{\pi }\right)-\frac{{\Omega }_1}{\pi }\text{sinc}\left(\frac{{\Omega }_{1}n}{\pi }\right)\hfill \end{array}$$

Band-reject:

$$\begin{array}{ll}{H}_d(\Omega )=\{\begin{array}{ll}1,\hfill & |\Omega |<{\Omega }_1\hfill \\ \hfill & \text{or }{\Omega }_2<|\Omega |<\pi \hfill \\ 0,\hfill & {\Omega }_1<|\Omega |<{\Omega }_2\hfill \end{array}\hfill & h_d\left[n\right]=\delta \left[n\right]-\frac{{\Omega }_2}{\pi }\text{sinc}\left(\frac{{\Omega }_{2}n}{\pi }\right)+\frac{{\Omega }_1}{\pi }\text{sinc}\left(\frac{{\Omega }_{1}n}{\pi }\right)\hfill \end{array}$$

Parks-McClellan technique for FIR filter design

FIR filter design technique due to Parks and McClellan is also referred to as the Remez exchange method. It is based on the idea of finding a polynomial approximation to the desired frequency response by minimizing the largest approximation error.

The Fourier series design method discussed in the previous section, when used with a rectangular window, provides the optimum solution in the sense of mean-squared error. If the error between the desired frequency response and the actual frequency response is E (Ω), then the length-N FIR filter designed using the Fourier series method leads to the mean-squared error

$$\begin{array}{ll}\text{MSE}=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{\left|E(\Omega )\right|}^{2}d\Omega \hfill & \left(10.185\right)\hfill \end{array}$$

that is the smallest possible with N coefficients. On the other hand, we have seen that the magnitude response of the filter suffers from oscillatory behavior near discontinuities of the spectrum.

Better approximations can be obtained if we minimize the maximum value of the error rather than its mean-squared value. Parks-McClellan algorithm begins with an initial guess of the set of frequencies for the extrema of E (Ω). The locations of the extrema are then iteratively shifted until the maximum deviation from the ideal frequency response is made as small as possible.

Computational details of the Parks-McClellan algorithm are outside the scope of this text and may be found in references. We will provide examples of using MATLAB for designing Parks-McClellan filters in the MATLAB Exercises section.

10.6 Further Reading

[1] A.V. Oppenheim and R.W. Schafer. Discrete-Time Signal Processing. Prentice Hall, 2010.

[2] R. Schaumann and M.E. Van Valkenburg. Design of Analog Filters. Oxford University Press, 2010.

[3] R.A. Schilling, R.J. Schilling, and P.D. Sandra L. Harris. Fundamentals of Digital Signal Processing Using MATLAB. Cengage Learning, 2010.

[4] L. Tan and J. Jiang. Digital Signal Processing: Fundamentals and Applications. Elsevier Science, 2013.

[5] F. Taylor. Digital Filters: Principles and Applications with MATLAB. IEEE Series on Digital and Mobile Communication. Wiley, 2011.

[6] L.D. Thede. Practical Analog and Digital Filter Design. Artech House, 2005.

[7] L. Wanhammar. Analog Filters Using MATLAB. Springer London, Limited, 2009.

MATLAB Exercises