ECE-UY 3054: Signals, Systems, and Transforms

Additional Notes for the MATLAB Lab

1
 Notes on Partial Fraction Expansion in Matlab

The Matlab command residue allows one to do partial fraction expansion.

>> help residue

RESIDUE Partial-fraction expansion (residues).

[R,P,K] = RESIDUE(B,A) finds the residues, poles and direct term of
a partial fraction expansion of the ratio of two polynomials B(s)/A(s).
If there are no multiple roots,
B(s) R(1) R(2) R(n)
---- = -------- + -------- + ... + -------- + K(s)
A(s) s - P(1) s - P(2) s - P(n)
Vectors B and A specify the coefficients of the numerator and
denominator polynomials in descending powers of s. The residues
are returned in the column vector R, the pole locations in column
vector P, and the direct terms in row vector K. The number of
poles is n = length(A)-1 = length(R) = length(P). The direct term
coefficient vector is empty if length(B) < length(A), otherwise
length(K) = length(B)-length(A)+1.

If P(j) = ... = P(j+m-1) is a pole of multiplicity m, then the

expansion includes terms of the form
R(j) R(j+1) R(j+m-1)
-------- + ------------ + ... + ------------
s - P(j) (s - P(j))^2 (s - P(j))^m

[B,A] = RESIDUE(R,P,K), with 3 input arguments and 2 output arguments,

converts the partial fraction expansion back to the polynomials with
coefficients in B and A.

See also POLY, ROOTS, DECONV.

2
To perform partial fraction expansion on T (z)

4z + 1
T (z) =
z2 −z−2
we enter the following Matlab commands.

>> num = [4 1];

>> den = [1 -1 -2];
>> [r,p,k] = residue(num,den)

r =
3
1

p =
2
-1

k =
[]

This tells us that,

3 1
T (z) = + .
z−2 z+1

3
The residue command also works when some poles are complex. For example, to find the partial
fraction expansion of X(z)

4z + 3
X(z) =
2 z3 − 3.4 z 2 + 1.98 z − 0.406
we can use the following Matlab commands.

>> num = [4 3];

>> den = [2 -3.4 1.98 -0.406];
>> [r,p,k] = residue(num,den)

r =
36.2500
-18.1250 +13.1250i
-18.1250 -13.1250i

p =
0.7000
0.5000 + 0.2000i
0.5000 - 0.2000i

k =
[]

This tells us that,

36.25 −18.125 + 13.125i −18.125 − 13.125i

X(z) = + + .
z − 0.7 z − (0.5 + 0.2i) z − (0.5 − 0.2i)
Notice that the residues and poles appear in complex-conjugate pairs. In fact, this is always the
case when the coefficients in X(z) are real.

4
Here is an example with a repeated pole. To perform partial fraction expansion on T (z)

2z + 1
T (z) =
z3 + 5z 2 + 8z + 4
we enter the following Matlab commands.

>> num = [2 1];

>> den = [1 5 8 4];
>> [r,p,k] = residue(num,den)

r =
1.0000
3.0000
-1.0000

p =
-2.0000
-2.0000
-1.0000

k =
[]

This tells us that,

1 3 −1
T (z) = + 2
+
z + 2 (z + 2) z+1

5
The residue command also works in the other direction. To write G(z),

5 6 7
G(z) = + −
z − 3 z + 4 z + 1/5
as the ratio of two polynomials we can use the following commands.

>> r = [5 6 -7];
>> p = [3 -4 -1/5];
>> [num,den] = residue(r,p,[])

num =

4.0000 -2.8000 84.4000

den =

1.0000 1.2000 -11.8000 -2.4000

This tells us that

4z 2 − 2.8 z + 84.4
G(z) = .
z 3 + 1.2z 2 − 11.8z − 2.4

The argument k is needed only when the degree of the numerator is greater than or equal to the
degree of the denominator, so it was not needed for these examples.

6
In this example, the degree of the numerator is greater than the degree of the denominator, so k is
required.
3 z 4 − 1.1 z 3 + 0.88 z 2 − 2.396 z + 1.348
X(z) =
z 3 − 0.7 z 2 − 0.14 z + 0.048
The partial fraction expansion of X(z) can be found using the following Matlab commands.

>> num = [3 -1.1 0.88 -2.396 1.348];

>> den = [1 -0.7 -0.14 0.048];
>> [r,p,k] = residue(num,den)

r =
1.0000
4.0000
-3.0000

p =
0.8000
-0.3000
0.2000

k =
3.0000 1.0000

This tells us that

1 4 3
X(z) = 3 z + 1 + + −
z − 0.8 z + 0.3 z − 0.2

7
 Notes on Recursive Digital Filters

A recursive digital filter is an LTI system based on a difference equation of the form:
N
X M
X
y(n) = − a(k) y(n − k) + b(k) x(n − k) (1)
k=1 k=0

A system described by this type of difference equation is called an IIR (Infinite Impulse Response)
filter, a recursive filter, or an autoregressive moving-average (ARMA) filter.
The output y(n) of the filter can be written as
∞
X
y(n) = h(k) x(n − k). (2)
k=0

The transfer function of the system H(z) can be written as

∞
X
H(z) = h(k) z −k . (3)
k=0

The frequency response is given by

∞
X
f
H (ω) = H(e ) = jω
h(n) e−jω . (4)
n=0

The transfer function can also be written in terms of the difference equation coefficients as

b0 + b1 z −1 + b2 z −2 + · · · + bM z −M
H(z) = (5)
1 + a1 z −1 + a2 z −2 + · · · + aN z −N
or
B(z)
H(z) = (6)
A(z)
where
M
X
B(z) = b(n) z −n (7)
n=0
N
X
A(z) = 1 + a(n) z −n . (8)
n=1

Computing h(n). Given the difference equation coefficients a(k) and b(k), the impulse response
h(n) can be obtained by taking the inverse Z-transform of H(z), but it is usually simpler to
calculate h(n) numerically by running the difference equation with the input x(n) = δ(n). The
Matlab function filter implements a recursive difference equation. The following code fragment
computes the first 30 values of h(n) from a(k) and b(k).

8
x = [1, zeros(1,29)];
h = filter(b,a,x);

where the vectors b and a contain the filter parameters

b = [b(0), b(1), . . . , b(M )] (9)

a = [a(0), a(1), . . . , a(N )] (10)

Computing H f (ω) [H(ejω )]. The frequency response can be calculated with the Matlab command
freqz, which uses the FFT (Fast Fourier Transform).
The following code, filt1.m, illustrates how to use these commands to display the frequency
response and impulse response of a digital filter. This code is available of the class web page.

b = [0.0931 0.0960 0.1801 0.1801 0.0960 0.0931];

a = [1.0000 -1.5757 2.2408 -1.5554 0.8123 -0.1837];

% ------ Plot frequency response magnitude ---

subplot(3,1,1)
[H,w]=freqz(b,a,2^10);
plot(w/pi,abs(H))
axis([0 1 0 1.2])
xlabel(’\omega/\pi’)
ylabel(’|H(e^{j \omega})|’)

% ------ Plot frequency response in dB -------

subplot(3,1,2)
plot(w/pi,20*log10(abs(H)))
axis([0 1 -60 10])
xlabel(’\omega/\pi’)
ylabel(’|H(e^{j \omega})| in dB’)

% ------ Plot impulse response h(n) ----------

subplot(3,1,3)
N = 30;
im = [1 zeros(1,N-1)];
h = filter(b,a,im);
stem(0:N-1,h,’.’)
xlabel(’n’)
ylabel(’h(n)’)
axis([-1 30 -0.2 0.6])

This code produces the following figure.

9
 1

|H(ej ω)|
0.5

0
0 0.2 0.4 0.6 0.8 1
ω/π
0

|H(ej ω)| in dB −20

−40

−60
0 0.2 0.4 0.6 0.8 1
ω/π

0.4
h(n)

0.2

−0.2
0 10 20 30
n

Because the ideal low-pass filter can not be implemented, it is necessary to settle for an approxima-
tion to it. The problem of digital filter design is choose coefficients a(k), b(k), so that the frequency
response of the filter H f (ω) matches a desired target function.
A good approximation to the ideal low-pass filter is the elliptic filter, also called the Cauer filter.
The coefficients a(k), b(k) of elliptic filters can be obtained with the Matlab command ellip. (It
is part of the signal processing toolbox.) Here is the beginning of the help file for this command.

ELLIP Elliptic or Cauer digital and analog filter design.

[B,A] = ELLIP(N,Rp,Rs,Wn) designs an Nth order lowpass digital
elliptic filter with Rp decibels of ripple in the passband and
a stopband Rs decibels down. ELLIP returns the filter
coefficients in length N+1 vectors B (numerator) and A (denominator).
The cut-off frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding
to half the sample rate. Use Rp = 0.5 and Rs = 20 as starting
points, if you are unsure about choosing them.

The meaning of Rp and Rs is given by

δp = 1 − 10−Rp /20 (11)

δs = 10−Rs /20 (12)

10
and

Rp = −20 log10 (1 − δp ) (13)

Rs = −20 log10 (δs ) (14)

where δp represents the size of the ripple in the pass-band, and δs represents the size of the ripple
in the stop-band. Rp , and Rs are just δp and δs in decibels. The pass-band ripple size δp is the
maximum deviation of the actual frequency response from 1 in the pass-band. The stop-band ripple
size δs is the maximum deviation of the actual frequency response from 0 in the stop-band.
The Matlab filter design commands use a somewhat non-standard convention for specifying the cut-
off frequency ωc . For example, if you want to design a discrete-time filter with cut-off frequency at
ωc = 0.4 π, then you should set Wn = 0.4 when using the ellip command. The same is true for
the other Matlab digital filter design commands.
By the way, this assumes that the sampling period is 1 second. That is, the values x(n) are
obtained from sampling a continuous-time signal once each second. If the sampling period (call it
Ts ) is different from 1 second, and you want to design a digital filter with cut-off frequency at Fo
Hz. Then you should take Wn = 2*Fo/Fs. Here Fs is the sampling frequency in Herz, Fs = 1/Ts .
For example, suppose we sample an analog signal with a sampling rate of 100 Hz and we want to
design a digital filter with a cut-off frequency of 30 Hz. Here are the Matlab commands to set up
the design parameters, call the function, and plot the frequency response versus frequency in Herz.

Ts = 0.01; % Sampling Period Ts = 0.01 seconds

Fs = 1/Ts % Sampling frequency Fs = 100 Hz
Fo = 30; % Cut-Off frequency Fo = 30 Hz

Wn = 2*Fo/Fs; % Normalize the cut-off frequency

dp = 0.01;
ds = 0.01;
N = 5;
Rp = -20*log10(1-dp);
Rs = -20*log10(ds);

[b,a] = ellip(N,Rp,Rs,Wn);
figure(1), clf

% Plot frequency response magnitude

[H,w] = freqz(b,a);
plot(w/pi*Fs/2,abs(H)) % rescale frequency axis appropriately.
axis([0 Fs/2 0 1.2])
xlabel(’FREQUENCY (Hz)’)
ylabel(’FREQUENCY RESPONSE’)
grid

11
orient portrait
print -deps filt2

1.2

0.8
FREQUENCY RESPONSE

0.6

0.4

0.2

0
0 5 10 15 20 25 30 35 40 45 50
FREQUENCY (Hz)

12
 Notes on Stable Inverses of Discrete-Time LTI Systems

Suppose the signal x(n) represents binary data that is to be transmitted.

1
x(n)

−1
0 20 40 60 80 100 120 140 160 180 200

However, the channel over which x(n) is transmitted introduces some distortion. The signal r(n)
that is received is changed.

1
r(n)

−1
0 20 40 60 80 100 120 140 160 180 200

In this example, the channel can be modeled as an LTI system described by the difference equation

r(n) − 11 r(n − 1)/20 − r(n − 2)/5 = x(n)/10 + 7 x(n − 1)/30 + 2 x(n − 2)/30

where r(n) represents the received signal. Equivalently, the channel can be modeled by the transfer
function
1 1 + 73 z −1 + 23 z −2
H(z) = · , (15)
10 1 − 11
20 z
−1 − 1 z −2
5

1 z 2 + 37 z + 23
H(z) = · 1, (16)
10 z 2 − 11
20 z − 5
or equivalently by its impulse response:

0.4

0.2
h(n)

−0.2
0 2 4 6 8 10 12 14 16 18 20

Define
1
G(z) := . (17)
H(z)

13
Then
z 2 − 11
20 z − 5
1
G(z) = 10 · . (18)
z 2 + 73 z + 23
A causal system with the transfer function G(z) can be implemented by the difference equation

y(n)/10 + 7 y(n − 1)/30 + 2y(n − 2)/30 = r(n) − 11 r(n − 1)/20 − r(n − 2)/5

where r(n) represents the input to the system and y(n) represents the output of the system. This
system can be implemented in a causal fashion by writing the difference equation as

y(n) = 10 r(n) − 11 r(n − 1)/2 − 2 r(n − 2) − 7 y(n − 1)/3 − 2 y(n − 2)/3.

However, when this is done, the signal y(n) that you obtain is:

1
y(n)

−1
0 20 40 60 80 100 120 140 160 180 200

Because the causal implementation of this difference is unstable, the signal blows up after a short
while. In order to obtain a stable inverse for the system H(z), we need to be more careful. Let us
find a partial fraction expansion of G(z)/z

G(z) 10z 2 − 11 z − 2
= 3 7 22 2 (19)
z z + 3z + 3z

We can use Matlab to help with the partial fraction expansion.

>> [r,p,k] = residue([10 -11/2 -2],[1 7/3 2/3 0])

r =
14.7000
-1.7000
-3.0000
p =
-2.0000
-0.3333
0
k =
[]

14
Therefore, we have
G(z) 14.7 −1.7 −3
= + + (20)
z z + 2 z + 1/3 z
or
14.7 z −1.7 z
G(z) = + −3 (21)
z+2 z + 1/3
and    


 14.7 (−2)n u(n) 




 −1.7 (−1/3)n u(n) 


g(n) = or + or − 3 δ(n).
   
−14.7 (−2)n u(−n − 1) 1.7 (−1/3)n u(−n − 1)

 
 
 


Consider the first component — the causal signal 14.7 (−2)n u(n) grows without bound:

14.7 (−2)n u(n)

100

50

−50

−100
−5 0 5

On the other hand, the noncausal signal −14.7 (−2)n u(−n − 1) decays:

−14.7 (−2)n u(−n−1)

10

−5

−10
−5 0 5

To obtain a stable system, we must choose the noncausal form for this component.
Consider the second component — the causal signal −1.7 (−1/3)n u(n) decays:

−1.7 (−1/3)n u(n)

2

−1

−2
−5 0 5

On the other hand, the noncausal signal 1.7 (−1/3)n u(−n − 1) grows without bound:

15
 1.7 (−1/3)n u(−n−1)
100

50

−50

−100
−5 0 5

To obtain a stable system, we must choose the causal form for this component.
Therefore, to obtain a stable system, the impulse response g(n) must be given by

g(n) = −14.7 (−2)n u(−n − 1) − 1.7 (−1/3)n u(n) − 3 δ(n). (22)

This impulse response g(n) is shown in the figure.

10
g(n) (STABLE)

−5

−10
−10 −8 −6 −4 −2 0 2 4 6 8 10

For an inverse system, it is required that h(n) ∗ g(n) = δ(n). We can verify this by numerically
convolving the two functions.

1
h(n) conv g(n)

0.5

−0.5
−10 −8 −6 −4 −2 0 2 4 6 8 10

We can now recover x(n) from r(n) by using the LTI system with impulse response g(n). Note
that even though g(n) is noncausal, for n < −10, the values of g(n) are very small and for practical
purposes can be set to zero. If they are set to zero, then the shifted signal g(n − 10) will be causal.
So you can define a causal system with impulse response gc (n):

 g(n − 10) n ≥ 0
gc (n) := (23)
 0 n<0

or more compactly as
gc (n) := g(n − 10) u(n). (24)

16
This is shown in the figure:

gc(n) (STABLE and CAUSAL)

10

−5

−10
0 2 4 6 8 10 12 14 16 18 20

The signal gc (n) is a good approximation to g(n) because the values of g(n) that are lost due to
truncation are close to zero.
Using the causal stable system gc (n) to filter the received signal r(n) gives:

1
x(n) and y2(n)

−1

−2

−3
0 20 40 60 80 100 120 140 160 180 200

The top signal is x(n). The lower signal is y2 (n), obtained by convolving r(n) with gc (n). Notice
that y2 (n) is almost the same as x(n) except for a delay of 10 samples.

Therefore, we can get a stable near-inverse of the channel if we allow for some delay.

The following Matlab code, used for this problem, can be found on the course webpage.

% generate a data signal (don’t worry about understanding this command)

x = kron(round(rand(1,40)),ones(1,5));

% ------ Plot data signal ------

figure(1)
subplot(4,1,1), plot(x), ylabel(’x(n)’)
axis([0 200 -1 2])

% ------ Define System I ------

b = [1 +7/3 2/3]*(1/10);
a = [1 -11/20 -4/20];

r = filter(b,a,x);

17
%% r(n) is the distorted version of x(n)

subplot(4,1,2), plot(r), ylabel(’r(n)’)

y = filter(a,b,r);

%%

% y(n) is what you get when you try the "simple" thing.
% The "simple" thing is to reverse the difference equation.

subplot(4,1,3), plot(y), ylabel(’y(n)’)

% --> the "simple" thing is unstable!

im = [1 zeros(1,100)]; % define impulse signal

h = filter(b,a,im); % impulse response of System I

% ------ Let G(z) = 1/H(z) ------

% ------ Compute partial fraction expansion of G(z)/z ------

[R,P,K] = residue([10 -11/2 -2],[1 7/3 2/3 0]);

%
% R =
% 14.7000
% -1.7000
% -3.0000
%
% P =
% -2.0000
% -0.3333
% 0
%
% K =
% []

% so g(n) = -14.7 (-2)^n u(-n-1) - 1.7 (-1/3)^n u(n) - 3 delta(n)

% ------ create stable signal g(n) corresponding to G(z) ------

n = -10:10;
g1 = -14.7 .* (-2).^n .* (-n-1>=0);
g2 = -1.7 .* (-1/3).^n .* (n>=0);
g3 = -3*(n==0);
g = g1 + g2 + g3;

% ------ plot h(n) and stable g(n) ------

% ------ also plot h(n) convolve g(n) ----

figure(2)
subplot(4,1,1), stem(h,’.’), ylabel(’h(n)’)
subplot(4,1,2), stem(g,’.’), ylabel(’stable g(n)’)
chck = conv(h,g);
subplot(4,1,3), stem(chck,’.’), ylabel(’h(n) conv g(n)’)

18
% ------ send r(n) through the system G(z) ------

y2 = conv(r,g);
figure(1)
subplot(4,1,4), plot(y2); ylabel(’y_2(n)’)
axis([0 200 -1 2])

% ------ plot x(n) and recovered signal y2(n) ----

% ------ on same graph ------

figure(3)
k = 1:200;
plot(k,x(k),k,y2(k)-2)
axis([0 200 -3 2])
ylabel(’x(n) and y_2(n)’)

% notice that y2(n) is a delayed version of x(n)

% ---> with a delay we can get x(n) back from y(n) using
% a STABLE system.

19