|
|
MATLAB Tutorial |
|
4. Discrete-Time System Analysis
Commands covered:
To perform discrete time convolution, x[n]*h[n], define the vectors x and h with elements in the sequences x[n] and h[n]. Then use the command
y = conv(x,h)
This command assumes that the first element in x and the first element in h correspond to n=0, so that the first element in the resulting output vector corresponds to n=0. If this is not the case, then the output vector will be computed correctly, but the index will have to be adjusted. For example,
yields y = [0 1 3 6 6 6 5 3]. If x is indexed as described above, then y[0] = 0, y[1] = 1, .... In general, total up the index of the first element in h and the index of the first element in x, this is the index of the first element in y. For example, if the first element in h corresponds to n = -2 and the first element in x corresponds to n = -3, then the first element in y corresponds to n = -5.
Care must be taken when computing the convolution of infinite duration signals. If the vector x has length q and the vector h has length r, then you must truncate the vector y to have length min(q,r). See the comments in Problem 3.7 of the textbook for additional information.
The command conv can also be used to multiply polynomials: suppose that the coefficients of a(s) are given in the vector a and the coefficients of b(s) are given in the vector b, then the coefficients of the polynomial a(s)b(s) can be found as the elements of the vector defined by ab = conv(a,b).
The command deconv is the inverse procedure to the convolution. In this text, it is used as a means of dividing polynomials. Given a(s) and b(s) with coefficients stored in a and b, then the coefficients of c(s) = b(s)/a(s) are found by using the command c = deconv(b,a).
B. Transfer Function Representation
For a discrete-time transfer function,
the coefficients are stored in descending powers of z or ascending
powers of 
.
For example,
then define the vectors as
Commands Covered:
There are three methods to compute the response of a system described by the following recursive relationship
The first method uses the command
recur and is useful when there are
nonzero initial conditions. This command is available from the
MathWorks ftp site and a shortened version is given in Figure
C.5 of the textbook. The inputs to the function are the coefficients
stored in the
vectors 
and 
,
the initial conditions on x and on y are stored in the vectors
x0 = [x[n0-M], x[n0-M+1],...,x[n0-1]]
and y0 = [y[n0-N], y[n0-N+1],...,y[n0-1]]],
and the time indices for which the solution needs to be calculated
are stored in the vector n
where n0 represents the first element in this vector. To use
recur, type
y = recur(a,b,n,x,x0,y0);
The output is a vector y with elements y[n]; the first element of y corresponds to the time index n0. For example, consider the system described by
y[n] - 0.6y[n-1] + 0.08y[n-2] = x[n-1]
where x[n] = u[n] and with initial conditions y[-1] = 2, y[-2] = 1, and x[-1] = x[-2] = 0. To compute the response y[n] for n = 0, 1,...,10, type
The vector y contains the values of y[n] for n = 0,1,...,10.
The second method to compute the response
uses convolution and is useful when the initial conditions on
y are zero. This method involves first finding the impulse response
of the system, h[n], and then convolving h[n] with x[n] as discussed
in Section 4.A. For example, consider the system described above
with zero initial conditions, that is, y[-1]=y[-2]=0. The impulse
response for this system is 
.
The commands to compute y[n] are
The vector y contains the values of y[n] for n = 0,1,...,10. Note that the vector was truncated to length(n) because both x[n] and h[n] are infinite duration signals. See the comments in Section 4.A regarding the convolution of infinite duration signals.
The third method of solving for the response
requires that the transfer function of the system be known. The
commands dstep and
dimpulse compute the unit step response
and the unit impulse response, respectively while the command
filter computes the response to initial
conditions and to arbitrary inputs. The denominator coefficients
are stored as 
and the numerator coefficients are stored as 
where there are N-M zeros padded on the end
of the coefficients. For example, consider the system given above
with initial conditions y[-1] = y[-2] = 0. To compute the step
response for n=0 to n=10, type the commands
The response can then be plotted using the stem plot. To compute the impulse response, simply replace dstep with dimpulse in the above commands.
To compute the response to an arbitrary input, store the input sequence in the vector x. The command
y = filter(num,den,x);
is used to compute the system
response. If the system has nonzero initial conditions, the initial
conditions can be stored in a vector v0. For a first order system
where N=M=1, define 
.
For a second order system where N=M=2, define 
. To compute the response with nonzero
initial conditions, type
y = filter(num,den,x,zi);
For example, consider the previous system with the initial conditions y[-1] = 2 and y[-2] = 1 and input x[n] = u[n]. Type the following commands to compute y[n].
Commands covered:
freqz
The DTFT of a system can be calculated from the transfer function using freqz. Define the numerator and the denominator of the transfer function in num and den. The command
[H,Omega] = freqz(num,den,n,'whole');
computes the DTFT for n points equally spaced around the unit circle at the frequencies contained in the vector Omega. The magnitude of H is found from abs(H) and the phase of H is found from angle(H). To customize the range for W, define a vector Omega of desired frequencies, for example Omega = -pi:2*pi/300:pi defines a vector of length 301 with values that range from -p to p. To get the DTFT at these frequencies, type
H = freqz(num,den,Omega);
back to list of contents
Commands covered:
The analog prototype method of designing IIR filters can be done by first designing an analog filter with the desired characteristics as shown in Section 3.D, then mapping the filter to the discrete-time domain. Store the numerator and denominator of the analog filter, H(s), in the vectors num and den, and let T be the sampling period. Then the numerator and denominator of the digital filter H(z) is found from the following command
[numd,dend] = bilinear(num,den,1/T)
Alternately, the commands butter and cheby1 automatically design the analog filter and then use the bilinear transformation to map the filter to the discrete-time domain. Lowpass, highpass, bandstop, and bandpass filters can be designed using this method. The digital cutoff frequencies must be specified; these should be normalized by pi. To design a digital lowpass filter based on the analog Butterworth filter, use the commands:
[num,den] = butter(n,Omegac)
where n
is the number of poles and Omegac
is the normalized digital cutoff frequency, 
. To design a highpass filter with cutoff
frequency Omegac, use the
commands
[num,den] = butter(n,Omegac,'high')
To design a bandpass filter with passband from Omega1 to Omega2, define Omega = [Omega1,Omega2] and use the command
[num,den] = butter(n,Omega)
To design a bandstop filter with stopband from Omega1 to Omega2, define Omega = [Omega1, Omega2] and use the command
[num,den] = butter(n,Omega,'stop')
The design for an nth order Type I Chebyshev filter is accomplished using the same methods as for butter except that "butter" is replaced by "cheby1":
If Omega has two elements,
- [num,den] = cheby1(n,Omega); % for a bandpass
- [num,den] = cheby1(n,Omega,'stop'); % for a bandstop
The windows used in FIR filter design are given by
These commands are used to truncate the infinite impulse response of an ideal digital filter with the result being an FIR filter with length N.
The Signal Processing Toolbox also provides commands for computing the FIR filter directly. To obtain an FIR filter with length N and cutoff frequency Omegac (normalized by pi) use the command
h = fir1(N-1,Omegac)
The vector h contains the impulse response of the FIR where h(1) is the value of h[0].
A length N highpass filter with normalized cutoff frequency Omegac is designed by using the command
h = fir1(N-1,Omegac,'high')
A bandpass with passband from Omega1 to Omega2 is obtained by typing
h = fir1(N-1,Omega)
where Omega = [Omega1,Omega2]. A bandstop filter with stopband from Omega1 to Omega2 is obtained by typing
h = fir1(N-1,Omega,'stop')
where Omega = [Omega1,Omega2]. The fir1 command uses the Hamming window by default. Other windows are obtained by adding an option of 'hanning' or 'boxcar' to the arguments; for example,
h = fir1(N-1,Omegac,'high',boxcar(N))
creates a highpass FIR filter with cutoff frequency Omegac using a rectangular window.
Commands covered:
An analog controller G(s) can be mapped to a digital controller D(z) using the bilinear transformation or the step response matching method. Store the numerator and denominator of G(s) in num and den. Then the numerator and denominator of D(z) is found from the bilinear transformation using the commands
[numd,dend] = bilinear(num,den,1/T)
where T is the sampling frequency. To use the step invariant method, use the commands
[numd,dend] = c2dm(num,den,T,'zoh')
To simulate the response of a continuous-time plant G(s) with a digital controller D(z), use the command hybrid. Consider the block diagram in Figure 11.25 of the text. The numerator and denominator coefficients of the plant are stored in numG and denG; the numerator and denominator coefficients of the controller are stored in numD and denD; the reference input signal is stored in r; and the sampling time is stored in T. The increments in the time vector should selected to be the sampling time divided by an integer, for example, t = 0:b:Tend where there is some integer m such that bm=T. The command is used as
[y,ud] = hybrid(numG,denG,numD,denD,T,t,r);
The outputs of the command are the system response, y, and the control signal that is input to the plant, ud. The M-file contains a loop which computes the discrete-time control and then simulates the continuous-time plant for T seconds with the constant control. The process repeats for the next T second interval.
Commands Covered:
Most of the commands for the continuous time state space representation also work for the discrete time state space. For example, ss2tf, tf2ss, and ss2ss for discrete time are used exactly the same way as for the continuous time case discussed in Section 3.F. There is a discrete time version of the command lsim, which is used as follows:
[y,x] = dlsim(A,B,C,D,u,n);
where the output is stored in y, the states are stored in x, the input is stored in u and the time index is stored in n.
|
|
MATLAB Tutorial Table of Contents |
|
This tutorial is available as a supplement to the textbook Fundamentals of Signals and Systems Using Matlab by Edward Kamen and Bonnie Heck, published by Prentice Hall. The tutorial is designed for students using either the professional version of MATLAB (ver. 5.0) with the Control Systems Toolbox (ver. 4.0) and the Signal Processing Toolbox (ver. 3.0), or using the Student Edition of MATLAB (ver. 4.0). For more information on MATLAB, contact The Mathworks, Inc.