[4.12 MATLAB: CONTINUOUS-TIME](#page-11-0) FILTERS

← Back to LINEAR SYSTEMS AND SIGNALS Overview

4.12 MATLAB: CONTINUOUS-TIME FILTERS

Continuous-time filters are essential to many if not most engineering systems, and MATLAB is an excellent assistant for filter design and analysis. Although a comprehensive treatment of continuous-time filter techniques is outside the scope of this book, quality filters can be designed and realized with minimal additional theory.

A simple yet practical example demonstrates basic filtering concepts. Telephone voice signals are often lowpass-filtered to eliminate frequencies above a cutoff of 3 kHz, or ω*^c* = 3000(2π ) ≈ 18,850 rad/s. Filtering maintains satisfactory speech quality and reduces signal bandwidth, thereby increasing the phone company’s call capacity. How, then, do we design and realize an acceptable 3 kHz lowpass filter?

4.12-1 Frequency Response and Polynomial Evaluation

Magnitude response plots help assess a filter’s performance and quality. The magnitude response of an ideal filter is a brick-wall function with unity passband gain and perfect stopband attenuation. For a lowpass filter with cutoff frequency ωc, the ideal magnitude response is

$$|H_{\text{ideal}}(j\omega)| = \begin{cases} 1 & |\omega| \le \omega_c \\ 0 & |\omega| > \omega_c \end{cases}$$

Unfortunately, ideal filters cannot be implemented in practice. Realizable filters require compromises, although good designs will closely approximate the desired brick-wall response.

A realizable LTIC system often has a rational transfer function that is represented in the s-domain as

$$H(s) = \frac{Y(s)}{X(s)} = \frac{B(s)}{A(s)} = \frac{\sum_{k=0}^{M} b_{k+N-M} s^{M-k}}{\sum_{k=0}^{N} a_k s^{N-k}}$$

Frequency response H(jω) is obtained by letting s = jω, where frequency ω is in radians per second.

MATLAB is ideally suited to evaluate frequency response functions. Defining a length-(N + 1) coefficient vector A = [a0,a1,…,aN] and a length-(M + 1) coefficient vector B = [bN−^M,bN−M+1,…, bN], program CH4MP1 computes H(jω) for each frequency in the input vector ω.

function [H] = CH4MP1(B,A,omega);
% CH4MP1.m : Chapter 4, MATLAB Program 1
% Function M-file computes frequency response for LTIC system
% INPUTS: B = vector of feedforward coefficients
% A = vector of feedback coefficients
% omega = vector of frequencies [rad/s].
% OUTPUTS: H = frequency response

H = polyval(B,j*omega)./polyval(A,j*omega);

The function polyval efficiently evaluates simple polynomials and makes the program nearly trivial. For example, when A is the vector of coefficients [a0,a1,…,aN], polyval (A,j*omega) computes

$$\sum_{k=0}^N a_k (j\omega)^{N-k}$$

for each value of the frequency vector omega. It is also possible to compute frequency responses by using the signal-processing toolbox function freqs.

DESIGN AND EVALUATION OF A SIMPLE RC FILTER

One of the simplest lowpass filters is realized by using an RC circuit, as shown in Fig. 4.61. This one-pole system has transfer function HRC(s) = (RCs + 1)−¹ and magnitude response |HRC(jω)|=|(jωRC + 1)−¹| = 1/ 1+(RCω)2. Independent of component values R and C, this circuit has many desirable characteristics, such as unity gain at ω = 0 and magnitude response that monotonically decreases to zero as ω → ∞.

Components R and C are chosen to set the desired 3 kHz cutoff frequency. For many filter types, the cutoff frequency corresponds to the half-power point, or |HRC(jωc)| = 1/ ^√2. Assign C a realistic capacitance of 1 nF, then the required resistance is computed by R = 1/ C²ω² c = 1/ (10−⁹)²(2π3000)2.

>> omega_c = 2*pi*3000; C = 1e-9; R = 1/sqrt(C^2*omega_c^2)
    R = 5.3052e+004

The root of this first-order RC filter is directly related to the cutoff frequency, λ = −1/RC = −18,850 = −ωc.

To evaluate the RC filter performance, the magnitude response is plotted over the mostly audible frequency range (0 ≤ f ≤ 20 kHz).

>> f = linspace(0,20000,200); Hmag_RC = abs(CH4MP1([1],[R*C 1],f*2*pi));

>> plot(f,abs(f*2*pi)<=omega_c,'k-',f,Hmag_RC,'k--');

>> axis([0 20000 -0.05 1.05]); xlabel('f [Hz]'); ylabel('|H(j2\pi f)|');

>> legend('Ideal','First-order RC','location','best');

The linspace(X1,X2,N) command generates an N-length vector of linearly spaced points between X1 and X2.

As shown in Fig. 4.62, the first-order RC response is indeed lowpass with a half-power cutoff frequency equal to 3 kHz. It rather poorly approximates the desired brick-wall response: the passband is not very flat, and stopband attenuation increases very slowly to less than 20 dB at 20 kHz.

Figure 4.62 Magnitude response |HRC(j2πf)| of a first-order RC filter.

Figure 4.63 A cascaded RC filter.

A CASCADED RC FILTER AND POLYNOMIAL EXPANSION

A first-order RC filter is destined for poor performance; one pole is simply insufficient to obtain good results. A cascade of RC circuits increases the number of poles and improves the filter response. To simplify the analysis and prevent loading between stages, we employ op-amp followers to buffer the output of each stage, as shown in Fig. 4.63. A cascade of N stages results in an Nth-order filter with transfer function given by

$$H_{\text{cascade}}(s) = [H_{RC}(s)]^N = (RCs + 1)^{-N}$$

Upon choosing a cascade of 10 stages and C = 1 nF, a 3 kHz cutoff frequency is obtained by setting ^{R = √ 21}/^{10 −}1/(Cωc) ^{= √} 21/¹⁰ −1/(6π(10)−⁶).

R = sqrt(2^(1/10)-1)/(C*omega_c) R = 1.4213e+004

This cascaded filter has a 10th-order pole at λ = −1/RC and no finite zeros. To compute the magnitude response, polynomial coefficient vectors A and B are needed. Setting B = [1] ensures there are no finite zeros or, equivalently, that all zeros are at infinity. The poly command, which expands a vector of roots into a corresponding vector of polynomial coefficients, is used to obtain A.

• >> B = 1; A = poly(-1/(R*C)*ones(10,1));A = A/A(end);

• Hmag_cascade = abs(CH4MP1(B,A,f*2*pi));

• plot(f,abs(f*2*pi)<=omega_c,‘k-’,f,Hmag_cascade,‘k—’);

• axis([0 20000 -0.05 1.05]); xlabel(‘f [Hz]’); ylabel(‘|H(j2\pi f)|’);

• legend(‘Ideal’,‘Tenth-order RC cascade’,‘location’,‘best’);

Notice that scaling a polynomial by a constant does not change its roots. Conversely, the roots of a polynomial specify a polynomial within a scale factor. The command A = A/A(end) properly scales the denominator polynomial to ensure unity gain at ω = 0.

The magnitude response plot of the tenth-order RC cascade is shown in Fig. 4.64. Compared with the simple RC response of Fig. 4.62, the passband remains relatively unchanged, but stopband attenuation is greatly improved to over 60 dB at 20 kHz.

Figure 4.64 Magnitude response |Hcascade(j2πf)| of a tenth-order RC cascade.

4.12-2 Butterworth Filters and the Find Command

The pole location of a first-order lowpass filter is necessarily fixed by the cutoff frequency. There is little reason, however, to place all the poles of a 10th-order filter at one location. Better pole placement will improve our filter’s magnitude response. One strategy, discussed in Sec. 4.10, is to place a wall of poles opposite the passband frequencies. A semicircular wall of poles leads to the Butterworth family of filters, and a semi-elliptical shape leads to the Chebyshev family of filters. Butterworth filters are considered first.

To begin, notice that a transfer function H(s) with real coefficients has a squared magnitude response given by |H(jω)| ^{2 = H}(jω)H∗(jω) ^{= H}(jω)H(−jω) ^{= H}(s)H(−s)|^s=jω. Thus, half the poles of |H(jω)| ² correspond to the filter H(s) and the other half correspond to H(−s). Filters that are both stable and causal require H(s) to include only left-half-plane poles.

The squared magnitude response of a Butterworth filter is

$$|H_{\text{BW}}(j\omega)|^2 = \frac{1}{1 + (j\omega/j\omega_c)^{2N}}$$

This function has the same appealing characteristics as the first-order RC filter: a gain that is unity at ω = 0 and monotonically decreases to zero as ω → ∞. By construction, the half-power gain

occurs at ωc. Perhaps most importantly, however, the first 2N − 1 derivatives of |HBW(jω)| with respect to ω are zero at ω = 0. Put another way, the passband is constrained to be very flat for low frequencies. For this reason, Butterworth filters are sometimes called maximally flat filters.

As discussed in Sec. B.7, the roots of minus 1 must lie equally spaced on a circle centered at the origin. Thus, the 2N poles of |HBW(jω)| ² naturally lie equally spaced on a circle of radius ω*^c* centered at the origin. Figure 4.65 displays the 20 poles corresponding to the case N = 10 and ω*^c* = 3000(2π ) rad/s. An Nth-order Butterworth filter that is both causal and stable uses the N left-half-plane poles of |HBW(jω)| 2.

To design a 10th-order Butterworth filter, we first compute the 20 poles of |HBW(jω)| 2:

N=10; poles = roots([(1j*omega_c)^(-2*N),zeros(1,2*N-1),1]);

The find command is a powerful and useful function that returns the indices of a vector’s nonzero elements. Combined with relational operators, the find command allows us to extract the 10 left-half-plane roots that correspond to the poles of our Butterworth filter.

BW_poles = poles(find(real(poles)<0));

To compute the magnitude response, these roots are converted to coefficient vector A.

>> A = poly(BW_poles); A = A/A(end); Hmag_BW = abs(CH4MP1(B,A,f*2*pi));

>> plot(f,abs(f*2*pi)<=omega_c,'k-',f,Hmag_BW,'k--');

>> axis([0 20000 -0.05 1.05]); xlabel('f [Hz]'); ylabel('|H(j2\pi f)|');

>> legend('Ideal','Tenth-order Butterworth','location','best');

The magnitude response plot of the Butterworth filter is shown in Fig. 4.66. The Butterworth response closely approximates the brick-wall function and provides excellent filter characteristics: flat passband, rapid transition to the stopband, and excellent stopband attenuation (>40 dB at 5 kHz).

Figure 4.66 Magnitude response |HBW(j2πf)| of a tenth-order Butterworth filter.

4.12-3 Using Cascaded Second-Order Sections for Butterworth Filter Realization

For our RC filters, realization preceded design. For our Butterworth filter, however, design has preceded realization. For our Butterworth filter to be useful, we must be able to implement it.

Since the transfer function HBW(s) is known, the differential equation is also known. Therefore, it is possible to try to implement the design by using op-amp integrators, summers, and scalar multipliers. Unfortunately, this approach will not work well. To understand why, consider the denominator coefficients a⁰ = 1.766×10−⁴³ and a¹⁰ = 1. The smallest coefficient is 43 orders of magnitude smaller than the largest coefficient! It is practically impossible to accurately realize such a broad range in scale values. To understand this, skeptics should try to find realistic resistors such that Rf /R = 1.766×10−43. Additionally, small component variations will cause large changes in actual pole location.

A better approach is to cascade five second-order sections, where each section implements one complex conjugate pair of poles. By pairing poles in complex conjugate pairs, each of the resulting second-order sections has real coefficients. With this approach, the smallest coefficients are only about nine orders of magnitude smaller than the largest coefficients. Furthermore, pole placement is typically less sensitive to component variations for cascaded structures.

The Sallen–Key circuit shown in Fig. 4.67 provides a good way to realize a pair of complex-conjugate poles.† The transfer function of this circuit is

$$H_{SK}(s) = \frac{\frac{1}{R_1R_2C_1C_2}}{s^2 + \left(\frac{1}{R_1C_1} + \frac{1}{R_2C_1}\right)s + \frac{1}{R_1R_2C_1C_2}} = \frac{\omega_0^2}{s^2 + \left(\frac{\omega_0}{Q}\right)s + \omega_0^2}$$

Geometrically, ω⁰ is the distance from the origin to the poles and Q = 1/2 cosψ, where ψ is the angle between the negative real axis and the pole. Termed the “quality factor” of a circuit, Q

^† A more general version of the Sallen–Key circuit has a resistor Ra from the negative terminal to ground and a resistor Rb between the negative terminal and the output. In Fig. 4.67, Ra = ∞ and Rb = 0.

provides a measure of the peakedness of the response. High-Q filters have poles close to the ω axis, which boost the magnitude response near those frequencies.

Although many ways exist to determine suitable component values, a simple method is to assign R¹ a realistic value and then let R² = R1, C¹ = 2Q/ω0R1, and C² = 1/2Qω0R2. Butterworth poles are a distance ω*^c* from the origin, so ω⁰ = ωc. For our 10th-order Butterworth filter, the angles ψ are regularly spaced at 9, 27, 45, 63, and 81 degrees. MATLAB program CH4MP2 automates the task of computing component values and magnitude responses for each stage.

% CH4MP2.m : Chapter 4, MATLAB Program 2
% Script M-file computes Sallen-Key component values and magnitude
% responses for each of the five cascaded second-order filter sections.
omega_0 = 3000*2*pi; % Filter cut-off frequency
psi = [9 27 45 63 81]*pi/180; % Butterworth pole angles
f = linspace(0,6000,200); % Frequency range for magnitude response calculations
Hmag_SK = zeros(5,200); % Pre-allocate array for magnitude responses
for stage = 1:5,
    Q = 1/(2*cos(psi(stage))); % Compute Q for current stage
    % Compute and display filter components to the screen:
    disp(['Stage ',num2str(stage),...
            ' (Q = ',num2str(Q),...
            '): R1 = R2 = ',num2str(56000),...
            ', C1 = ',num2str(2*Q/(omega_0*56000)),...
            ', C2 = ',num2str(1/(2*Q*omega_0*56000))]);
    B = omega_0^2; A = [1 omega_0/Q omega_0^2]; % Compute filter coefficients
    Hmag_SK(stage,:) = abs(CH4MP1(B,A,2*pi*f)); % Compute magnitude response
end
plot(f,Hmag_SK,'k',f,prod(Hmag_SK),'k:')
xlabel('f [Hz]'); ylabel('Magnitude Response')

The disp command displays a character string to the screen. Character strings must be enclosed in single quotation marks. The num2str command converts numbers to character strings and facilitates the formatted display of information. The prod command multiplies along the columns of a matrix; it computes the total magnitude response as the product of the magnitude responses of the five stages.

Executing the program produces the following output:

>> CH4MP2
   Stage 1 (Q = 0.50623): R1 = R2 = 56000, C1 = 9.5916e-10, C2 = 9.3569e-10
   Stage 2 (Q = 0.56116): R1 = R2 = 56000, C1 = 1.0632e-09, C2 = 8.441e-10
   Stage 3 (Q = 0.70711): R1 = R2 = 56000, C1 = 1.3398e-09, C2 = 6.6988e-10
   Stage 4 (Q = 1.1013): R1 = R2 = 56000, C1 = 2.0867e-09, C2 = 4.3009e-10
   Stage 5 (Q = 3.1962): R1 = R2 = 56000, C1 = 6.0559e-09, C2 = 1.482e-10

Figure 4.68 Magnitude responses for Sallen–Key filter stages.

Since all the component values are practical, this filter is possible to implement. Figure 4.68 displays the magnitude responses for all five stages (solid lines). The total response (dotted line) confirms a 10th-order Butterworth response. Stage 5, which has the largest Q and implements the pair of conjugate poles nearest the ω axis, is the most peaked response. Stage 1, which has the smallest Q and implements the pair of conjugate poles furthest from the ω axis, is the least peaked response. In practice, it is best to order high-Q stages last; this reduces the risk that the high gains will saturate the filter hardware.

4.12-4 Chebyshev Filters

Like an order-N Butterworth lowpass filter (LPF), an order-N Chebyshev LPF is an all-pole filter that possesses many desirable characteristics. Compared with an equal-order Butterworth filter, the Chebyshev filter achieves better stopband attenuation and reduced transition bandwidth by allowing an adjustable amount of ripple within the passband.

The squared magnitude response of a Chebyshev filter is

$$|H_{\rm C}(j\omega)|^2 = \frac{1}{1 + \epsilon^2 C_N^2(\omega/\omega_c)}$$

where controls the passband ripple, CN(ω/ωc) is a degree-N Chebyshev polynomial, and ω*^c* is the radian cutoff frequency. Several characteristics of Chebyshev LPFs are noteworthy:

• An order-N Chebyshev LPF is equi-ripple in the passband (|ω| ≤ ωc), has a total of N maxima and minima over (0 ≤ ω ≤ ωc), and is monotonic decreasing in the stopband (|ω| > ωc).

464 CHAPTER 4 CONTINUOUS-TIME SYSTEM ANALYSIS

• In the passband, the maximum gain is 1 and the minimum gain is 1/ ^√1⁺ 2. For odd-valued N, |H(j0)| = 1. For even-valued N, |HC(j0)| = 1/ ^√1⁺ 2.
• Ripple is controlled by setting ^{= √} 10*^R*/¹⁰ −1, where R is the allowable passband ripple expressed in decibels. Reducing adversely affects filter performance (see Prob. 4.12-10).
• Unlike Butterworth filters, the cutoff frequency ω*^c* rarely specifies the 3 dB point. For = 1, |HC(jωc)| ² = 1/(1+ ²) = 0.5. The cutoff frequency ω*^c* simply indicates the frequency after which |HC(jω)| < 1/ √ 1+ 2.

The Chebyshev polynomial CN(x) is defined as

$$C_N(x) = \cos[N\cos^{-1}(x)] = \cosh[N\cosh^{-1}(x)]$$

In this form, it is difficult to verify that CN(x) is a degree-N polynomial in x. A recursive form of CN(x) makes this fact more clear (see Prob. 4.12-13).

$$C_N(x) = 2xC_{N-1}(x) - C_{N-2}(x)$$

With C0(x) = 1 and C1(x) = x, the recursive form shows that any CN is a linear combination of degree-N polynomials and is therefore a degree-N polynomial itself. For N ≥ 2, MATLAB program CH4MP3 generates the (N +1) coefficients of Chebyshev polynomial CN(x).

function [C_N] = CH4MP3(N);
% CH4MP3.m : Chapter 4, MATLAB Program 3
% Function M-file computes Chebyshev polynomial coefficients
% using the recursion relation C_N(x) = 2xC_{N-1}(x) - C_{N-2}(x)
% INPUTS: N = degree of Chebyshev polynomial
% OUTPUTS: C_N = vector of Chebyshev polynomial coefficients
C_Nm2 = 1; C_Nm1 = [1 0]; % Initial polynomial coefficients:
for t = 2:N;
   C_N = 2*conv([1 0],C_Nm1)-[zeros(1,length(C_Nm1)-length(C_Nm2)+1),C_Nm2];
   C_Nm2 = C_Nm1; C_Nm1 = C_N;

end

As examples, consider C2(x) = 2xC1(x) − C0(x) = 2x(x) − 1 = 2x2 − 1 and C3(x) = 2xC2(x) −
C1(x) = 2x(2x2 −1)−x = 4x3 −3x. CH4MP3 easily confirms these cases.

>> CH4MP3(2)
  ans = 2 0 -1
>> CH4MP3(3)
  ans = 4 0 -3 0

Since CN(ω/ωc) is a degree-N polynomial, |HC(jω)| ² is an all-pole rational function with 2N finite poles. Similar to the Butterworth case, the N poles specifying a causal and stable Chebyshev filter can be found by selecting the N left-half-plane roots of 1+ ²C^{2 N}[s/(jωc)].

Root locations and dc gain are sufficient to specify a Chebyshev filter for a given N and . To demonstrate, consider the design of an order-8 Chebyshev filter with cutoff frequency fc = 1 kHz and allowable passband ripple R = 1 dB. First, filter parameters are specified.

>> omega_c = 2*pi*1000; R = 1; N = 8;
>> epsilon = sqrt(10^(R/10)-1);

The coefficients of CN[s/(jωc)] are obtained with the help of CH4MP3, and then the coefficients of [1+ ²C^{2 N}(s/(jωc))] are computed by using convolution to perform polynomial multiplication.

>> CN = CH4MP3(N).*((1/(1j*omega_c)).^[N:-1:0]);
>> CP = epsilon^2*conv(CN,CN); CP(end) = CP(end)+1;

Next, the polynomial roots are found, and the left-half-plane poles are retained and plotted.

>> poles = roots(CP); i = find(real(poles)<0); C_poles = poles(i);
>> plot(real(C_poles),imag(C_poles),'kx'); axis equal;
>> axis(omega_c*[-1.1 1.1 -1.1 1.1]);
>> xlabel('Real'); ylabel('Imaginary');

As shown in Fig. 4.69, the roots of a Chebyshev filter lie on an ellipse† (see Prob. 4.12-14).

Figure 4.69 Pole-zero plot for an order-8 Chebyshev LPF with fc = 1 kHz and R = 1 dB.

To compute the filter’s magnitude response, the poles are expanded into a polynomial, the dc gain is set based on the even value of N, and CH4MP1 is used.

>> A = poly(C_poles); B = A(end)/sqrt(1+epsilon^2);

omega = linspace(0,2*pi*2000,2001); H_C = CH4MP1(B,A,omega);

>> plot(omega/2/pi,abs(H_C),'k'); axis([0 2000 0 1.1]);

>> xlabel('f [Hz]'); ylabel('|H_C(j2\pi f)|');

^† E. A. Guillemin demonstrates a wonderful relationship between the Chebyshev ellipse and the Butterworth circle in his book Synthesis of Passive Networks (Wiley, New York, 1957).

Figure 4.70 Magnitude responses for an order-8 Chebyshev LPF with fc = 1 kHz and R = 1 dB.

As seen in Fig. 4.70, the magnitude response exhibits correct Chebyshev filter characteristics: passband ripples are equal in height and never exceed R = 1 dB; there are a total of N = 8 maxima and minima in the passband; and the gain rapidly and monotonically decreases after the cutoff frequency of fc = 1 kHz.

For higher-order filters, polynomial rooting may not provide reliable results. Fortunately, Chebyshev roots can also be determined analytically. For

$$\phi_k = \frac{2k+1}{2N}\pi \quad \text{and} \quad \xi = \frac{1}{N}\sinh^{-1}\left(\frac{1}{\epsilon}\right)$$

the Chebyshev poles are

$$p_k = \omega_c \sinh(\xi) \sin(\phi_k) + j\omega_c \cosh(\xi) \cos(\phi_k)$$

Continuing the same example, the poles are recomputed and again plotted. The result is identical to Fig. 4.69.

>> k = [1:N]; xi = 1/N*asinh(1/epsilon); phi = (k*2-1)/(2*N)*pi;
>> C_poles = omega_c*(-sinh(xi)*sin(phi)+1j*cosh(xi)*cos(phi));
>> plot(real(C_poles),imag(C_poles),'kx'); axis equal;
>> axis(omega_c*[-1.1 1.1 -1.1 1.1]);
>> xlabel('Real'); ylabel('Imaginary');

As in the case of high-order Butterworth filters, a cascade of second-order filter sections facilitates practical implementation of Chebyshev filters. Problems 4.12-5 and 4.12-8 use second-order Sallen–Key circuit stages to investigate such implementations.