[3.11 MATLAB: DISCRETE-TIME](#page-9-0) SIGNALS AND SYSTEMS

← Back to LINEAR SYSTEMS AND SIGNALS Overview

3.11 MATLAB: DISCRETE-TIME SIGNALS AND SYSTEMS

MATLAB is naturally and ideally suited to discrete-time signals and systems. Many special functions are available for discrete-time data operations, including the stem, filter, and conv commands. In this section, we investigate and apply these and other commands.

3.11-1 Discrete-Time Functions and Stem Plots

Consider the discrete-time function f[n] = e−n/⁵ cos(πn/5)u[n]. In MATLAB, there are many ways to represent f[n] including M-files or, for particular n, explicit command line evaluation. In this example, however, we use an anonymous function.

f = @(n) exp(-n/5).*cos(pi*n/5).*(n>=0);

A true discrete-time function is undefined (or zero) for noninteger n. Although anonymous function f is intended as a discrete-time function, its present construction does not restrict n to be integer, and it can therefore be misused. For example, MATLAB dutifully returns 0.8606 to f(0.5) when a NaN (not-a-number) or zero is more appropriate. The user is responsible for appropriate function use.

Next, consider plotting the discrete-time function f[n] over (−10 ≤ n ≤ 10). The stem command simplifies this task.

>> n = (-10:10)';
>> stem(n,f(n),'k');
>> xlabel('n'); ylabel('f[n]');

Here, stem operates much like the plot command: dependent variable f(n) is plotted against independent variable n with black lines. The stem command emphasizes the discrete-time nature of the data, as Fig. 3.31 illustrates.

For discrete-time functions, the operations of shifting, inversion, and scaling can have surprising results. Compare f[−2n] with f[−2n + 1]. Contrary to the continuous case, the second is not a shifted version of the first. We can use separate subplots, each over (−10 ≤ n ≤ 10), to help illustrate this fact. Notice that unlike the plot command, the stem command cannot simultaneously plot multiple functions on a single axis; overlapping stem lines would make such plots difficult to read anyway.

>> subplot(2,1,1); stem(n,f(-2*n),'k'); ylabel('f[-2n]');
>> subplot(2,1,2); stem(n,f(-2*n+1),'k'); ylabel('f[-2n+1]'); xlabel('n');

The results are shown in Fig. 3.32. Interestingly, the original function f[n] can be recovered by interleaving samples of f[−2n] and f[−2n+1] and then time-reflecting the result.

Care must always be taken to ensure that MATLAB performs the desired computations. Our anonymous function f is a case in point: although it correctly downsamples, it does not properly upsample (see Prob. 3.11-2). MATLAB does what it is told, but it is not always told how to do everything correctly!

Figure 3.31 f[n] over (−10 ≤ n ≤ 10).

Figure 3.32 f[−2n] and f[−2n+1] over (−10 ≤ n ≤ 10).

3.11-2 System Responses Through Filtering

MATLAB’s filter command provides an efficient way to evaluate the system response of a constant coefficient linear difference equation represented in delay form as

$$\sum_{k=0}^{N} a_k y[n-k] = \sum_{k=0}^{N} b_k x[n-k]$$

\n(3.44)

In the simplest form, filter requires three input arguments: a length-(N + 1) vector of feedforward coefficients [b0,b1,…,bN], a length-(N + 1) vector of feedback coefficients [a0,a1,…,aN], and an input vector.† Since no initial conditions are specified, the output corresponds to the system’s zero-state response.

To serve as an example, consider a system described by y[n]−y[n−1]+y[n−2] = x[n]. When x[n] = δ[n], the zero-state response is equal to the impulse response h[n], which we compute over (0 ≤ n ≤ 30).

b = [1 0 0]; a = [1 -1 1]; >> n = (0:30)’; delta = @(n) 1.0.*(n==0); >> h = filter(b,a,delta(n)); >> clf; stem(n,h,‘k’); axis([-.5 30.5 -1.1 1.1]); >> xlabel(‘n’); ylabel(‘h[n]’);

^† It is important to pay close attention to the inevitable notational differences found throughout engineering documents. In MATLAB help documents, coefficient subscripts begin at 1 rather than 0 to better conform with MATLAB indexing conventions. That is, MATLAB labels a⁰ as a(1), b⁰ as b(1), and so forth.

Figure 3.33 h[n] for y[n] −y[n−1] +y[n−2] = x[n].

Figure 3.34 Resonant zero-state response y[n] for x[n] = cos(2πn/6)u[n].

As shown in Fig. 3.33, h[n] appears to be (N⁰ = 6)-periodic for n ≥ 0. Since periodic signals are not absolutely summable, %^{∞ n}=−∞ |h[n]| is not finite and the system is not BIBO-stable. Furthermore, the sinusoidal input x[n] = cos(2πn/6)u[n], which is (N⁰ = 6)-periodic for n ≥ 0, should generate a resonant zero-state response.

>> x = @(n) cos(2*pi*n/6).*(n>=0);
>> y = filter(b,a,x(n));
>> stem(n,y,'k'); xlabel('n'); ylabel('y[n]');

The response’s linear envelope, shown in Fig. 3.34, confirms a resonant response. The characteristic equation of the system is γ ² − γ + 1, which has roots γ = e±jπ/3. Since the input x[n] = cos(2πn/6)u[n] = (1/2)(ejπn/³ + e−jπn/³)u[n] coincides with the characteristic roots, a resonant response is guaranteed.

By adding initial conditions, the filter command can also compute a system’s zero-input response and total response. Continuing the preceding example, consider finding the zero-input response for y[−1] = 1 and y[−2] = 2 over (0 ≤ n ≤ 30).

>> z_i = filtic(b,a,[1 2]);
>> y_0 = filter(b,a,zeros(size(n)),z_i);
>> stem(n,y_0,'k'); xlabel('n'); ylabel('y_{0} [n]');
>> axis([-0.5 30.5 -2.1 2.1]);

Figure 3.35 Zero-input response y0[n] for y[−1] = 1 and y[−2] = 2.

There are many physical ways to implement a particular equation. MATLAB implements Eq. (3.44) by using the popular direct form II transposed structure.† Consequently, initial conditions must be compatible with this implementation structure. The signal-processing toolbox function filtic converts the traditional y[−1], y[−2], …, y[−N] initial conditions for use with the filter command. An input of zero is created with the zeros command. The dimensions of this zero input are made to match the vector n by using the size command. Finally, _{ } forces subscript text in the graphics window, and ^{ } forces superscript text. The results are shown in Fig. 3.35.

Given y[−1] = 1 and y[−2] = 2 and an input x[n] = cos(2πn/6)u[n], the total response is easy to obtain with the filter command.

$$\Rightarrow y\_total = filter(b,a,x(n),z_i);$$

Summing the zero-state and zero-input response gives the same result. Computing the total absolute error provides a check.

>> sum(abs(y_total-(y + y_0)))
    ans = 1.8430e-014

Within computer round-off, both methods return the same sequence.

3.11-3 A Custom Filter Function

The filtic command is available only if the signal-processing toolbox is installed. To accommodate installations without the signal-processing toolbox and to help develop your MATLAB skills, consider writing a function similar in syntax to filter that directly uses the ICs y[−1], y[−2], …, y[−N]. Normalizing a⁰ = 1 and solving Eq. (3.44) for y[n] yield

$$y[n] = \sum_{k=0}^{N} b_k x[n-k] - \sum_{k=1}^{N} a_k y[n-k]$$

This recursive form provides a good basis for our custom filter function.

^† Implementation structures, such as direct form II transposed, are discussed in Ch. 4.

function [y] = CH3MP1(b,a,x,yi);
% CH3MP1.m : Chapter 3, MATLAB Program 1
% Function M-file filters data x to create y
% INPUTS: b = vector of feedforward coefficients
% a = vector of feedback coefficients
% x = input data vector
% yi = vector of initial conditions [y[-1], y[-2], ...]
% OUTPUTS: y = vector of filtered output data
yi = flipud(yi(:)); % Properly format IC's.
y = [yi;zeros(length(x),1)]; % Preinitialize y, beginning with IC's.
x = [zeros(length(yi),1);x(:)]; % Append x with zeros to match size of y.
b = b/a(1);a = a/a(1); % Normalize coefficients.
for n = length(yi)+1:length(y),
   for nb = 0:length(b)-1,
       y(n) = y(n) + b(nb+1)*x(n-nb); % Feedforward terms.
   end
   for na = 1:length(a)-1,
       y(n) = y(n) - a(na+1)*y(n-na); % Feedback terms.
   end
end
y = y(length(yi)+1:end); % Strip off IC's for final output.

Most instructions in CH3MP1 have been discussed; now we turn to the flipud instruction. The flip up-down command flipud reverses the order of elements in a column vector. Although not used here, the flip left-right command fliplr reverses the order of elements in a row vector. Note that typing help filename displays the first contiguous set of comment lines in an M-file. Thus, it is good programming practice to document M-files, as in CH3MP1, with an initial block of clear comment lines.

As an exercise, the reader should verify that CH3MP1 correctly computes the impulse response h[n], the zero-state response y[n], the zero-input response y0[n], and the total response y[n]+y0[n].

3.11-4 Discrete-Time Convolution

Convolution of two finite-duration discrete-time signals is accomplished by using the conv command. For example, the discrete-time convolution of two length-4 rectangular pulses, g[n] = (u[n]−u[n−4])∗(u[n]−u[n−4]), is a length-(4+4−1=7) triangle. Representing u[n]−u[n−4] by the vector [1, 1, 1, 1], the convolution is computed by

conv([1 1 1 1],[1 1 1 1]) ans = 1 2 3 4 3 2 1

Notice that (u[n+4] −u[n]) ∗ (u[n] −u[n−4]) is also computed by conv([1 1 1 1],[1 1 1 1]) and obviously yields the same result. The difference between these two cases is the regions of support: (0 ≤ n ≤ 6) for the first and (−4 ≤ n ≤ 2) for the second. Although the conv command

312 CHAPTER 3 TIME-DOMAIN ANALYSIS OF DISCRETE-TIME SYSTEMS

does not compute the region of support, it is relatively easy to obtain. If vector w begins at n = nw and vector v begins at n = nv, then conv(w,v) begins at n = nw +nv.

In general, the conv command cannot properly convolve infinite-duration signals. This is not too surprising, since computers themselves cannot store an infinite-duration signal. For special cases, however, conv can correctly compute a portion of such convolution problems. Consider the common case of convolving two causal signals. By passing the first N samples of each, conv returns a length-(2N − 1) sequence. The first N samples of this sequence are valid; the remaining N −1 samples are not.

To illustrate this point, reconsider the zero-state response y[n] over (0 ≤ n ≤ 30) for system y[n]−y[n−1]+y[n−2] = x[n] given input x[n] = cos(2πn/6)u[n]. The results obtained by using a filtering approach are shown in Fig. 3.34.

The response can also be computed using convolution according to y[n] = h[n] ∗ x[n]. The impulse response of this system is†

$$h[n] = \left\{ \cos\left(\frac{\pi n}{3}\right) + \frac{1}{\sqrt{3}} \sin\left(\frac{\pi n}{3}\right) \right\} u[n]$$

Both h[n] and x[n] are causal and have infinite duration, so conv can be used to obtain a portion of the convolution.

>> u = @(n) 1.0.*(n>=0); h = @(n) (cos(pi*n/3)+sin(pi*n/3)/sqrt(3)).*u(n);
>> y = conv(h(n),x(n));
>> stem([0:60],y,'k'); xlabel('n'); ylabel('y[n]');

The conv output is fully displayed in Fig. 3.36. As expected, the results are correct over (0 ≤ n ≤ 30). The remaining values are clearly incorrect; the output envelope should continue to grow, not decay. Normally, these incorrect values are not displayed.

stem(n,y(1:31),‘k’); xlabel(‘n’); ylabel(‘y[n]’);

The resulting plot is identical to Fig. 3.34.

Figure 3.36 y[n] for x[n] = cos(2πn/6)u[n] computed with conv.

^† Techniques to analytically determine ^h[n] are presented in Ch. 5.