Solution Manual For Additional Problems For Signals and Systems Using Matlab Luis F. Chaparro and Aydin Akan
Solution Manual For Additional Problems For Signals and Systems Using Matlab Luis F. Chaparro and Aydin Akan
Solution Manual For Additional Problems For Signals and Systems Using Matlab Luis F. Chaparro and Aydin Akan
Problems for
SIGNALS AND SYSTEMS
USING MATLAB
0.1 Consider the following problems about trigonometric and polar forms.
(a) Let z = 6ejπ/4 find (i) Re(z), (ii) Im(z)
(b) If z = 8 + j3 and v = 9 − j2, is it true that
1
Chaparro-Akan — Signals and Systems using MATLAB 0.2
0.2 Using the vectorial representation of complex numbers it is possible to get some interesting inequalities.
(a) Is it true that for a complex number z = x + jy we have that |x| ≤ |z|? Show it geometrically by
representing z as a vector.
(b) The so called triangle inequality says that for any complex (or real) numbers z and v we have that
|z + v| ≤ |z| + |v|. Show this geometrically.
(c) If z = 1 + j and v = 2 + j is it true that
Answer: (a) |x| = |z|| cos(θ)| and since | cos(θ)| ≤ 1 then |x| ≤ |z|; (c) yes to both.
Solution
(a) Representing the complex number z = x + jy = |z|ejθ then |x| = |z|| cos(θ)| and since | cos(θ)| ≤ 1
then |x| ≤ |z|, the equality holds when θ = 0 or when z = x, i.e., it is real.
(b) Adding two complex numbers is equivalent to adding two vectors to create a triangle with two sides
the two vectors being added and the other side the vector resulting from the addition. Unless the two
vector being added have the same angle, in which case |z| + |v| = |z + v|, it holds that |z| + |v| > |z + v|.
!z + !v !v
!z
Solution
(a) i. 2 cos(α + β) = ej(α+β) + e−j(α+β) = (ejα ejβ ) + (ejα ejβ )∗ = 2Re(ejα ejβ ) and
Re[ejα ejβ ] = Re[(cos(α) + j sin(α))(cos(β) + j sin(β))] = cos(α) cos(β) − sin(α) sin(β) so that
ii. 2j sin(α + β) = ejα ejβ − (ejα ejβ )∗ = 2jIm[ejα ejβ ], and the imaginary is
sin(α) cos(β) + cos(α) sin(β) = sin(α + β)
(b)
1
ej2πt 1 ej2π − 1
Z
ej2πt dt = |0 = =0
0 j2π j2π
also Z 1 Z 1 Z 1
ej2πt dt = cos(2πt)dt + j sin(2πt)dt = 0 + j0
0 0 0
since the integrals of the sinusoids are over a period.
(c) i. Yes, (−1)n = (ejπ )n = ejnπ = cos(nπ) + j sin(nπ) = cos(nπ) since sin(nπ) = 0 for any integer
n.
ii. Yes, ej0 = −ejπ and ejπ/2 = −ej3π/2 so they add to zero.
0.4 Consider the following problems related to computation with complex numbers.
(a) Find and plot all roots of
(b) Suppose you want to find the natural log of a complex number z = |z|ejθ , which is
If z is negative it can be written as z = |z|ejπ and we can find log(z) by using the above derivation.
The natural log of any complex number can be obtained this way also. Justify each one of the steps
in the above equation and find
(i) find z 2 /4, (ii) is it true that (cos(π/4) + j sin(π/4))2 = cos(π/2) + j sin(π/2) = j?
Answers: (a) (i) zk = ejπ(2k+1)/3 ; (b) log(2ejπ/4 ) = log(2) + jπ/4; (c) elog(z) = −1; (d) (i) z 2 /4 = j.
Solution
(a) We have
i. z 3 = −1 = ejπ(2k+1) for k = 0, 1, 2, so roots are zk = ejπ(2k+1)/3 , k = 0, 1, 2, i.e., on a circle of
unit radius and separated π/3 with one at −1.
ii. z 2 = 1 = ej2πk for k = 0, ±1, ±2, · · · , so roots are zk = ejπk =, i.e., on a circle of unit radius
and separated π , one at zero and the other at −1
p
iii. z 2 + 3z + 1 = (z + 1.5)2 + (1 − 1.52 ) = 0 so the roots are z1,2 = ± (1.52 − 1) − 1.5
(b) The log (using the Naperian base) of a product is the sum of the logs of the terms in the product, and
the log and the exponential are the opposite of each other so the last term in the equation.
Using the given expression for the log of a complex number (log of real numbers is a special case):
Determine for which of the given values of Ts the discrete-time signal has lost the information in the
continuous-time signal. Use MATLAB to plot x(t) (use Ts = 10−4 and the plot function) and the resulting
discrete-time signals (use the stem function). Superimpose the analog and the discrete-time signals for
0 ≤ t ≤ 3; use subplot to plot the four figures as one figure. You also need to figure out how to label the
different axes and how to have the same scales and units.
Answers: For Ts = 1 the discrete signal has lost information.
Solution
As we will see later, the sampling period of x(t) with a frequency of Ωmax = 2πfmax = 2π should satisfy
the Nyquist sampling condition
1
fs = ≥ 2fmax = 2 samples/sec
Ts
so Ts ≤ 1/2 (sec/sample). Thus when Ts = 0.1 the continuous-time and the discrete-time signals look
very much like each other, indicating the signals have the same information — such a statement will be
justified in the chapter on sampling where we will show that the continuous-time signal can be recovered
from the sampled signal. It is clear that when Ts = 1 the information is lost. Although it is not clear
from the figure that when we let Ts = 0.5 the discrete-time signal keeps the information, this sampling
period satisfies the Nyquist sampling condition and as such the original signal can be recovered from the
sampled signal. The following MATLAB script is used.
% Pr. 0.5
clear all; clf
T=3; Tss= 0.0001; t=[0:Tss:T];
xa=4*cos(2*pi*t); % continuous-time signal
xamin=min(xa);xamax=max(xa);
figure(1)
subplot(221)
plot(t,xa); grid
title(’Continuous-time Signal’); ylabel(’x(t)’); xlabel(’t sec’)
axis([0 T 1.5*xamin 1.5*xamax])
N=length(t);
for k=1:3,
if k==1,Ts= 0.1; subplot(222)
t1=[0:Ts:T]; n=1:Ts/Tss: N; xd=zeros(1,N); xd(n)=4*cos(2*pi*t1);
plot(t,xa); hold on; stem(t,xd);grid;hold off
axis([0 T 1.5*xamin 1.5*xamax]); ylabel(’x(0.1 n)’); xlabel(’t’)
elseif k==2, Ts=0.5; subplot(223)
t2=[0:Ts:T]; n=1:Ts/Tss: N; xd=zeros(1,N); xd(n)=4*cos(2*pi*t2);
plot(t,xa); hold on; stem(t,xd); grid; hold off
axis([0 T 1.5*xamin 1.5*xamax]); ylabel(’x(0.5 n)’); xlabel(’t’)
else,Ts=1; subplot(224)
t3=[0:Ts:T]; n=1:Ts/Tss: N; xd=zeros(1,N); xd(n)=4*cos(2*pi*t3);
plot(t,xa); hold on; stem(t,xd); grid; hold off
axis([0 T 1.5*xamin 1.5*xamax]); ylabel(’x(n)’); xlabel(’t’)
end
end
Analog Signal
6 6
4 4
2 2
x(0.1 n)
x(t)
0 0
−2 −2
−4 −4
−6 −6
0 1 2 3 0 1 2 3
t sec t
6 6
4 4
2 2
x(0.5 n)
x(n)
0 0
−2 −2
−4 −4
−6 −6
0 1 2 3 0 1 2 3
t t
Figure 2: Problem 5: Analog continuous-time signal (top left); continuous-time and discrete-time
signals superposed for Ts = 0.1 sec (top right) and Ts = 0.5 sec and Ts = 1 sec (bottom left to right).
0.6 Differential and difference equations — Find the ordinary differential equation relating a
current source is (t) = cos(Ω0 t) with the current iL (t) in an inductor, with inductance L = 1
Henry, connected in parallel with a resistor of R = 1 Ω (see Fig. 3). Assume a zero initial
current in the inductor.
iL (t)
is (t) 1Ω 1H
(a) Obtain a discrete equation from the ordinary differential equation using the trapezoidal
approximation of an integral.
(b) Create a MATLAB script to solve the difference equation for Ts = 0.01 and for three fre-
quencies for is (t), Ω0 = 0.005π, 0.05π and 0.5π. Plot the input current source is (t) and the
approximate solution iL (nTs ) in the same figure. Use the MATLAB function plot. Use the
MATLAB function filter to solve the difference equation (use help to learn about filter).
(c) Solve the ordinary differential equation using symbolic MATLAB when the input fre-
quency is Ω0 = 0.5π. Compare is (t) and iL (t).
Answer: diL (t)/dt + iL (t) = is (t)
vL (t)
is (t) = iR (t) + iL (t) = + iL (t)
R
but vL (t) = LdiL (t)/dt so that the ordinary differential equation relating the input is (t) to the
output current in the inductor iL (t) is
diL (t)
+ iL (t) = is (t)
dt
after replacing L = 1 and R = 1. Notice that this d.e. is the dual of the one given in the Chapter,
so that the difference equation is
Ts 2 − Ts
iL (nTs ) = [is (nTs ) + is ((n − 1)Ts )] + iL ((n − 1)Ts ) n≥1
2 + Ts 2 + Ts
iL (0) = 0
(b)(c) The scripts to solve the difference and ordinary differential equations are the following.
% Pr. 0.6
clear all
cos(1/2 π t)
1
input current 1
is(t),iL(t) output current
0.5 filter gain 0.5
0
0
0 10 20 30 40 50 60 70 80 90 100
t −0.5
1 −1
input current
is(t),iL(t)
output current
0 filter gain
−0.5
−1
0 10 20 30 40 50 60 70 80 90 100 −1
t 0 10 20 30 40 50 60 70 80 90 100
t
Figure 4: Problem 6: Left (top to bottom): solution of difference equation for Ω0 = 0.005, 0.05, 0.5
(rad/sec). Right: input (top), solution of ordinary differential equation (bottom).
0.7 More integrals and sums — Although sums behave like integrals, because of the discrete
nature of sums one needs to be careful with the upper and lower limits more than in the integral
case. To illustrate this consider the separation of an integral into two integrals and compare it
with the separation of a sum into two sums. For the integral we have that
Z 1 Z 0.5 Z 1
tdt = tdt + tdt
0 0 0.5
Show that this this true by computing the three integrals. Then consider the sum
100
X
S= n
n=0
find this sum and determine which of the following is equal to this sum
50
X 100
X 50
X 100
X
(i) S1 = n+ n, (ii) S2 = n+ n
n=0 n=50 n=0 n=51
Solution
The indefinite integral equals 0.5t2 . Computing it in [0, 1] gives the same value as the sum of
the integrals computed between [0, 0.5] and [0.5, 1].
As seen before, the sum
100
X 100(101)
S= n= = 5050
n=0
2
while
S1 = S + 50 = 5100
S2 = S
the first sum has an extra term when n = 50 while the other does not. To verify this use the
following script:
% Pr. 0.7
clear all
N=100;
syms n,N
S=symsum(n,0,N)
S1=symsum(n,0,N/2)+symsum(n,N/2,N)
S2=symsum(n,0,N/2)+symsum(n,N/2+1,N)
giving
S = 5050
S1 = 5100
S2 = 5050
0.8 Complex functions of time — Consider the complex function x(t) = (1 + jt)2 for −∞ < t < ∞.
(a) Find the real and the imaginary parts of x(t) and carefully plot them with MATLAB. Try
to make MATLAB plot x(t) directly, what do you get? Does MATLAB warn you? Does it
make sense?.
(b) Compute the derivative y(t) = dx(t)/dt and plot its real and imaginary parts, how do
these relate to the real and the imaginary parts of x(t)?
(c) Compute the integral
Z 1
x(t)dt
0
(d) Would the following statement be true? (remember ∗ indicates complex conjugate)
Z 1 ∗ Z 1
x(t)dt = x∗ (t)dt
0 0
R1
Answers: 0
x(t)dt = 2/3 + j1; statement is true.
Solution
(a)(b) Since
x(t) = (1 + jt)2 = 1 + j2t + j 2 t2 = 1 − t2 +j |{z}
2t
| {z }
real imag.
its derivative with respect to t is
When plotting the complex function x(t) as function of t, MATLAB ignores the imaginary part.
One should not plot complex functions as functions of time as the results are not clear when
using MATLAB. See Fig. 5 for plots.
(c) Using the rectangular expression of x(t) we have
Z 1 Z 1 Z 1 Z 1
t3 2t2 1 2
x(t)dt = (1 − t2 + 2jt)dt = (1 − t2 )dt + 2j t dt = t − +j 0 = + j1
0 0 0 0 3 2 3
which is the complex conjugate of the integral calculated in (c). So yes, the expression is true.
−5
0
−10
−15
−20 −5
−25
−5 −4 −3 −2 −1 0 1 2 3 4 5
−10
Imaginary part of x(t)
10
−15
5
0
−20
−5
−10 −25
−5 −4 −3 −2 −1 0 1 2 3 4 5 −5 −4 −3 −2 −1 0 1 2 3 4 5
Time Time
Figure 5: Problem 8: Real and imaginary parts of x(t) (left); complex signal x(t) ignoring imaginary
part.
0.9 Hyperbolic sinusoids — In filter design you will be asked to use hyperbolic functions. In this
problem we relate these functions to sinusoids and obtain a definition of these functions so that
we can actually plot them.
(a) Consider computing the cosine of an imaginary number, i.e., use
ejx + e−jx
cos(x) =
2
let x = jθ and find cos(x). The resulting function is called the hyperbolic cosine or
cos(jθ) = cosh(θ)
(b) Consider then the computation of the hyperbolic sine sinh(θ), how would you do it? Care-
fully plot it as a function of θ.
(c) Show that the hyperbolic cosine is always positive and bigger than 1 for all values of θ
(d) Show that sinh(θ) = − sinh(−θ)
(e) Write a MATLAB script to compute and plot these functions between −10 and 10.
Answers: sinh(θ) = −j sin(jθ) and odd.
Solution
(a) If x = jθ
1 −θ
cos(jθ) = (e + eθ ) = cosh(θ)
2
(b) The hyperbolic sine is defined as
1 θ
sinh(θ) = (e − e−θ )
2
which is connected with the circular sine as follows
1 −θ
sin(jθ) = (e − eθ ) = j sinh(θ) ⇒ sinh(θ) = −j sin(jθ)
2j
(c) Since e±θ > 0 then cosh(θ) = cosh(−θ) > 0, the smallest value is for θ = 0 which gives
cosh(0) = 1
(d) Indeed,
1
sinh(−θ) = (e−θ − eθ ) = − sinh(θ)
2
% Pr. 0.9
clear all
theta=sym(’theta’);
x= 0.5*(exp(-theta)+exp(theta));
y= 0.5*(exp(theta)-exp(-theta));
figure(9)
subplot(211)
ezplot(x,[-10,10])
grid
subplot(212)
ezplot(y,[-10,10])
grid
2000
1000
−1000
−2000
−10 −8 −6 −4 −2 0 2 4 6 8 10
θ
Continuous–time Signals
Solution
(a) If x(t) = t for 0 ≤ t ≤ 1, then x(t + 1) is x(t) advanced by 1, i.e., shifted to the left by 1 so that x(0) = 0
occurs at t = −1 and x(1) = 1 occurs at t = 0.
x(t) x(t + 1)
1 1
t t
0 1 −1 0
x(−t) x(−t + 1)
1 1
t t
−1 0 0 1 2
Figure 1.1: Problem 1: Original signal x(t), shifted versions x(t + 1), x(−t) and x(−t + 1).
1
Chaparro-Akan — Signals and Systems using MATLAB 1.2
The signal x(−t) is the reversal of x(t) and x(−t + 1) would be x(−t) advanced to the right by 1. Indeed,
t x(−t + 1)
1 x(0)
0 x(1)
−1 x(2)
The sum y(t) = x(t + 1) + x(−t + 1) is such that at t = 0 it is y(0) = 2; y(t) = x(t + 1) for t < 0; and
y(t) = x(−t + 1) for t > 0. Thus,
or
t+1 −1 ≤ t < 0
y(t) = 2 t=0
−t + 1 0<t≤1
(b) Except for the discontinuity at t = 0, y(t) looks like the even triangle signal Λ(t), their integrals are
y(t)
−1 1 t
Figure 1.2: Problem 1: Triangular signal y(t) with discontinuity at the origin.
Solution
(a) x(t) is called causal because it is zero for t < 0, it repeats every 0.5 sec. for t ≥ 0.
x(t)
1
···
0.5 1
xe (t)
0.5
··· ···
−0.5 0.5 1
xo (t)
0.5
−0.5 ···
··· 0.5 1
−0.5
(b) Even component xe (t) = 0.5[x(t) + x(−t)] is periodic of fundamental period Te = 0.5. The odd
component is xo (t) = 0.5[x(t) − x(−t)] is not periodic.
(c) xe (t) and xo (t) are non-causal signals as they are different from zero for negative times.
Answers: (a) and (b) are true; cos(2πt)δ(t − 1) = δ(t − 1); (g) yes, Pxe = 0.5Px .
Solution
(a) Yes, expressing ej2πt = cos(2πt) + j sin(2πt), periodic of fundamental period T0 = 1, then the
integral is the area under the cosine and sine in one or more periods (which is zero) when k 6= 0 and
integer. If k = 0, the integral is also zero.
(b) Yes, whether t0 = 0 (first equation) or a value different from zero, the two integrals are equal as the
area under a period is the same. In the case x(t) = cos(2πt), both integrals are zero.
(c) It is not true, cos(2πt)δ(t − 1) = cos(2π)δ(t − 1) = δ(t − 1).
(d) It is true, considering x(t) the product of cos(t) and u(t) its derivative is
(e) Yes,
Z ∞ Z ∞
e−t u(t) δ(t − 2)dτ
−2
= e δ(t − 2)dτ
−∞ 0
= e−2
(f) Yes,
dx(t)
= 0.5[et u(t) + et δ(t)] + 0.5[−e−t u(t) + e−t δ(t)]
dt
= 0.5[et − e−t ]u(t) + δ(t) = sinh(t)u(t) + δ(t)
(g) The even component xe (t) is a periodic full-wave rectified signal of amplitude 1/2 and fundamental
period T1 = π.
Power of x(t)
Z π
1
Px = 0.5 x2 (t)dt
π 0
Power of xe (t)
Z π
1
Pxe = (0.5x(t))2 dt = 0.5Px
π 0
Answers: y(t) = 2δ(t + 2) − u(t + 2) + 2u(t) − u(t − 2) − 2δ(t − 2); yes, it is true.
Solution
(a) See Fig. 4a
x(t) = |t| [u(t + 2) − u(t − 2)] Derivative
| {z }
p(t)
x(t)
t
−2 2
y(t)
1
(2)
t
−2 2
−1 (−2)
dx(t)
y(t) = = 2δ(t + 2) − u(t + 2) + 2u(t) − u(t − 2) − 2δ(t − 2)
dt
(b) Integral
0 t < −2
t
Z
−t −2 ≤ t < 0
0 0
y(t )dt =
−∞
t 0≤t<2
0 t≥2
and use MATLAB to plot it for 0 ≤ t ≤ 10. In a simple way this problem illustrates the operation of a
discrete to analog converter which converts a discrete-time into a continuous-time signal (its cousin
is the digital to analog converter or DAC).
Rt
Answers: (a) ss(t) = ∞
P P
k=0 u(t − kT ); (c) −∞ xs (t)dt = k x(kTs )u(t − kTs ).
Solution
(a) The sampling signal is a sequence of unit impulses at uniform times kT , for t > 0. The integral
Z t X ∞ ∞ Z t
X ∞
X
ssT (t) = δ(t − kT )dt = δ(t − kT )dt = u(t − kT )
−∞ k=0 k=0 −∞ k=0
This signal is called the “stairway to the stars” (ssT (t)) for obvious reasons.
(b)(c) The following script will display the signal ss(t) and the conversion to analog (just uncomment the
desired signal and comment the not desired signal).
% Pr. 1_5 parts (b) and (c)
clear all; clf
T=0.1; t=0:T:10;
%f=t; % part (b)
f=cos(2*pi*t); % part (c)
figure(6)
stairs(t,f);grid
%axis([0 10 0 10]) % part (b)
axis([0 10 -1.1 1.1]) % part (c)
hold on
plot(t,f,’r’)
hold off
When f (t) = t + 1 the output looks like a digitized line with unit slope and cut at 1 (see figure
below), similarly when f (t) = cos(2πt) the output looks like a digitized sinusoid.
10
1
9
0.8
8
0.6
7
0.4
6 0.2
5 0
4 −0.2
−0.4
3
−0.6
2
−0.8
1
−1
0
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
Figure 1.5: Problem 5: ‘Digitized’ ramp and cosine signals using ’stairway to stars’
Continuous–time Systems
(a) Let the input be x(t) = sin(2πt)u(t), plot the corresponding output y(t). What is the
output of the system if we double the input? That is, what is the output y1 (t) if the input
x1 (t) = 2x(t) = 2 sin(2πt)u(t)? Is the system linear or not?
(b) Is the given system time-invariant?
P∞
Answers: output y(t) = k=0 p(t − k) with p(t) = u(t) − 2u(t − 0.5) + u(t − 1).
Solution
(a) x(t) = sin(2πt)u(t) is zero for t < 0 and repeats periodically every T0 = 1. Thus,
∞
X
y(t) = p(t − k), p(t) = u(t) − 2u(t − 0.5) + u(t − 1)
k=0
If the input is 2x(t) the output is the same as before, so the system is non-linear.
(b) The system is time-invariant: if input x(t − τ ) the output is y(t − τ ) for any value of τ .
1
Chaparro-Akan — Signals and Systems using MATLAB 2.2
(a) A system is represented by the ordinary differential equation dz(t)/dt = w(t) − w(t − 1)
where w(t) is the input and z(t) the output.
i. How is this system related to an averager having an input/output equation
Z t
z(t) = w(τ )dτ + 2 ?
t−1
ii. Is the system represented by the given ordinary differential equation LTI?
(b) If we consider the current i(t) the input of a capacitor, with a zero initial voltage across it,
the voltage across the capacitor –or the output— is given by
Z t
vC (t) = i(τ )dτ
0
i. Is this capacitor time-invariant? If not, what conditions would you impose to make it
time-invariant?
ii. If we let i(t) = u(t) what is the corresponding voltage vC (t)? If we delay the current
source so that the input is i(t − 1), find the corresponding voltage across the capacitor
and indicate if this result shows the capacitor is time invariant.
(c) An amplitude modulation system has an input/output equation
Answer: (a) z(t) solution of differential equation with initial condition of 2; (b) If i(t) = 0 for
t < 0, system is TI; (c) If x(t − 0.5) = u(t − 0.5) the output is y1 (t) = sin(2πt)u(t − 0.5).
Solution
(a) Derivative
dz(t)
= w(t) − w(t − 1)
dt
which excludes the initial condition of 2. System is LTI if initial condition is zero.
that is, provided that i(t) = 0 for t < 0, the system is time-invariant.
Rt
ii. If i(t) = u(t) then vc (t) = 0
u(τ )dτ = r(t). If we shift the inputs i1 (t) = i(t − 1) =
u(t − 1) the previous output is shifted, so system is time-invariant.
(c) If x(t) = u(t) then y(t) = sin(2πt)u(t) while corresponding to x(t − 0.5) = u(t − 0.5) is
y1 (t) = sin(2πt)u(t−0.5) indicating the system is not time-invariant as y1 (t) is not y(t−0.5).
sin(2 π t) heaviside(t)
2.3 The following problems relate to linearity, time–invariance and causality of systems.
(a) A system is represented by the equation z(t) = v(t)f (t) + B where v(t) is the input, z(t)
the output, f (t) a function and B a constant.
where T > 0 is the averaging interval, x(t) and y(t) are the system input and output.
Answer: (a) If f (t) = v(t) = u(t) − u(t − 1) then z(t) = u(t) − u(t − 1), and if input is
v1 (t) = u(t−2)−u(t−3) output is z1 (t) = v1 (t)f (t) = 0; (b) If x(t) = u(t) then y(t) = r(t)−r(t−1)
and if x1 (t) = u(t − 2) the corresponding output is y1 (t) = y(t − 2).
Solution
ii. If x(t) = u(t) then y(t) = r(t) − r(t − 1) and if x1 (t) = u(t − 2) the corresponding
output is y1 (t) = y(t − 2). System is time-invariant.
In general, if x1 (t) = x(t − λ) the output is
Z t Z t−λ
1 1
x1 (τ ) dτ = x(ν)dν
T t−T | {z } T t−T −λ
x(t−λ)
(a) Find the impulse response h(t) of the averager. Is this system causal?
(b) Let x(t) = u(t), find the output of the averager.
Solution
(a) Letting x(t) = δ(t) the impulse response is
Z t+T /2
1
h(t) = δ(τ )dτ
T t−T /2
Z t Z t+T /2
1 1
= δ(τ )dτ + δ(τ )dτ
T t−T /2 T t
If t > 0, and t − T /2 < 0 the first integral includes 0, while the second does not. Thus
Z
1 t 1
h(t) = δ(τ )dτ + 0 = t > 0 and t − T /2 < 0, or 0 < t < T /2
T t−T /2 T
Likewise when t < 0 then t − T /2 < −T /2 and t + T /2 < T /2 the reverse of the previous case
happens and so
Z t+T /2
1 1
h(t) = 0 + δ(τ )dτ = t < 0 and t + T /2 > 0, or − T /2 < t < 0
T t T
so that
1
[u(t + T /2) − u(t − T /2)]
h(t) =
T
indicating that the system is non-causal as h(t) 6= 0 for t < 0.
(b) If x(t) = u(t) then the output of the averager is
Z t+T /2
1
y(t) = u(τ )dτ
T t−T /2
If t + T /2 < 0 then y(t) = 0 since the argument of the unit step signal is negative. If t + T /2 ≥ 0
and t − T /2 < 0 then
Z t+T /2
1
y(t) = u(τ )dτ = (t + T /2)
0 T
and finally when t − T /2 ≥ 0 then
Z t+T /2
1
y(t) = u(τ )dτ = 1
T t−T /2
2.5 The input-output equation characterizing a system of input x(t) and output y(t) is
Z t
−2t
y(t) = e y(0) + 2 e−2(t−τ ) x(τ )dτ t≥0
0
(a) Find the ordinary differential equation that also characterizes this system.
(b) Suppose x(t) = u(t) and any value of y(0),we wish to determine the steady state response
of the system/ Is the value of y(0) of any significance, i.e., do we get the same steady-state
response if y(0) = 0 or y(0) = 1? Explain.
(c) Compute the steady-state response when y(0) = 0, and x(t) = u(t). To do so first find the
impulse response of the system, h(t) using the given input–output equation, and then find
y(t) by computing the convolution graphically to determine the steady–state response.
(d) Suppose the input is zero, is the system depending on the initial condition BIBO stable?
Find the zero-input response y(t) when y(0) = 1. Is it bounded?
Solution
(a) To find the differential equation let y(0) = 0, so that
Z t
y(t) = 2e−2t e2τ x(τ )dτ
0
e2τ t
y(t) = y(0)e−2t + 2e−2t = y(0)e−t + (1 − e−2t )
2 0
so that in the steady state, i.e., when t → ∞, the output is 1, independent of the value of the
initial condition.
(c) From the given input-output equation, letting x(t) = δ(t) and y(0) = 0 the output is the
impulse response of the system
Z t Z t
h(t) = 2 e−2(t−τ ) δ(τ )dτ = 2e−2t δ(τ )dτ = 2e−2t u(t)
0 0
it can be obtained graphically by reflecting the input x(t) = u(t) and shifting it linearly from
−∞ to ∞ to get
and it is bounded for any finite value of the initial condition y(0), in particular for y(0) = 1,
therefore the system that depends on the initial condition is BIBO.
2.6 The input x(t) and the corresponding output y(t) of a linear time-invariant (LTI) system are
(i) x1 (t) = u(t) − u(t − 1) − u(t − 2) + u(t − 3), (ii) x2 (t) = u(t + 1) − 2u(t) + u(t − 1),
(iii) x3 (t) = δ(t) − δ(t − 1)
Solution
(a) x1 (t) = x(t)−x(t−2) so y1 (t) = y(t)−y(t−2), two triangular pulses, the second multiplied
by −1.
x1 (t) y1 (t)
1
1
2 3 t
2 4
t
1
1
x2 (t) y2 (t)
1
1
1 2
t t
1 1
1
1
(b) x2 (t) = x(t + 1) − x(t) then y2 (t) = y(t + 1) − y(t) (they overlap between 0 and 1).
(c) x3 (t) = δ(t) − δ(t − 1) so y3 (t) = dy(t)/dt = u(t) − 2u(t − 1) + u(t − 2). Considering that
the output of x(t) is y(t), i.e., y(t) = S[x(t)], and that the integrator and the differentiator
are LTI systems Fig. 2.3 shows how to visualize the result in this problem by considering
that you can change the order of the cascading of LTI systems.
dy(t) dy(t)
Z x(t) y(t) dt Z dt
d d
dx(t)
S
dt
⌘ S
dt
x3 (t) = dx(t) dy(t)
dt x3 (t) = dt
dt | {z }
1
(c) Another way to show the above expression is the solution of the ordinary differential
equation is to multiply both sides of the ordinary differential equation by eat
i. show that the left term is d(eat y(t))/dt,
ii. integrate the two terms to obtain the solution.
Solution
f (t)h
= lim = f (t)
h→0 h
(b) y(t) satisfies the ordinary differential equation, indeed
Z
dy(t) −at d −at t aτ
= −ay(0)e + e e x(τ )dτ
dt dt 0
Z t
= −a[y(0)e−at + e−at eaτ x(τ )dτ ] + x(t) = −ay(t) + x(t)
0
| {z }
y(t)
(c) Multiplying the two terms of the ordinary differential equation by eat we get
2.8 The bounded-input bounded-output stability assumes that the input is always bounded, lim-
ited in amplitude. If that is not the case, even a stable system would provide an unbounded
output. Consider the analog averager, with an input-output relationship
Z t
1
y(t) = x(τ )dτ
T t−T
(a) Suppose that the input to the averager is a bounded signal x(t), i.e., there is a finite value
M such that |x(t)| < M . Find the value for the bound of the output y(t) and determine
whether the averager is BIBO stable or not.
(b) Let the input to the averager be x(t) = tu(t), i.e., a ramp signal, compute the output y(t)
and determine if it is bounded or not. If y(t) is not bounded, does that mean that the
averager is an unstable system? Explain.
Answers: h(t) = (1/T )[u(t) − u(t − T )]; if input is ramp signal the output is not bounded.
Solution
(a) We can either find the impulse response or show that the output is bounded for any input
signal that is bounded. To find the impulse response change the variable to σ = t − τ which
gives
Z 0 Z T
−1 1
y(t) = x(t − σ)dσ = x(t − σ)dσ
T T 0 T
which is a convolution integral with h(t) = (1/T )[u(t) − u(t − T )]. This impulse response is
absolutely integrable as
Z ∞ Z T
1
|h(t)|dt = dt = 1
0 0 T
Likewise, if we assume that x(t) is bounded, i.e., there is a value M < ∞ such that |x(t)| < M ,
then Z Z
t t
1 M
|y(t)| ≤ |x(τ )|dτ ≤ dτ = M < ∞
T t−T T t−T
2.9 An ideal low-pass filter — The impulse response of an ideal low-pass filter is h(t) = sin(t)/t.
(a) Given that the impulse response is the response of the system to an input x(t) = δ(t)
with zero initial conditions, can an ideal low-pass filter be used for real-time processing?
Explain.
(b) Is the ideal low-pass filtering bounded-input bounded-output stable? Use MATLAB to
check if the impulse response satisfies the condition for BIBO stability.
Solution
(a) The system with the sinc impulse response is non-causal because h(t) 6= 0 for t < 0. As
such an ideal low-pass filter cannot be used for real-time processing as it would require future
inputs, not available in real-time. This can be seen by looking at the convolution with a causal
signal x(t):
Z ∞
y(t) = x(τ )h(t − τ )dτ
0
Z t Z ∞
= x(τ )h(t − τ )dτ + x(τ )h(t − τ )dτ
0 t
where the upper limit of the top integral is due to h(τ ) having an infinite support, and the
lower limit because x(t) is causal. The second bottom integral shows that the output depends
on future values of the input.
abs(sin(t))/abs(t)
1
0.8
0.6
0.4
0.2
0
−30 −20 −10 0 10 20 30
t
6
Integral
0
0 5 10 15 20 25 30
n
Figure 2.4: Problem 9: absolute value of sinc function and its integral values for t = n sec.
(b) For the low-pass filter to be BIBO stable the impulse response h(t) must be absolutely inte-
grable, i.e., Z ∞
| sin(t)/t|dt < ∞
−∞
To approximtely compute this integral we use the following script. The integral is bounded, so
the filter is BIBO stable.
%% Pr 2.9
clear all; clf
syms x t x1
x=abs(sinc(t/pi));
for k=1:30,
x1=int(x,t,0,k);
xx(k)=subs(2*x1);
end
n=1:30;
figure(1)
subplot(211)
ezplot(x,[-30,30]); axis([-30 30 -0.1 1.2]);grid
subplot(212)
stem(n,xx); axis([1 30 0 8]);grid
Answers: (a) Y (s) = 2 sinh(s), ROC the whole s-plane; W (s) = 2s sinh(s)/(s2 + 4π 2 ), ROC the
whole s-plane; (b) Y1 (s) = e−(s+1) /(s + 1), ROC σ > −1.
Solution
1
Chaparro-Akan — Signals and Systems using MATLAB 3.2
iv. Notice that because cos(2π(t ± 1)) = cos(2πt ± 2π) = cos(2πt) then
so that
es s e−s s (es − e−s )s 2s sinh(s)
W (s) = − = = 2
s2 + 4π 2 s2 + 4π 2 s2 + 4π 2 s + 4π 2
The poles p1,2 = ±j2π are cancelled by zeros z1,2 = ±j2π ( indeed, es − e−s |s=±j2π =
1 − 1 = 0) so the ROC is the whole s-plane.
(b) Causal signals.
i. X1 (s) = L[x1 (t)] = 1/(s + 1), ROC σ > −1
ii. Notice the exponential is not delayed, so the shifting property cannot be applied.
Instead,
e−(s+1)
Y1 (s) = L[y1 (t)] = L[e−1 e−(t−1) u(t − 1)] = , ROC σ > −1
s+1
iii. L[z1 (t)] = L[e−(t−1) u(t − 1)] = e−s /(s + 1), ROC σ > −1
iv. Since w1 (t) = x1 (t) − y1 (t), where x1 (t) and y1 (t) are given before then
1 − e−(s+1)
L[w1 (t)] = X1 (s) − Y1 (s) = ,
s+1
ROC whole s-plane because pole s = −1 is cancelled by zero, i.e., 1−e−(s+1) |s=−1 = 0.
3.2 Consider the pulse x(t) = u(t) − u(t − 1). Find the zeros and poles of X(s) and plot them.
(a) Suppose x(t) is the input of a LTI system with a transfer function H(s) = Y (s)/X(s) =
1/(s2 + 4π 2 ), find and plot the poles and zeros of Y (s) = L[y(t)]), where y(t) is the output
of the system.
(b) If the transfer function of the LTI system is
∞
Y
Z(s) 1
G(s) = =
X(s) s + (2kπ)2
2
k=1
and the input is the above signal x(t), compute the output z(t).
Solution
The transform is X(s) = (1 − e−s )/s. The pole of X(s) is s = 0, while the zeros of X(s) are
values of s that make 1 − e−s = 0 or sk = j2kπ for −∞ < k < ∞ and integer. The pole is
cancelled with the k = 0 zero, so X(s) has an infinite number of zeros on the jΩ-axis.
(a) The output of the LTI system has Laplace transform
∞
(s + j2π)(s − j2π) Y
Y (s) = X(s)H(s) = (s + j2πk)(s − j2πk)
s2 + 4π 2
k=2
∞
Y
= (s2 + (2πk)2 )
k=2
There are no poles, but an infinite number of zeros at ±j2πk for integer k ≥ 2.
(b) This is similar to the above problem, but here
Z(s) = X(s)G(s) = 1
which gives
z(t) = δ(t)
3.3 Find the Laplace transform of the following signals and their region of convergence
(a) the reflection of the unit step signal, i.e., u(−t). And then use the result together with the
Laplace transform of u(t) to see if you can obtain the Laplace transform of a constant or
x(t) = u(t) + u(−t) (assume u(0) = 0.5 so there is no discontinuity at t = 0).
(b) the non-causal signal y(t) = e−|t| u(t + 1). Carefully plot it, and find its Laplace transform
Y (s) by separating y(t) into a causal and a non-causal signals, yc (t) and yac (t), and plot
them separately. Find the ROC of Y (s), Yc (s) and Yac (s).
Answers: (a) x(t) has no Laplace transform; (b) Y (s) = (es−1 (s + 1) − 2)/(s2 − 1).
Solution
1 es−1 − 1
Y (s) = Yc (s) + Yac (s) = +
s+1 s−1
ses−1 + es−1 − 2
= ROC σ > −1
s2 − 1
The poles of Y (s) are s = ±1, and a zero is s = 1. Thus the zero cancels the pole at s = 1
leaving the pole at s = −1, and so the ROC is σ > −1.
3.4 You are given the following Laplace transform of the output y(t) of a system with input x(t)
and Laplace transform X(s):
(s − 1)X(s) 1
Y (s) = +
(s + 2)2 + 1 (s + 2)2 + 1
Solution
(s − 1)X(s)
(s + 2)2 + 1
1
.
(s + 2)2 + 1
s−1 A B + Cs
Yzs (s) = 2
= +
s((s + 2) + 1) s (s + 2)2 + 1
−1/5 9/5 + s/5 −1/5 1 (s + 2) + 7
= + 2
= +
s (s + 2) + 1 s 5 (s + 2)2 + 1
−1/5 1 (s + 2) 7 1
= + +
s 5 (s + 2)2 + 1 5 (s + 2)2 + 1
1
Yzi (s) = ⇒ yzi (t) = e−2t sin(t)u(t)
(s + 2)2 + 1
(a) Use the Laplace transform to find the unit-step response s(t) = (h ∗ x)(t).
(b) Find the response due to the following inputs
(i) x1 (t) = u(t) − u(t − 1), (ii) x2 (t) = δ(t) − δ(t − 1), (iii) x3 (t) = r(t)
Solution
(a) The Laplace transform of s(t) = (h ∗ x)(t) is S(s) = H(s)X(s) which gives
1
S(s) = ⇒ s(t) = e−t sin(t)u(t)
(s + 1)2 + 1
(b) By LTI
3.6 Consider the impulse response of a LTI system h(t) = e−at [u(t) − u(t − 1)] a > 0.
Answers: H(s) = (1 − e−(s+a) )/(s + a), ROC: whole s plane, pole cancelled by zero.
Solution
Y (s) 1
H(s) = =
X(s) 1 − e−s
(a) Find the poles and zeros of H(s), and from this determine if the filter is BIBO stable or not.
(b) Draw a block diagram for such a system.
(c) Find the impulse response h(t) of the system and use it to verify your stability result.
Solution
(a) Poles satisfy the identity e−s = 1 = ej2πk for k = 0, ±1, ..., so there are an infinite number
of poles, sk = −j2πk. No zeros. The system is unstable because of poles on jΩ-axis.
x(t) + y(t)
Delay
y(t − 1)
(b) Y (s)(1 − e−s ) = X(s) gives y(t) = y(t − 1) + x(t). A positive feedback system with a delay
in the feedback represents the system. See Fig. 3.1.
(c) The transfer function is
X∞
1
H(s) = = (e−s )n
1 − e−s n=0
with inverse
∞
X
h(t) = δ(t − n)
n=0
(a) The impulse response of a LTI system is h(t) = e−2t u(t) and the system input is a pulse
x(t) = u(t) − u(t − 3). Find the output of the system y(t) by means of the convolution
integral graphically and by means of the Laplace transform.
(b) It is known that the impulse response of an analog averager is h1 (t) = u(t) − u(t − 1),
consider the input to the averager x1 (t) = u(t) − u(t − 1), determine graphically as well
as by means of the Laplace transform the corresponding output of the averager y1 (t) =
[h1 ∗ x1 ](t). Is y1 (t) smoother than the input signal x1 (t)? provide an argument for your
answer.
(c) Suppose we cascade 3 analog averagers each with the same impulse response hi (t) =
u(t) − u(t − 1), i = 1, 2, 3, determine the transfer function of this system. If the duration of
the input to the first averager is M sec., what would be the duration of the output of the
3th averager?
Answers: (a) Y (s) = (1 − e−3s )/(s(s + 2)); (c) the support of the output will be 3M .
Solution
(a) Using the Laplace transform
1
H(s) =
s+2
1 − e−3s
X(s) =
s
1 − e−3s
Y (s) = H(s)X(s) =
s(s + 2)
Letting
1 0.5 0.5
F (s) = = − ⇒ f (t) = 0.5(1 − e−2t )u(t)
s(s + 2) s s+2
then
y(t) = f (t) − f (t − 3) = 0.5(1 − e−2t )u(t) − 0.5(1 − e−2(t−3) )u(t − 3)
Graphically, we plot h(τ ) and x(t − τ ) as functions of τ and shift x(t − τ ) from the left to the
right and integrate the overlapping areas to get
0 t<0
R
t −2τ
y(t) = e dτ = 0.5(1 − e−2t ) 0≤t<3
0
Rt
t−3
e−2τ dτ = 0.5(e−2(t−3) − e−2t ) t ≥ 3
1 − e−s
H1 (s) = X1 (s) = H(s)
s
1 − 2e−s + e−2s
Y1 (s) =
s2
so that
y1 (t) = r(t) − 2r(t − 1) + r(t − 2)
a triangular pulse.
Graphically, letting x1 (τ ) and h1 (t − τ ) (moving from left to right) be functions of τ , and inte-
grating their overlap gives
0 t<0
Rt
0 dτ = t 0≤t<1
y1 (t) =
R1
dτ = 1 − t + 1 = 2 − t 1 ≤ t < 2
t−1
0 t≥2
coinciding with the result using the Laplace transform. The signal y1 (t) is smoother than x1 (t)
since it is the output of an averager filter.
(c) Using the Laplace transform, letting h(t) = hi (t), i = 1, 2, 3
1 − e−s
H(s) =
s
(1 − e−s )3 1 − 3e−s + 3e−2s − e−3s
(H(s))3 = =
s3 s3
The support of the output will be 3M .
3.9 To see the effect of the zeros on the complete response of a system, suppose you have a system
with a transfer function
Y (s) s2 + 4
H(s) = =
X(s) s((s + 1)2 + 1)
(a) Find and plot the poles and zeros of H(s). Is this BIBO stable?
(b) Find the frequency Ω0 of the input x(t) = 2 cos(Ω0 t)u(t) such that the output of the given
system is zero in the steady state. Why do you think this happens?
(c) If the input is x(t) = 2 sin(Ω0 t)u(t) would you get the same result as above? Explain why
or why not.
Answers: h(t) not absolutely integrable; cosine as input gives Y (s) = 2/((s + 1)2 + 1).
Solution
(a) The poles of H(s) are at s = 0 and s = −1 ± j1, and the zeros are at s = ±j2. This system
is not BIBO stable because of the pole at the origin, the impulse response will not be absolutely
integrable.
(b) If x(t) = 2 cos(2t)u(t) then
s
X(s) = 2
s2 + 4
and thus
2
Y (s) =
(s + 1)2 + 1
i.e., the pole at s = 0 and the zeros are cancelled by the transform of the cosine. The remaining
poles are in the left-hand s-plane and so the steady state is zero.
This is a very special case because of the cancellation of the pole at s = 0 and the numerator
due to the used input. The system is not stable, and the frequency response does not exist for
all frequencies. Indeed it is infinite whenever the input frequency is 0.
(c) The input cannot be a sine of the same frequency because its Laplace transform does not
cancel the pole at zero and so the steady state is not zero.
3.10 Inverse Laplace transform— Consider the following inverse Laplace transform problems.
Answers: y(t) = −4δ(t) + dδ(t)/dt + · · · ; z(t) = δ(t) − 2e−t u(t) − e−2t u(t).
% Pr. 3_10
clear all; clf
syms s t w
disp(’>>>>> Inverse Laplace <<<<<’)
x=ilaplace((3*s-4)/(sˆ3+3*sˆ2+2*s))
figure(1)
ezplot(x,[0,10])
axis([0 10 -2.5 .5])
grid
−0.5
−1
−1.5
−2
−2.5
0 1 2 3 4 5 6 7 8 9 10
t
3.11 Solving differential equations— One of the uses of the Laplace transform is the solution of
differential equations.
(a) Suppose you are given the ordinary differential equation that represents a LTI system,
where y(t) is the output and x(t) is the input of the system, y (1) (t) and y (2) (t) are first
and second order derivatives with respect to t. The input is causal, i.e., x(t) = 0 t < 0.
What should the initial conditions be for the system to be LTI? Find Y (s) for those initial
conditions.
(b) If y (1) (0) = 1 and y(0) = 1 are the initial conditions for the above ordinary differential
equation, find Y (s). If the input to the system is doubled, i.e., the input is 2x(t) is Y (s)
doubled so that its inverse Laplace transform y(t) is doubled? Is the system linear?
(c) Use MATLAB to find the solutions of the ordinary differential equation when the input is
u(t) and 2u(t) with the initial conditions given above. Compare the solutions and verify
your response in (b).
Answers: (b) Y (s) = X(s)/(s2 + 0.5s + 0.15) + (s + 1.5)/(s2 + 0.5s + 0.15), not LTI.
Solution
(a) If the initial conditions are zero the system is LTI. The Laplace transform of the differential
equations is then
Y (s)[s2 + 0.5s + 0.15] = X(s)
and so
X(s)
Y (s) =
s2 + 0.5s + 0.15
(b) If the initial conditions are dy(0)/dt = 1 and y(0) = 1 the Laplace transform of the differen-
tial equation is
[Y (s)s2 − s − 1] + 0.5[sY (s) − 1] + 0.15Y (s) = X(s)
which gives
X(s) s + 1.5
Y (s) = + 2
s2 + 0.5s + 0.15 s + 0.5s + 0.15
which can be written, with the appropriate definitions, as Y (s) = X(s)/A(s) + I(s)/A(s), so
that when we double the input the Laplace transform of the corresponding output is
2X(s) I(s)
+ 6= 2Y (s)
A(s) A(s)
s2 + 1.5s + 1
Y (s) =
s3 + 0.5s2 + 0.15s
when we double the input the Laplace transform of the output (again including the initial
conditions) is
s2 + 1.5s + 2
Y1 (s) = .
s3 + 0.5s2 + 0.15s
The following script is used to obtain the two responses. The second one is not the double of
the first, so the system is not LTI — due to the initial conditions not being zero.
%% Pr. 3_11
clear all; clf
% response to u(t)
den=[1 0.5 0.15 0];
num=[0 1 1.5 1];
syms s t y y1
[p,r]=pfeLaplace(num,den)
disp(’>>>>> Inverse Laplace <<<<<’)
y=ilaplace((num(2)*sˆ2+num(3)*s+num(4))/(den(1)*sˆ3 +den(2)*sˆ2+den(3)*s))
% response to 2u(t)
num=[0 1 1.5 2];
pfeLaplace(num,den)
disp(’>>>>> Inverse Laplace <<<<<’)
y1=ilaplace((num(2)*sˆ2+num(3)*s+num(4))/(den(1)*sˆ3 +den(2)*sˆ2+den(3)*s))
figure(1)
subplot(311)
ezplot(y,[0,40])
axis([0 40 -0.01 11]);grid;title(’Response to u(t)’)
subplot(312)
ezplot(y1,[0,40])
axis([0 40 -0.01 15]);grid;title(’Response to 2 u(t)’)
subplot(313)
ezplot(y1-2*y,[0,40]);grid;title(’Difference of two responses’)
axis([0 40 -10 10])
Response to u(t)
10
0
0 5 10 15 20 25 30 35 40
t
Response to 2 u(t)
15
10
0
0 5 10 15 20 25 30 35 40
t
Difference of two responses
10
−10
0 5 10 15 20 25 30 35 40
t
Figure 3.3: Problem 11: Responses y(t) and y1 (t) corresponding to inputs u(t) and 2u(t), the differ-
ence y1 (t) − 2y(t) 6= 0 indicates the system is not linear.
3.12 Zero steady-state response of analog averager— An analog averager can be represented by
the differential equation
dy(t) 1
= [x(t) − x(t − T )]
dt T
where y(t) is its output and x(t) the input.
Show how to obtain the above differential equation and that y(t) is the solution of the
differential equation.
(b) If x(t) = cos(πt)u(t), choose the value of T in the averager so that the output y(t) = 0 in
the steady state. Graphically show how this is possible for your choice of T . Is there a
unique value for T that makes this possible? How does it relate to the frequency Ω0 = π
of the sinusoid?
(c) Use the impulse response h(t) of the averager found before, to show using Laplace that
the steady state is zero when x(t) = cos(πt)u(t) and T is the above chosen value. Use
MATLAB to solve the differential equation and to plot the response for the value of T you
chose. [HINT: consider x(t)/T the input and use superposition and time-invariance to
find y(t) due to (x(t) − x(t − T ))/T .]
Solution
(a) The input/output equation for the analog averager is
Z
1 t
y(t) = x(τ )dτ
T t−T
1 − e−2s
Y (s) = 0.5 = Y1 (s) − Y1 (s)e−2s
s2 + π 2
where Y1 (s) = 0.5/(s2 + π 2 ) so that
0.5
y(t) = (sin(πt)u(t) − sin(π(t − 2)u(t − 2))
π
%% Prob 3_12
clear all; clf
T=2;
num=[0 0 1/T]
den=[1 0 piˆ2]
syms s t y
figure(1)
subplot(221)
[r,p]=pfeLaplace(num,den);
disp(’>>>>> Inverse Laplace <<<<<’)
y1=ilaplace(0.5/(sˆ2+piˆ2))
subplot(222)
ezplot(y1,[0,8])
axis([0 8 -0.2 0.2]); title(’y_1(t)’)
grid
y=y1-y1*heaviside(t-2)
subplot(223)
ezplot(y,[0,8])
y1(t)
4 0.2
2 0.1
0 0
jΩ
−2 −0.1
−4 −0.2
−1 0 1 2 3 0 2 4 6 8
σ t
y(t)
0.2
0.1
−0.1
−0.2
0 2 4 6 8
t
Figure 3.4: Problem 12: Poles/zeros of Y1 (s) (top–left) , y1 (t) and y(t) (bottom).
3.13 Polynomial multiplication— When the numerator or denominator are given in a factorized
form, we need to multiply polynomials. Although this can be done by hand, MATLAB provides
the function conv that computes the coefficients of the polynomial resulting from the product
of two polynomials.
(a) Use help in MATLAB to find how conv can be used, and then consider two polynomials
do the multiplication of these polynomials by hand to find Z(s) = P (s)Q(s) and use conv
to verify your results.
(b) The output of a system has a Laplace transform
N (s) (s + 2)
Y (s) = = 2
D(s) s (s + 1)((s + 4)2 ) + 9)
use conv to find the denominator polynomial and then find the inverse Laplace transform
using ilaplace.
Solution
(a)(b) The following script shows how to do polynomial multiplication using conv function.
The multiplication of the polynomials
gives
Z(s) = P (s)Q(s) = 2s5 + 5s4 + 6s3 + 5s2 + 2s + 1
The script shows how to obtain this result in the Poly Multiplication part of the script.
In the Application part of the script we show how to find the numerator N1 (s) and the denom-
inator D1 (s) of Y (s) by multiplying X(s) and H(s) to give
X(s)N (s) (s + 2)
Y (s) = = 2
D(s) s (s + 1)((s + 4)2 + 9)
so that the poles of Y (s) are s = 0 (double), and s = −4 ± j3, and no zeros. The double pole
gives a ramp and the other poles a modulated sinusoid by a decaying exponential.
% Pr. 3_13
% Poly multiplication
% coefficients of higher to lower orders of s
P=[1 1 1];Q=[2 3 1 1];
Z=conv(P,Q);
% Application
d1=[1 1]; d2=[1 8 25]; D=conv(d1,d2);
d3=[1 0 0]; den=conv(D,d3)
N=length(den)
num=[ zeros(1,N-2) 1 2]
syms s t y
figure(1)
subplot(121)
[r,p]=pfeLaplace(num,den);
disp(’>>>>> Inverse Laplace <<<<<’)
y=ilaplace((num(N-1)*s+num(N))/(den(1)*sˆ5+den(2)*sˆ4+den(3)*sˆ3+den(4)*sˆ2+den(5)*s+d
subplot(122)
ezplot(y,[0,50])
axis([0 50 -0.1 1.5]); grid
den = 1 9 33 25 0 0
num = 0 0 0 0 1 2
>>>>> Zeros <<<<<
z = -2
>>>>> Poles <<<<<
p = -4.0000 + 3.0000i
-4.0000 - 3.0000i
-1.0000
0
0
>>>>> Residues <<<<<
r = 0.0050 - 0.0026i
0.0050 + 0.0026i
0.0556
-0.0656
0.0800
>>>>> Inverse Laplace <<<<<
y = -41/625+1/18*exp(-t)+2/25*t+113/11250*exp(-4*t)*cos(3*t)
+59/11250*exp(-4*t)*sin(3*t)
3 4
3.5
2
3
1
2.5
0
jΩ
−1
1.5
−2 1
−3 0.5
0
−4
−4 −2 0 2 0 10 20 30 40 50
σ t
4.1 Consider the following problems related to periodicity and Fourier series:
i. Determine which of these signals is periodic and the corresponding fundamental pe-
riod.
ii. For the periodic signal, find a trigonometric Fourier series representation.
iii. Indicate why it is not possible to obtain a trigonometric Fourier series representation
for the non-periodic signal.
(b) Consider the signal
sin(2t + π)
x(t) =
sin(t)
i. Is this signal periodic? If so, what is its fundamental period T0 ?
ii. Find the Fourier series of x(t).
Answers: (a) x1 (t) periodic of fundamental period T = 1; x2 (t) does not have Fourier series;
(b) periodic of fundamental period T0 = 2π; x(t) = −2 cos(t).
Solution
1
Chaparro-Akan — Signals and Systems using MATLAB 4.2
(a) i. The sinusoidal components of x1 (t) have periods T1 = 1 and T2 = 1/3, with ratio
T1 /T2 = 3 so that x1 (t) is periodic of fundamental period T1 = 3T2 = 1.
The periods of the sinusoidal components of x2 (t) are T1 = 1 and T2 = 2π/6, their
ratio is T1 /T2 = 3/π, not a rational number so x2 (t) is not periodic.
ii. Trigonometric Fourier series of x1 (t), with fundamental period T0 = 1 and fundamen-
tal frequency Ω0 = 2π, is
sin(2(t + T0 ) + π)
x(t + T0 ) = = x(t)
sin(t + T0 )
this is possible if T0 = 2π
ii. It can be shown that x(t) = −2 cos(t), indeed sin(2t + π) = − sin(2t) and
so
x(t) = −2 cos(t) = −ejt − e−jt
4.2 Find the complex exponential Fourier series for the following signals. In each case plot the
magnitude and phase line spectra for k ≥ 0.
(i) x1 (t) = cos(5t + 45o ), (ii) x2 (t) = sin2 (t), (iii) x3 (t) = cos(3t) + cos(5t)
∗
Answers: The Fourier series coefficients for x1 (t) are X1 = X−1 = 0.5ejπ/4 .
Solution
(i) Fundamental frequency Ω0 = 5 so
∗
and X1 = X−1 = 0.5ejπ/4 .
(ii) We have
2
1 jt
x2 (t) = (e − e−jt ) = 0.5 − 0.25ej2t − 0.25e−j2t
2j
∗
and so X0 = 0.5, X1 = X−1 = −0.25 and Ω0 = 2 is the fundamental frequency.
(iii) Fundamental frequency of x3 (t) is Ω0 = 1 (ratio of fundamental periods of sinusoidal
components T1 /T2 = (2π/3)/(2π/5) = 5/3 and Ω0 = 2π/3T1 = 1).
0.4 0.6
|Xk|
∠Xk
0.4
0.2
0.2
0 0
0 2 4 6 0 2 4 6
x2(t)
0.6 3
0.4 2
|Xk|
∠Xk
0.2 1
0 0
0 2 4 6 8 10 0 2 4 6 8 10
x3(t)
0.6 3
0.4 2
|Xk|
∠Xk
0.2 1
0 0
0 1 2 3 4 5 0 1 2 3 4 5
kΩ0 kΩ0
4.3 Let the Fourier series coefficients of a periodic signal x(t) of fundamental frequency Ω0 =
2π/T0 ) be {Xk }. Consider the following functions of x(t).
Determine if y(t), z(t) and w(t) are periodic, and if so give the corresponding Fourier coeffi-
cients in terms of those for x(t).
Answers: Zk = Xk (1 + e−2jΩ0 k ); Wk = Xk/2 for k even and 0 otherwise.
Solution
(a) If the periodic signal
X
x(t) = Xk ejΩ0 kt
k
then
X
y(t) = 2x(t) − 3 = (2X0 − 3) + 2Xk ejΩ0 kt
k6=0
The signal
4.4 The period of a periodic signal x(t) with fundamental period T0 = 2 is x1 (t) = cos(t)[u(t) −
u(t − 2)].
(a) Plot the signal x(t) and find the Fourier series coefficients Xk for x(t) using the integral
equation.
(b) Use the Laplace transform to find the Fourier series coefficients {Xk }.
(c) After considering the symmetry of the zero-mean signal x(t) − X0 , can you determine
whether the Xk should be real, purely imaginary or complex? Indicate how you make
this determination.
(d) To simplify computing the Laplace transform of x1 (t) consider the following expression:
Find the values of A and B and show how to compute L[cos(t)u(t − 2)].
Solution
x(t)
1
··· ···
2 4 6
t
cos(2)
Z 2 2
1 1 jt 1 ejt(1−kπ) e−jt(1+kπ)
Xk = (e + e−jt )e−jkπt dt = +
2 0 2 4 j(1 − kπ) −j(1 + kπ) 0
j2 −j2
1 e −1 e −1
= −
4 j(1 − kπ) j(1 + kπ)
1 j2
= 2 2
j(e − e−j2 ) + jπk(ej2 + e−j2 ) − j(1 + kπ) + j(1 − kπ)
4(k π − 1)
jkπ(1 − cos(2)) + sin(2)
=
2(1 − k 2 π 2 )
Z ∞ Z ∞
L[cos(t)u(t − 2)] = cos(t)e−st dt = cos(ρ + 2)e−s(ρ+2) dρ
2 0
Z ∞ Z ∞
−2s
= e cos(2) cos(ρ)e−sρ dρ − sin(2) sin(ρ)e−sρ dρ
0 0
= e−2s [cos(2)L[cos(t)u(t)] − sin(2)L[sin(t)u(t)]]
cos(2)s sin(2)
= e−2s −
s2 + 1 s2 + 1
We then have
s e−2s (cos(2)s − sin(2))
Xk = 0.5 2
−
s +1 s2 + 1 s=jkπ
jkπ(1 − cos(2)) + sin(2)
=
2(1 − k 2 π 2 )
which gives
p
A2 + B 2 = 1
tan−1 (B/A) = −2
4.5 Consider the following problems related to steady state and frequency responses.
(a) The input x(t) and the output y(t), −∞ < t < ∞, of a filter with transfer function H(s) are
Determine as much as possible about the frequency response H(jΩ) of the filter using the
above input and output.
(b) Given a LTI system with frequency response
(
−π/2 Ω≥0
|H(jΩ)| = u(Ω + 2) − u(Ω − 2), ∠H(jΩ) =
π/2 Ω<0
Answers: (a) H(j2π) = −0.5, H(j3π) = 0; (b) yss (t) = 2 cos(3t/2 − π/2).
Solution
so H(j2π) = 0.5ejπ = −0.5, H(j3π) = 0. Nothing else can be learned about the filter from
the input/output.
(b)
∞
X 2
yss (t) = |H(j3k/2)| cos(3kt/2 + ∠H(j3k/2))
k2
k=1
= 2|H(j3/2)| cos(3t/2 + ∠H(j3/2))
= 2 cos(3t/2 − π/2)
(a) A periodic signal x(t) of fundamental frequency Ω0 = π/4 is the input of an ideal band-
pass filter with the following frequency response
1 π ≤ Ω ≤ 3π/2 −Ω
π ≤ Ω ≤ 3π/2
|H(jΩ)| = 1 −3π/2 ≤ Ω ≤ −π ∠[H(jΩ)] = Ω −3π/2 ≤ Ω ≤ −π
0 otherwise 0 otherwise
∗ ∗
The non-zero Fourier series coefficients of x(t) are X1 = X−1 = j, X5 = X−5 = 2.
i. Express x(t) in the form
∞
X
x(t) = Ak cos(Ωk + θk ).
k=0
i. Carefully plot the frequency response, magnitude and phase, of the filter and deter-
mine the type of filter it is.
ii. Calculate the steady-state response yss (t) of the filter for the given input.
Answers: (a) y(t) = 4 cos(5π(t − 1)/4); (b) yss (t) = 2 cos(3πt − π/2).
Solution
(a)
|H(jΩ)|
1
Ω
−4.5π −1.5π 1.5π 4.5π
∠H(jΩ)
π/2
Ω
−π/2
ii. According to the frequency response of the band-pass filter y(t) can have a second,
third and fourth harmonics, but since the magnitude of the second and fourth Fourier
coefficients are X2 = X4 = 0, the output has only the third harmonic giving
|H(jΩ)| ∠H(jΩ)
2
π/2
Ω Ω
−1.5π 1.5π
−π/2
(a) Determine the Fourier series coefficients needed to find the output y(t) of the filter.
(b) Is the output signal y(t) periodic? If so, determine its fundamental period T0 , and its dc
value.
(c) Provide the constants A, B and C in the output: y(t) = A + B cos(πt + C).
Solution
(a) The fundamental frequency is Ω0 = π (T0 = 2), so only the Fourier coefficients corre-
sponding to k = 0 and k = ±1, i.e., frequencies Ω = 0, ±Ω0 , are needed as the filter deletes
any other frequency component.
es − 2 + e−s
x1 (t) = r(t + 1) − 2r(t) + r(t − 1) → X1 (s) =
s2
1 −1 1 − cos(πk)
Xk = (2 cos(πk) − 2) = k 6= 0
2 k2 π2 k2 π2
2
X0 = 0.5 (by observation), X1 = X−1 = 2
π
(b) Yes, y(t) is periodic of period 2π and dc value 2X0 = 1.
(c) The steady state output is
4.8 Given the Fourier series representation for a periodic signal x(t), we can compute derivatives
of it just like with any other signal.
(a) Consider the periodic train of pulses shown in Fig. 4.5, compute its derivative y(t) =
dx(t)/dt and carefully plot it. Find the Fourier series of y(t).
(b) Use the Fourier series representation of x(t) and find its derivative to obtain the Fourier
series of y(t). How doe it compare to the Fourier series obtained above.
x(t)
2
··· ···
T0 = 1
P
Answer: y(t) = k 4j sin(πk/2)ej2πkt .
Solution
(a) The derivative of a period of x(t) in −0.5 ≤ t ≤ 0.5 gives a period
dx1 (t)
y1 (t) = = 2δ(t + 0.25) − 2δ(t − 0.25) − 0.5 ≤ t ≤ 0.5
dt
of the periodic signal y(t). The Laplace transform of y1 (t) is Y1 (s) = 2e0.25s − 2e−0.25s so that
the Fourier series of y(t) of period T0 = 1, and Ω0 = 2π (just like x(t)) is
1
Yk = Y1 (s)|s=j2πk = 4j sin(πk/2)
T0
1 sin(πk/2)
Xk = X1 (s) |s=j2kπ = 2
T0 πk
so that
X
x(t) = Xk ej2πkt
k
where
sin(πk/2)
Yk = j2πkXk = 4jkπ = 4j sin(πk/2)
πk
which coincides with the result in part (a).
4.9 Consider the integral of the Fourier series of the pulse signal p(t) = x(t) − 1 of period T0 = 1,
where x(t) is given in Fig. 4.5.
(a) Given that an integral of p(t) is the area under the curve, find and plot the function
Z t
s(t) = p(t)dt t≤1
−∞
using the Fourier series of p(t). What is the integral equal to?
(d) You can also compute the integral from the plot of p(t)
Z T0 /2
p(t)dt
−T0 /2
What is it? Does it coincide with the result obtained using the Fourier series? Explain
P
Answer: p(t) = k,odd
4ej2πkt /(kπ).
Solution
(a) The pulse p(t) has an average of zero, so its integral is zero in each period, for kT0 ≤ t ≤
(k + 1)T0 , k an integer, and T0 = 1. If s1 (t) is the integral of p1 (t) for 0 ≤ t ≤ 1 we have that
Z t Z t Z t
s1 (t) = p1 (τ )dτ = [x1 (τ ) − 1]dτ = x1 (τ )dτ − t 0≤t≤1
0 0 0
Notice that x1 (t) = 2 for 0 ≤ t ≤ 0.25 and for 0.75 ≤ t ≤ 1, and zero for 0.25 ≤ t ≤ 0.75, so for
the given times,
Z 0
s(0) = 2dτ − 0 = 0
0
Z 0.25
s(0.25) = 2dτ − 0.25 = 0.25
0
Z 0.25
s(0.5) = 2dτ − 0.5 = 0
0
Z 0.25
s(0.75) = 2dτ − 0.75 = −0.25
0
Z 0.25 Z 1
s(1) = 2dτ + 2dτ − 1 = 0
0 0.75
p1 (t) = x1 (t) − [u(t) − u(t − 1)] = [2u(t) − 2u(t − 0.25) + 2u(t − 0.75) − 2u(t − 1)] − [u(t) − u(t − 1)]
= u(t) − 2u(t − 0.25) + 2u(t − 0.75) − u(t − 1)
so that the Fourier series coefficients of p(t) are given by (T0 = 1, Ω0 = 2π)
1 −0.5s 0.5s
Pk = e (e − 2e0.25s + 2e−0.25s − e−0.5s )|s=j2πk
T0 s
1 −jπk
= e 2(sin(πk) − 2 sin(πk/2))
πk
4
= (−1)(k+1) sin(πk/2)
πk
where we have replaced sin(πk) = 0. The value of Pk = 0 when k is even since sin(πk/2) = 0
for k even, and for odd k since sin(πk/2) = (−1)(k+1) then Pk = 4/(πk). Since the average of
p(t) is zero, we have that
X 4 j2πkt
p(t) = e
kπ
k, odd
Since
0.25 −0.5s 0.5s 0.25
S(s) = e (e − 2e0.25s + 2e−0.25s − e−0.5s ) = P (s)
s2 s
then (T1 = 0.5 and Ω1 = 4π):
1 0.25
Sk = S(s)|s=jΩ1 k = 2 P (s)|s=j4πk
T1 s
(c)(d) The integral over a period [−T0 /2, T0 /2] = [−0.5, 0.5] is
Z 0.5 X 4 Z 0.5
p(τ )dτ = ej2πkτ dτ
−0.5 kπ −0.5
k, odd
X 4 ej2πkτ
= |0.5
kπ j2πk −0.5
k, odd
X 4 sin(πk)
= =0
kπ πk
k, odd
since sin(πk) = 0 for k odd. Since p(t) has zero mean, the integral over a period is zero, coin-
ciding with the above.
4.10 Addition of Periodic Signals — Consider a saw-tooth signal x(t) with period T0 = 2 and
period (
t 0≤t<1
x1 (t) =
0 otherwise
(a) Find the Fourier coefficients Xk using the Laplace transform. Consider the cases when k
is odd and even (k 6= 0). You need to compute X0 directly from the signal.
(b) Let y(t) = x(−t), find the Fourier coefficients Yk .
(c) The sum z(t) = x(t) + y(t) is a triangular function. Find the Fourier coefficients Zk and
compare them to Xk + Yk .
(d) Use MATLAB to plot x(t), y(t), z(t) and their corresponding magnitude line spectra.
Solution
(a) The period x1 (t) = r(t) − r(t − 1) − u(t − 1), and T0 = 2 so that Ω0 = π. The Fourier
coefficients are
1 1 − e−s e−s 1 − e−s (1 + s)
Xk = − =
2 s2 s s=jπk 2s2 s=jπk
(
1 1 j/(2πk) k 6= 0, even
= 2 2
(−1)k (1 + jπk) − 1 = 2 2
2 π k −(2 + jπk)/(2π k ) k odd
The average value which cannot be obtained from above is
Z
1 1
X0 = tdt = 0.25
2 0
(b) The Fourier series of y(t) is
X
y(t) = x(−t) = X−k ejkΩ0 t
k
so that
−j/(2πk)
k 6= 0, even
Yk = X−k = 0.25 k=0
−(2 − jπk)/(2π 2 k 2 ) k odd
(c) Adding the above signals z(t) = x(t)+y(t) we get triangular pulses. The Fourier coefficients
of z(t) are
0 k 6= 0, even
Zk = Xk + Yk = 0.5 k=0
2 2
−2/(π k ) k odd
Notice that the Fourier coefficients of x(t) and y(t) which are neither even nor odd are complex,
while the Fourier coefficients of z(t) are real given that it is an even signal.
A period of z(t) is
z1 (t) = r(t) − 2r(t − 1) + r(t − 2)
and by inspection Z0 = 0.5. which coincide with the ones obtained before. The following script
gives an approximation of z(t) using 41 harmonics, and shown in Fig. 4.6.
% Pr. 4_10
clear all; clf
% magnitude line spectra
X1(1)=0.25;
for k=1:2:40,
X1(k+1)=abs(-(2+j*pi*k)/(2*piˆ2*kˆ2));
end
for k=2:2:40,
X1(k+1)=abs(1/(2*pi*k));
end
kk=-40:40;
X=[fliplr(X1) X1(2:41)];
figure(1)
subplot(311)
stem(kk,X); grid;axis([-40 40 0 1]);ylabel(’|X_k|’)
Y1(1)=0.25;
for k=1:2:40,
Y1(k+1)=abs((2+j*pi*k)/(2*piˆ2*kˆ2));
end
for k=2:2:40,
Y1(k+1)=abs(-1/(2*pi*k));
end
Y=[fliplr(Y1) Y1(2:41)];
subplot(312)
stem(kk,Y); grid; axis([-40 40 0 1]);ylabel(’|Y_k|’)
Z=[fliplr(Y1+X1) Y1(2:41)+X1(2:41)];
subplot(313)
stem(kk,Z); grid; axis([-40 40 0 1])
xlabel(’k’);ylabel(’|Z_k|’)
% approximate of z(t)
t=0:0.001:10;
z=0.5*ones(1,length(t));
for k=1:2:10,
z=z-4*cos(k*pi*t)/(pi*k)ˆ2;
end
figure(2)
plot(t,z);grid
xlabel(’t (sec)’); ylabel(’z_1(t)’)
0.4 0.9
|Xk|
0.2
0.8
0
−40 −30 −20 −10 0 10 20 30 40 0.7
0.6
0.6
0.4
z1(t)
|Yk|
0.5
0.2
0.4
0
−40 −30 −20 −10 0 10 20 30 40
0.3
0.6
0.2
0.4
|Zk|
0.2 0.1
0 0
−40 −30 −20 −10 0 10 20 30 40 0 1 2 3 4 5 6 7 8 9 10
k t (sec)
Figure 4.6: Problem 10: Line spectra for x(t), y(t) and z(t), (right) approximate ẑ(t).
4.11 Fourier series coefficients via Laplace — The computation of the Fourier series coefficients is helped
by the relation between the formula for these coefficients and the Laplace transform of a period of the
periodic signal.
(a) A periodic signal x(t) of period T0 = 2 sec., has as a period the signal x1 (t) = u(t) − u(t − 1). Use
the Laplace transform of x1 (t) to obtain the Fourier coefficients of x(t).
(b) Use MATLAB to approximate x(t) with 40 harmonics, plot the approximate signal x̂(t) and its mag-
nitude line spectrum.
Solution
(a) x(t) is a train of rectangular pulses.The Laplace transform of the period x1 (t) is
(b) The following script is used to find and plot the magnitude line spectrum of x(t), and to approximate
the signal with 40 harmonics.
1.2
1
0.5
0.8
0.4
0.6
xa(t)
|Xk|
0.3
0.4
0.2
0.2
0.1
0
0 −0.2
−40 −30 −20 −10 0 10 20 30 40 0 1 2 3 4 5 6 7 8 9 10
t (sec)
Figure 4.7: Problem 11: Magnitude line spectrum and approximate of train of pulses.
% Pr. 4_11
clear all; clf
% magnitude line spectra
X(1)=0.5;
for k=1:40,
X(k+1)=abs(0.5*sin(pi*k/2)/(pi*k/2));
end
kk=-40:40;
X=[fliplr(X) X(2:41)];
figure(1)
stem(kk,X); grid;axis([-40 40 0 0.6]);ylabel(’|X_k|’)
% approximate of z(t)
t=0:0.001:10;
x=0.5*ones(1,length(t));
for k=1:40,
x=x+(sin(pi*k/2)/(pi*k/2))*cos(pi*k*t-k*pi/2);
end
figure(2)
plot(t,x);grid
xlabel(’t (sec)’); ylabel(’x_a(t)’)
4.12 Time-support and frequency content — The support of a period of a periodic signal relates inversely to
the support of the line spectrum. Consider two periodic signals: x(t) of period T0 = 2 and y(t) of period
T1 = 1 and with periods
(a) Find the Fourier series coefficients for x(t) and y(t).
(b) Use MATLAB to plot the magnitude line spectra of the two signals from 0 to 20π rad/sec. Plot them
on the same figure so you can determine which has a broader support. Indicate which signal is
smoother and explain how it relates to its line spectrum.
Answers: y1 (t) = x1 (2t), line spectrum of x(t) is compressed when compared to that of y(t).
Solution
The periods of the two signals are:
% Pr. 4_12
clear all; clf
% magnitude line spectra
X(1)=0.5;Y(1)=X(1);
for k=1:40,
X(k+1)=abs(0.5 *sin(pi*k/2)/(pi*k*0.5));
end
Y=X;
kk=-40:40;
X=[fliplr(X) X(2:41)];
Y=[fliplr(Y) Y(2:41)];
figure(1)
subplot(211)
stem(pi*kk,X); grid;axis([-40*pi 40*pi 0 0.6]);ylabel(’|X_k|’)
subplot(212)
stem(2*pi*kk,Y); grid;axis([-40*2*pi 40*2*pi 0 0.6]);ylabel(’|Y_k|’)
0.5
0.4
|Xk|
0.3
0.2
0.1
0
−100 −50 0 50 100
0.5
0.4
|Yk|
0.3
0.2
0.1
0
−250 −200 −150 −100 −50 0 50 100 150 200 250
Figure 4.8: Problem 12: The line spectrum of Yk is an expansion of that of Xk (notice the frequencies
for each).
4.13 Windowing and music sounds — In the computer generation of musical sounds, pure tones
need to be windowed to make them more interesting. Windowing mimics the way a musician
would approach the generation of a certain sound. Increasing the richness of the harmonic fre-
quencies is the result of the windowing as we will see in this problem. Consider the generation
of a musical note with frequencies around fA = 880 Hz. Assume our “musician” while playing
this note uses three strokes corresponding to a window w1 (t) = r(t)−r(t−T1 )−r(t−T2 )+r(t−
T0 ), so that the resulting sound would be the multiplication, or windowing, of a pure sinusoid
cos(2πfA t) by a periodic signal w(t) with w1 (t) a period that repeats every T0 = 5T where T is
the period of the sinusoid. Let T1 = T0 /4, and T2 = 3T0 /4.
(a) Analytically determine the Fourier series of the window w(t) and plot its line spectrum
using MATLAB. Indicate how you would choose the number of harmonics needed to
obtain a good approximation to w(t).
(b) Use the modulation or the convolution properties of the Fourier series to obtain the coeffi-
cients of the product s(t) = cos(2πfA t)w(t). Use MATLAB to plot the line spectrum of this
periodic signal and again determine how many harmonic frequencies you would need to
obtain a good approximation to s(t).
(c) The line spectrum of the pure tone p(t) = cos(2πfA t) only displays one harmonic, the one
corresponding to the fA = 880 Hz frequency, how many more harmonics does s(t) have?
To listen to the richness in harmonics use the function sound to play the sinusoid p(t) and
s(t) (use Fs = 2 × 880 Hz to play both).
(d) Consider a combination of notes in a certain scale, for instance let
Use the same windowing w(t), and let s(t) = p(t)w(t). Use to plot p(t) and s(t) and to
compute and plot their corresponding line spectra. Use sound to play p(nTs ) and s(nTs )
using Fs = 1000.
Answers: Since the fundamental period of the sinusoid is T = 1/880 let T0 = 5T = 5/880 sec.,
T1 = T0 /4 and T2 = 3T0 /4.
Solution
(a) The fundamental period of the sinusoid is T = 1/880 so that T0 = 5T = 5/880 sec.,
T1 = T0 /4 and T2 = 3T0 /4. The Fourier series coefficients of ω(t) are
1 1
Wk = W1 (s)|s=jkΩ0 = 2 (1 − e−sT0 /4 − e−s3T0 /4 + e−sT0 )|s=jkΩ0
T0 s
e−sT0 /2 sT0 /2
= (e − esT0 /4 − e−sT0 /4 + e−sT0 /2 )|s=jkΩ0
s2
2e−jkπ 2π
= 2 2
(cos(πk) − cos(πk/2)) Ω0 = = 352π rad/sec
−Ω0 k T0
(b) The signal s(t) = w(t) cos(2πfA t) shifts the harmonics of w(t) to a central frequency 2πfA .
The number of harmonic frequencies is increased by windowing the sinusoid or sum of sinu-
soids.
%Pr. 4_13
clear all; clc
% parameters
Fs = 880; %frequency
T = 1/Fs; %sample period
T0 = 5*T; T1 = T0/4; T2 = T0*3/4;
% a)----------------------------------------------------------
% define one period of w signal -> w1
Ts=1/17600; % overall sample period
t=0:Ts:T0-Ts; Nn=length(t);
x1=[t(1:ceil(Nn/4)) t(ceil(Nn/4)+1)*ones(1,ceil(Nn/4))]; x2=fliplr(x1);
w1=[x1 x2];
N=length(w1); %number of samples in w1
% plot line spectrum
figure(1)
FSeries(w1,Ts,T0,N);
% harmonics H to approximate w1?
w0=2*pi/T0;
DC=0.001065; H=3;
for k=1:H,
a(k)=2*(cos(k*w0*T0/4)-cos(k*w0*T0/2))/T0/kˆ2/w0ˆ2;
end
a=[DC 2*a];
figure(2)
for N1=2:H+1,
x=a(1)*ones(1,length(t));
for k=2:N1,
x=x+a(k)*cos((k-1)*w0.*(t+T0/2));
end
plot(t,x);
axis([0 max(t) 1.1*min(x) 1.1*max(x)]);
hold on;
plot(t,w1,’r’)
ylabel(’Amplitude’);xlabel(’t (sec)’);
legend(’aprox w1(t)’,’w1(t)’);grid;
hold off
pause(0.1);
clear x;
end
clear H a x1
% b)----------------------------------------------------------
s = cos(2*pi*Fs*t).*w1; %modulated signal s
figure(3);
FSeries(s,Ts,T0,N); % line spectrum of s
% how many harmonics H are necessary to approximate s?
w0=2*pi/T0;
DC=0.001065;H=3;
for k=1:H,
a(k)=2*(cos(k*w0*T0/4)-cos(k*w0*T0/2))/T0/kˆ2/w0ˆ2;
end
a=[DC 2*a];
figure(4)
for N1=2:H+1,
x=a(1)*ones(1,length(t));
for k=2:N1,
x=x+a(k)*cos((k-1)*w0.*(t+T0/2));
x1=x.*cos(5*w0.*t);
end
plot(t,x1);
axis([0 max(t) 1.1*min(s) 1.1*max(s)]);
hold on;
plot(t,s,’r’)
ylabel(’Amplitude’); xlabel(’t (sec)’);
legend(’aprox s(t)’,’s(t)’); grid;
hold off
pause(0.1);
clear x;
end
% c)----------------------------------------------------------
p = cos(2*pi*Fs*t); %signal p
figure(5);
FSeries(p,Ts,T0,N); %line spectrum of p
%sound(p); pause(1); sound(w0*s);pause(2);
% d)----------------------------------------------------------
p = sin(2*pi*440*t)+sin(2*pi*550*t)+sin(2*pi*660*t); %signal p
s = p.*w1; %modulated signal s
figure(6);
FSeries(s,Ts,T0,N); %line spectrum of s
% sound(100*p); pause(1); sound(100*w0*s);
−3 −4 −3
x 10 Periodic signal x 10 x 10 Periodic signal
x(t)
0
1 12
−1
0.5 −2
0 10
Amplitude
<Xk|
|Xk|
|Xk|
0.6 0 3 0
4
0.4 2
−0.5
0.2 1
2
0 −1 0 −5
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 0 1 2 3 4 5 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
w (rad/sec) 4 w (rad/sec) 4 t (sec) −3 w (rad/sec) 4 w (rad/sec) 4
x 10 x 10 x 10 x 10 x 10
Figure 4.9: Problem 13. Right to left: window and line spectra, approximation of window, win-
dowed sinusoid and its line spectra.
−3 −3
x 10 Periodic signal x 10 Periodic signal
2
1.5
aprox s(t) 2
s(t) 1
0
1
x(t)
x(t)
0
−2
−1 −4
0.5
−2
Amplitude
<Xk|
|Xk|
|Xk|
0 0
−1
0.2 2 −0.5
0.1 1 −1
−1.5 −5
0 0 −1.5
0 1 2 3 4 5 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
t (sec) −3 w (rad/sec) 4 w (rad/sec) 4 w (rad/sec) 4 w (rad/sec) 4
x 10 x 10 x 10 x 10 x 10
Figure 4.10: Problem 13: from right to left: approximation of windowed signal, line spectra of
sinusoid, windowed sum of sinusoids and line spectra.
4.14 Computation of π — As you know, π is an irrational number that can only be approximated
by a number with a finite number of decimals. How to compute this value recursively is a
problem of theoretical interest. In this problem we show that the Fourier series can provide
that formulation.
and period T0 = 1. Plot the periodic signal and find its trigonometric Fourier series.
(b) Use the above Fourier series to find an infinite sum for π.
(c) If πN is an approximation of the infinite sum with N coefficients, and π is the value given
by , find the value of N so that πN is 95% of the value of π given by MATLAB.
P∞
Answer: π = 4 k=1 sin(πk/2)/k
Solution
(a) The signal x(t) is zero-mean. The Fourier series coefficients are obtained from (Ω0 =
2π/T0 = 2π)
2 0.25s sin(πk/2)
X1 (s) = (e − e−0.25s )|s=j2kπ =
s πk/2
so that
∞
X sin(πk/2) j2kπt
x(t) = e
πk/2
k=−∞,6=0
X∞
sin(πk/2)
= 2 cos(2πkt)
πk/2
k=1
Since x(0) = 1, letting t = 0 in both sides and then multiplying by π we get that
∞
X sin(πk/2)
π=4
k
k=1
If we let
N
X sin(πk/2)
πN = 4 = 0.95π
k
k=1
% Pr. 4_14
clear all; clc
N=100; x = 0; %DC value
for k=1:N
x = x + 4*sin(k*pi/2)/k;
xx(k) =x;
end
format long
display(’Error of approximation for’)
N
e = pi-x
ppi=pi*ones(1,length(xx)); n=1:N;
figure(1)
subplot(211); plot(n,xx); title(’Approximation of PI’); axis([1 N 0 4.5])
hold on; plot(n,ppi,’r’); grid; hold off; legend(’approximate’, ’pi’)
subplot(212); plot(n,xx-pi); title(’Error of approximation’); grid
% N for 95 percent of pi
sum = 0; %initial approximation
for i = 1:N,
sum = sum + 4*cos(i*pi/2-pi/2)/i;
if (sum >=(0.95*pi))&&(sum<=pi) break;
end
end
display([num2str(i) ’ coeficients needed of Fourier series for 95% approx. of PI ’])
Approximation of PI
2 approximate
pi
1
0
10 20 30 40 50 60 70 80 90 100
n
Error of approximation
1
0.5
−0.5
0 10 20 30 40 50 60 70 80 90 100
n
4.15 Square error approximation of periodic signals — To understand the Fourier series consider
a more general problem, where a periodic signal x(t), of period T0 , is approximated as a finite
sum of terms
N
X
x̂(t) = X̂k φk (t)
k=−N
where {φk (t)} are ortho-normal functions. To pose the problem as an optimization problem,
consider the square error Z
ε= |x(t) − x̂(t)|2 dt
T0
(a) Assume that x(t) as well as x̂(t) are real valued, and that x(t) is even so that the Fourier
series coefficients Xk are real. Show that the error can be expressed as
Z N
X Z N
X
2
ε= x (t)dt − 2 X̂k x(t)φk (t)dt + |X̂` |2 T0 .
T0 k=−N T0 `=−N
(b) Compute the derivative of ε with respect to X̂n and set it to zero to minimize the error.
Find X̂n .
(c) In the Fourier series the {φk (t)} are the complex exponentials and the {X̂n } coincide with
the Fourier series coefficients. To illustrate the above procedure consider the case of the
pulse signal x(t), of period T0 = 1 and a period
Use MATLAB to compute and plot the approximation x̂(t) and the error ε for increasing
values of N from 1 to 100.
(d) Concentrate your plot of x̂(t) around one of the discontinuities, and observe the Gibb’s
phenomenon. Does it disappear when N is very large. Plot x̂(t) around the discontinuity
for N = 1000.
R
Answer: dε/dX̂n = −2 T0
x(t)φn (t)dt + 2T0 X̂n = 0.
Solution
(a) Assuming x(t) and x̂(t) being real and even, so that X̂k are real we have
Z Z
ε = (x(t) − x̂(t))2 dt = x2 (t) − 2x(t)x̂(t) + x̂2 (t)dt
T T0
Z 0 Z X Z XX
= 2
x (t)dt − 2 x(t) X̂k φk (t)dt + X̂k X̂`∗ φk (t)φ∗` (t)
T0 T0 k T0 k `
Z XZ XX Z
= x2 (t)dt − 2 x(t)X̂k φk (t)dt + X̂k X̂` φk (t)φ∗` (t)
T0 k T0 k ` T0
To minimize the error we set its derivative w.r.t. X̂n to zero, i.e.,
Z Z
dε 1
= −2 x(t)φn (t)dt + 2T0 X̂n = 0 so that X̂n = x(t)φn (t)dt
dX̂n T0 T0 T0
% Pr. 4_15
clear all; clc
% c)d)--------------------------------------------------------
% Fourier Series of a train of pulses of period T0=1 and dc value 1
w0 = 2*pi; T0 = 1; DC=1; N=100;
for k=1:N,
a(k)=sin(k*pi/2)/(k*pi/2);
end
a=[DC 2*a];
figure(3)
plot(t,x); axis([0.175 0.275 1.75 2.25]); hold on
plot(t,x1,’r’); title(’Gibbs Phenomenon’);
ylabel(’x(t), x_a(t)’); xlabel(’t (sec)’);
pause(0.05)
%compute error
e=(x1-x).ˆ2;
error(N1-1)=0;
for k = 1:length(e)
error(N1-1) = error(N1-1) + e(k);
end
clear x
end
figure(2)
axis([1 N 0 1.1*max(error)]); hold on
plot(error,’linewidth’,2);
title(’e = int((x-x1)ˆ2)’); ylabel(’Amplitude of error’); xlabel(’t (samples)’)
hold off e = int((x−x1)2)
200
180
2
160
1.5 140
Amplitude of error
120
x(t), xa(t)
100
1
80
60
0.5
40
20
0
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 20 30 40 50 60 70 80 90 100
t (sec) t (samples)
Gibbs Phenomenon
2.25
2.2
2.15
2.1
2.05
x(t), xa(t)
1.95
1.9
1.85
1.8
1.75
0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27
t (sec)
Figure 4.12: Problem 15: Gibb’s phenomenon: (top–left) square pulse and approximate, (bottom)
detail of approximation around the discontinuity; (top-right) approximation error.
Solution
1
Chaparro-Akan — Signals and Systems using MATLAB 5.2
(c) Ω = s/j
2 −2 1 1
X(s) = = = + then
1 + (s/j)2 (s − 1)(s + 1) 1−s s+1
x(t) = et u(−t) + e−t u(t)
5.2 Find the Fourier transform of δ(t − τ ) and use it to find the Fourier transforms of
Answers: F[δ(t − 1) + δ(t + 1)] = 2 cos(Ω) and apply duality for the others.
(b) By duality
by duality
(a) Find the Fourier series of z(t) = d2 x(t)/dt2 using the Laplace transform. Then use the
derivative property to find the Fourier transform of x(t).
(b) To check your results figure out what the dc value of x(t) and z(t) should be, and whether
the Fourier coefficients should be real, purely imaginary or complex.
(c) Show that the following equation can be used to compute the Fourier transform of the
periodic signal z(t):
∞
2π X
Z(Ω) = Z1 (Ω) δ(Ω − kΩ0 )
T0
k=−∞
where Z1 (Ω) is the Fourier transform of z1 (t) = d2 x1 (t)/dt2 . Find its inverse z(t).
P
Answer: z(t) = m 2ej(2m+1)πt .
Solution
(a) A period of z(t) is z1 (t) = d2 x1 (t)/dt2 = δ(t) − 2δ(t − 1) + δ(t − 2), with period T0 = 2
(Ω0 = π), its Fourier series coefficients are
1 1
Zk = Z1 (s)|s=jπk = 2e−s (cosh(s) − 1)|s=jπk
2 2 (
k+1 k 2 k odd
= (−1) (1 − (−1) ) =
0 k even
then
∞
X ∞
X
z(t) = 2ejkπt = 2ej(2m+1)πt
k=−∞,odd m=−∞
Zk −2
Xk = = 2 2 , k odd
(jkΩ0 )2 k Ω0
and
∞
X ∞
X
2(−1)k jkπt 2(−1)2m+1 j(2m+1)πt
x(t) = e = e
k 2 Ω20 m=−∞
(2m + 1)2 Ω20
k=−∞,odd
(b) The average X0 = 0.5, and x(t) − 0.5 is even so Xk , k 6= 0 are real and
z(t) is even and zero mean, so Zk are real and Z0 = 0
5.4 A sinc signal x(t) = sin(0.5t)/(πt) is passed through an ideal low-pass filter with a frequency
response H(jΩ) = u(Ω + 0.5) − u(Ω − 0.5)
(a) Find the Fourier transform X(Ω) and carefully plot it.
(b) Find the output y(t) of the filter by calculating first Y (Ω). What is the convolution integral
of x(t) with x(t)? Explain.
Solution
|Ω|
X(Ω) = [u(Ω + π) − u(Ω − π)]
π
(a) Carefully plot X(Ω) as function of Ω.
(b) Determine the value of x(0)
Solution
X(Ω)
Ω
−π π
(b)
Z π Z π
1 |Ω| 2 Ω 1
x(0) = dΩ = dΩ =
2π −π π 2π 0 π 2
5.6 The Fourier series coefficients of a periodic signal x(t), with fundamental frequency Ω0 = π/4
∗
are X1 = X−1 = j, X5 = X5∗ = 2 and the rest are zero. Suppose x(t) is the input to a band-pass
filter with the following magnitude and phase responses
1 π ≤ Ω ≤ 1.5π
|H(jΩ)| = 1 −1.5π ≤ Ω ≤ −π ∠H(jΩ) = −Ω.
0 otherwise
Solution
(a) We have
X(Ω) = 2π [X1 δ(Ω − π/4) + X−1 δ(Ω + π/4) + X5 δ(Ω − 5π/4) + X−5 δ(Ω + 5π/4)]
= 2π [jδ(Ω − π/4) − jδ(Ω + π/4) + 2δ(Ω − 5π/4) + 2δ(Ω + 5π/4)]
Y (Ω) = 4πe−j5π/4 δ(Ω − 5π/4) + 4πej5π/4 δ(Ω + 5π/4)
dy(t)
= 0.5[x(t) − x(t − 2)]
dt
where x(t) is the input and y(t) the output. If x(t) = u(t) − 2u(t − 1) + u(t − 2)
Solution
(a) The FT of the equation characterizing the averager is
so that
sin(Ω)e−jΩ
Y (Ω) = X(Ω)
Ω
where the FT of x(t) is
(b) Calling Z(Ω) = X(Ω)/(jΩ) we have that y(t) = 0.5z(t) − 0.5z(t − 2) and
2
X(Ω) sin(Ω/2)
Z(Ω) = = e−jΩ
jΩ Ω/2
so that z(t) is the integral of x(t) = u(t) − 2u(t − 1) + u(t − 2), or z(t) = r(t) − 2r(t − 1) + r(t − 2),
replacing it in the above expression for y(t) we get
5.8 Fourier transform from Laplace transform of infinite support signals — For signals with
infinite support, their Fourier transforms cannot be derived from the Laplace transform unless
they are absolutely integrable or the region of convergence of the Laplace transform contains
the jΩ axis. Consider the signal x(t) = 2e−2|t|
Solution
(a) This is a decaying exponential for both negative and positive times.
(b) The signal is absolutely integrable,
Z ∞ Z ∞
e−2t ∞
|x(t)|dt = 2 2e−2t dt = 4 | =2
−∞ 0 −2 t=0
(c) Using the integral definition of the Fourier transform
Z 0 Z ∞
X(Ω) = 2 e2t e−jΩt dt + 2 e−2t e−jΩt dt
−∞ 0
t(2−jΩ) −t(2+jΩ)
2e 2e
= |0 − |∞
2 − jΩ −∞ 2 + jΩ 0
2 2
= +
2 − jΩ 2 + jΩ
8
=
4 + Ω2
(d) The Laplace transform of X(s) is
2 2 8
X(s) = + =
−s + 2 s + 2 4 − s2
with ROC: −2 ≤ σ ≤ 2, which includes the jΩ-axis. so that
8
X(Ω) = X(s)|s=jΩ =
4 + Ω2
%% Pr. 5_8
syms t w x
x=2*exp(-2*abs(t)); X=fourier(x)
figure(1)
subplot(211)
ezplot(x,[-2,2]); grid; axis([-2 2 0 2.2])
subplot(212)
ezplot(abs(X),[-30,30]);grid;axis([-20 20 0 2.2]);
xlabel(’\Omega (rad/sec)’); ylabel(’|X(\Omega)|’)
2/exp(2 abs(t))
2
1.5
0.5
0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
t
8/abs(w2 + 4)
2
1.5
|X(Ω)|
0.5
0
−20 −15 −10 −5 0 5 10 15 20
Ω (rad/sec)
5.9 Smoothness and frequency — Let the signal x(t) = r(t+1)−2r(t)+r(t−1) and y(t) = dx(t)/dt.
Solution
(a) The triangular pulse x(t) is even
1 s es es
X(s) = 2
[e − 2 + e−s ] = 2 [1 − 2e−s + e−2s ] = 2 (1 − e−s )2
s s s
which cancels out the poles at s = 0 so that the region of convergence is the whole s-plane. By
letting s = jΩ we have that
2 sin2 (Ω/2)
X(Ω) = [cos(Ω) − 1] =
−Ω2 (Ω/2)2
sin2 (Ω/2)
Y (Ω) = jΩX(Ω) = jΩ
(Ω/2)2
|Y (Ω)| = |X(Ω)||Ω|
which is the product of the sinc square and the absolute value of Ω. Thus x(t) is smoother than
its derivative y(t).
%% Pr. 5_9
syms x t w y X Y
x=(t+1)*(heaviside(t+1)-heaviside(t))+(1-t)*(heaviside(t)-heaviside(t-1));
y=heaviside(t+1)-2*heaviside(t)+heaviside(t-1);
figure(1)
subplot(211)
ezplot(x,[-1.5,1.5]);grid
subplot(212)
ezplot(y,[-1.5,1.5]);grid
X=int(x*exp(-j*w*t),t,-Inf,-Inf);
Y=int(y*exp(-j*w*t),t,-Inf,-Inf);
figure(2)
subplot(211)
XX=X*conj(X);
ezplot(10*log10(XX),[-100,100]);grid;axis([-100 100 -80 1]);ylabel(’X(\Omega)’)
subplot(212)
YY=Y*conj(Y);
ezplot(10*log10(YY),[-100,100]);grid;axis([-100 100 -80 1]);ylabel(’Y(\Omega)’)
0.6 −20
X(Ω)
0.4 −40
0.2
−60
0
−80
−1.5 −1 −0.5 0 0.5 1 1.5 −100 −80 −60 −40 −20 0 20 40 60 80 100
t w
heaviside(t+1)−2 heaviside(t)+heaviside(t−1) 10 log(16/w sin(1/2 w)2 conj(1/w sin(1/2 w)2))/log(10)
1 0
0.5 −20
Y(Ω)
0 −40
−0.5 −60
−1 −80
−100 −80 −60 −40 −20 0 20 40 60 80 100
−1.5 −1 −0.5 0 0.5 1 1.5 w
t
Figure 5.3: Problem 9: signals and their magnitude spectra (in dB).
6.1 Consider a series RC circuit with input a voltage source vi (t) and output the voltage across the
capacitor vo (t).
(a) Draw a negative feedback system for the circuit using an integrator, a constant multiplier
and an adder.
(b) Let the input be a battery, i.e., vi (t) = Au(t), find the steady state error e(t) = vi (t) − vo (t).
Answer: Vo (s)/Vi (s) = G(s)/(1 + G(s)) with G(s) = 1/(RCs).
Solution
(a) The transfer function of the RC circuit is
Vo (s) 1/RC
=
Vi (s) s + 1/(RC)
An equivalent negative feedback with a feed–forward transfer function G(s), and feedback
transfer function H(s) = 1 gives
Vo (s) G(s)
=
Vi (s) 1 + G(s)
Comparing the two we get that G(s) = 1/(RCs). Considering an integrator has a transfer
function 1/s then we have the equivalent representations for the RC circuit shown in Fig. 6.1.
(b) If vi (t) = Au(t), Vi (s) = A/s, then
A/RC B D
Vo (s) = = +
s(s + 1/(RC)) s s + 1/RC
1
Chaparro-Akan — Signals and Systems using MATLAB 6.2
+ e(t) 1 !
vi (t) (.)dt vo (t)
RC
−
Since 1/RC > 0, i.e., the pole of Vo (s) is in the left-hand s-plane, then in the steady state
vo (t) = B, where
A/RC
B = Vo (s)s|s=0 = =A
1/(RC)
so that in the steady–state the error signal e(t) = vi (t) − vo (t) will be zero.
6.2 Suppose the feed–forward transfer function of a negative feedback system is G(s) = N (s)/D(s),
and the feedback transfer function is unity. Let X(s) be the Laplace transform of the input x(t)
of the feedback system.
(a) Given that the Laplace transform of the error is E(s) = X(s)[1 − F (s)] where F (s) =
G(s)/(1 + G(s)) is the overall transfer function of the feedback system, find an expres-
sion for the error in terms of X(s), N (s) and D(s). Use this equation to determine the
conditions under which the steady state error is zero for x(t) = u(t).
(b) If the input is x(t) = u(t), N (s) = 1, and D(s) = (s + 1)(s + 2) find an expression for E(s)
and from it determine the initial value e(0) and the final value limt→∞ e(t) of the error.
Answers: X(s) = 1/s, zero steady–state error: roots of N (s) + D(s) in left–hand s-plane, D(s)
of the form sD1 (s).
Solution
(a) Replacing F (s) we have
G(s) 1
E(s) = X(s) 1 − = X(s)
1 + G(s) 1 + N (s)
D(s)
X(s)D(s)
=
D(s) + N (s)
If X(s) = 1/s then for the error to go to zero in the steady state, the roots of D(s) + N (s) should
be in the left–hand s-plane and D(s) must cancel the pole of E(s) at zero contributed by X(s),
i.e., D(s) must be of the form D(s) = sD1 (s).
(b) For G(s) = 1/(s + 1)(s + 2) and X(s) = 1/s the Laplace transform of the error is
D(s) (s + 1)(s + 2)
E(s) = =
s(1 + D(s)) s(s2 + 3s + 3)
(1 + 1/s)(1 + 2/s)
e(0) = lim sE(s) = lim =1
s→∞ s→∞ 1 + 3/s + 3/s2
p
The poles of E(s) are s = 0 and s1,2 = −3/2 ± j 3/4 since s2 + 3s + 3 = (s + 3/2)2 + 3/4. A
partial fraction expansion for E(s) is
A Bs + C
E(s) = + 2
s s + 3s + 3
where
D(0) 2
A = sE(s)|s=0 = =
1 + D(0) 3
in this case the steady state error is 2/3. The error is not zero in the steady–state because D(s)
does not cancel the pole due to X(s).
6.3 Let H(s) = Y (s)/X(s) be the transfer function of the feedback system in Fig. 6.2. The impulse
response of the plant (with transfer function Hp (s)) is hp (t) = sin(t)u(t).
(a) If we want the impulse response of the feedback system (with input x(t) and output y(t))
to be h(t) = hp (t)e−t u(t), what should be the transfer function Hc (s) of the controller?
(b) Find the unit-step response s(t) of the feedback system.
x(t) + y(t)
Hc (s) Hp (s)
−
Answers: Hc (s) = (s2 + 1)/(s + 1)2 ; s(t) = [0.5 + 0.707 cos(t + 3π/4)]u(t)
Solution
(a) h(t) = hp (t)e−t = e−t sin(t)u(t) so
1
H(s) =
(s + 1)2 + 1
Hc (s) 1 s2 + 1
= ⇒ Hc (s) =
s2 + 1 + Hc (s) 2
(s + 1) + 1 (s + 1)2
H(s) 1
S(s) = =
s s((s + 1)2 + 1)
A B B∗
= + +
s s − (−1 + j) s − (−1 − j)
where
A = sS(s)|s=0 = 0.5
1 1
B = |s=−1+j = √ ej3π/4
s(s + 1 + j) 2 2
thus
s(t) = [0.5 + 0.707 cos(t + 3π/4)]u(t).
6.4 Consider the cascade connection of two continuous-time systems shown in Fig. 6.3 where
dw(t) dy(t)
+ w(t) = x(t), + 2y(t) = w(t)
dt dt
w(t)
x(t) Syst. A Syst. B y(t)
(a) Determine the input/output differential equation for the overall cascade connection.
(b) Suppose that w(0) = 1, y(0) = 0 and x(t) = 0 for t ≥ 0,
i. Compute w(t) for t ≥ 0.
ii. Compute then y(t) for t ≥ 0.
Answers: (a) d2 y(t)/dt2 + 3dy(t)/dt + 2y(t) = x(t); (b) y(t) = (e−t − e−2t )u(t)
Solution
(a) Laplace transforms
X(s) Y (s) 1 1
(s + 2)Y (s) = ⇒ = = 2
s+1 X(s) (s + 1)(s + 2) s + 3s + 2
so that
d2 y(t) dy(t)
+3 + 2y(t) = x(t)
dt2 dt
(b) i. For the first system, the Laplace transform of its differential equation is
1
sW (s) − w(0) + W (s) = 0 ⇒ W (s) = ⇒ w(t) = e−t u(t)
s+1
ii. For the second system with input w(t) = e−t u(t)
1 1
sY (s) − y(0) + 2Y (s) = ⇒ Y (s) =
s+1 (s + 1)(s + 2)
so that
1 1
Y (s) = − ⇒ y(t) = (e−t − e−2t )u(t)
s+1 s+2
Y (s) 1
H(s) = = 2
X(s) s + s/Q + 1
where Y (s) and X(s) are the Laplace transforms of output y(t) and the input x(t) of the system.
Q is called the quality factor.
(a) If the feedback gain is unity, determine the feedforward transfer function G(s) .
(b) Find the poles of H(s), and plot them expressing the poles as p1,2 = re±jφ , where φ is
an angle measured with respect to the negative axis, and r is a positive radius, give the
values of r and φ. Show that Q = 1/(2 cos(φ)).
√
(c) Consider the cases when Q = 0.5, Q = 2/2 and Q → ∞, and determine the correspond-
ing poles. Find the impulse response h(t) for these cases. What happens as Q increases?
Explain.
Solution:
(a) Let G(s) = 1/(s(s + 1/Q)) be the Laplace transform in the feedforward loop and with unit
gain in the feedback then
1 p
p1,2 = − ± j 1 − 1/(4Q2 )
2Q
which have a magnitude of one and the angle φ measured with respect to the negative
real–axis is p
1 − 1/(4Q2 ) p
−1
φ = tan = tan−1 ( 4Q2 − 1)
1/2Q
The real part
1 1
= 1 cos(φ) → Q =
2Q 2 cos(φ))
(c) When Q = 0.5, 0.707, ∞ the angle φ = 0, π/4, π/2. For
1
H(s) = when
(s + 1/(2Q))2 + (1 − 1/(4Q2 ))
1
Q = 0.5 → H(s) =
(s + 1)2
1
Q = 0.707 → H(s) = √
(s + 1/ 2)2 + 1/2
1
Q → ∞ → H(s) = 2
s +1
6.6 Given the two realizations in Fig. 6 obtain the corresponding transfer functions
Y1 (s) Y2 (s)
H1 (s) = , and H2 (s) = .
X1 (s) X2 (s)
Realization 1 +
y1 (t)
−2
� v1 (t) � v2 (t)
+
x1 (t)
2
Realization 2
x2 (t)
−2
w2 (t) w1 (t)
� �
+ +
y2 (t)
2
Replacing V2 (s) = V1 (s)/s from equation (iii) in equation (ii) we have that
so that V1 (s) = X1 (s)/(s − 1 − 2/s) and and using (iii) 2V2 (s) = 2V1 (s)/s = X1 (s)(2/s)/(s − 1 −
2/s) so that equation (i) becomes
1 − 2/s
Y1 (s) = X1 (s)
s − 1 − 2/s
so that
Y1 (s) s−2
H1 (s) = = 2
X1 (s) s −s−2
y2 (t) = w1 (t)
ẇ1 (t) = x2 (t) + y2 (t) + w2 (t)
ẇ2 (t) = −2x2 (t) + 2y2 (t)
Y2 (s) s−2
H2 (s) = = 2
X2 (s) s −s−2
Indicating that the two realizations are different realizations of the same system.
To find a transformation that diagonalizes Ao use MATLAB function eigs which calculates the
eigenvalues and eigenvectors corresponding to the matrix and allows us to express
λ1 0
V−1 Ao V =
0 λ2
where V is a matrix created with the eigenvectors and {λi }, i = 1, 2, are the eigenvalues.
(a) Find the characteristic equation
det(sI − Ao )
corresponding to the state variable representation, and show it is the same as the denom-
inator of the transfer function (use function ss2tf to obtain the transfer function from the
state variable representation).
(b) Use the matrix V as an invertible transform to obtain a new set of state variables with a
diagonal matrix A and vectors b and cT . Obtain these matrix and vectors.
(c) Suppose that you find the controller form by letting Ac = ATo , bc = cTo and cTc = bT0 , repeat
the above diagonalization and comment on your results.
Answers: H(s) = (s − 1)/((s − 1)(s + 2)) = 1/(s + 2)
Solution
(a)(b) To find the transfer function let the input v(t) = δ(t), y(t) = h(t) and initial conditions
equal zero, then in the Laplace domain
1 s 1 1 s−1
H(s) = cTo (sI − Ao )−1 b = 2 1 0 = 2
s +s−2 2 s+1 −1 s +s−2
% Pr. 6.7
% Observer
Ao=[-1 1; 2 0]
bo=[1 ;-1]
co=[1 0]’
% Transfer function
[num,den]=ss2tf(Ao,bo,co’,0)
roots(den)
% Diagonalization
[V,D]=eigs(Ao)
A=(inv(V))*Ao*V
b=(inv(V))*bo
c=V’*co
[num,den]= ss2tf(A,b,c’,0)
roots(den)
Ao = -1 1
2 0
bo = 1
-1
co = 1
0
num = 0 1.0000 -1.0000
den = 1 1 -2
ans = -2
1
V =-0.7071 -0.4472
0.7071 -0.8944
D = -2 0
0 1
A = -2.0000 0.0000
0 1.0000
b = -1.4142
0
c = -0.7071
-0.4472
num = 0 1 -1
den = 1 1 -2
ans =-2
1
Showing that the two realizations correspond to the same transfer function.
(c)
%Controller
Ac=Ao’
bc=co’
cc=bo’
[V,D]=eigs(Ao)
A=(inv(V))*Ao*V
b=(inv(V))*bo
c=V’*co
[num,den]= ss2tf(A,b,c’,0)
roots(den)
Ac = -1 2
1 0
bc = 1 0
cc =1 -1
V =-0.7071 -0.4472
0.7071 -0.8944
D =-2 0
0 1
A =-2.0000 0.0000
0 1.0000
b =-1.4142
0
c =-0.7071
-0.4472
num = 0 1 -1
den =1 1 -2
ans =-2
1
Again showing that the two realizations correspond to the same transfer function and that the
observer and the controller forms represent the same system.
s + z1
H(s) = K
s + p1
(a) If we want unity dc gain, i.e., |H(j0)| = 1, what should be the value of K.
(b) What types of filters (lowpass, bandpass, band-eliminating, all-pass and high-pass) can
you obtain with this first-order system?
(c) Let H1 (s) = 1/(s + 1), corresponding to a low-pass filter, determine the dc gain of this
filter and use it to find the half-power frequency of the filter.
(d) Determine the magnitude response of a filter with transfer function
H2 (s) = (s − 1)/(s + 1). Plot the poles and zeros of this filter. What type of filter is it?
(e) What type of filter has the transfer function H3 (s) = s/(s + 1) ?
Answers: Low-, high- and all-pass filters can be obtained with H(s); H2 (s) is all-pass filter.
Solution
1
Chaparro-Akan — Signals and Systems using MATLAB 7.2
(c) For the low-pass filter, |H1 (j0)| = 1, the half-power frequency Ωhp is such that
so Ωhp = 1.
(d) The filter with H2 (s) = (s − 1)/(s + 1) is an all-pass filter. Its pole is s = −1 and zero s = 1.
Given this symmetry the magnitude response is unity for all values of Ω.
(e) The filter with H3 (s) = s/(s + 1) is a high-pass filter, |H(j0)| = 0 and when Ω → ∞ the
magnitude response grows and becomes 1 at ∞.
7.2 The loss at a frequency Ω = 2000 (rad/sec) is α(2000) = 19.4 dBs for a fifth-order lowpass
Butterworth filter with unity dc gain. If we let α(Ωp ) = αmax = 0.35 dBs, determine
• the half-power frequency Ωhp , and
• the passband frequency Ωp
of the filter.
Answers: Ωhp = 1280.9 rad/sec; Ωp = 999.82 (rad/sec)
Solution
Using the loss function
2N !
Ω Ω
α(Ω) = 10 log10 1+ ⇒ Ωhp =
Ωhp (100.1α(Ω) − 1)1/2N
2000
Ωhp = = 1280.9 rad/sec.
(101.94 − 1)0.1
which gives after replacing the values on the right term Ωp = 999.82 (rad/sec).
Solution
Since the dc loss is not zero, the normalized loss specifications are
αmax = α1 − α(0) = 0.5
αmin = α2 − α(0) = 30
with a dc loss of 20 dB. The following script is used to find the answers
% Pr 7_3
clear all; clf
alphamax=0.5 ;alphamin=30; Wp=1500; Ws=3500;
alpha0=20;
% Butterworth
K=10ˆ(alpha0/20)
D=(10ˆ(0.1*alphamin)-1)/(10ˆ(0.1*alphamax)-1);
E=Ws/Wp;
N=ceil(log10(D)/(2*log10(E)))
Whp=Wp/(10ˆ(0.1*alphamax)-1)ˆ(1/(2*N))
alpha_p=10*log10(1+(Wp/Whp)ˆ(2*N))
alpha_s=10*log10(1+(Ws/Whp)ˆ(2*N))
% Chebyshev
N=ceil(acosh(Dˆ(0.5))/acosh(E))
eps=sqrt(10ˆ(0.1*alphamax)-1)
Whp1=Wp*cosh(acosh(1/eps)/N)
alpha_p= 10*log10(1+(epsˆ2)*(cos(N*acos(1)))ˆ2)
alpha_s=10*log10(1+(epsˆ2)*cosh(N*acosh(Ws/Wp))ˆ2)
alpha1=alpha0+10*log10(1+(epsˆ2)*cosh(N*acosh(0))ˆ2);
Kc=10ˆ(alpha1/20)
% Butterworth
K = 10 % dc gain
N =6 % minimum order
Whp =1.7874e+03 % half-power frequency
alpha_p = 0.5000 % loss at Wp
alpha_s =35.0228 % loss at Ws
% Chebyshev
N = 4 % min order
eps = 0.3493 % ripple factor
Whp1 = 1.6397e+03 % half-power frq
alpha_p = 0.5000 % loss at Wp
alpha_s = 36.6472 % loss at Ws
Kc =10.5925 % dc gain
Notice the computation of the dc gain in the Chebyshev filter. In this case the dc loss depends
on the order of the filter and so it is not necessarily 0 dB, so to get the dc gain Kc we use
K2 2 2
α(0) = 10 log10 2 (0) = 20 log10 K − 10 log10 (1 + CN (0))
1 + 2 C N
7.4 Getting rid of 60 Hz hum with different filters — A desirable signal x(t) = cos(100πt) −
2 cos(50πt) is recorded as y(t) = x(t) + cos(120πt), i.e., as the desired signal but with a a 60
Hz hum. We would like to get rid of the hum and recover the desired signal. Use symbolic
MATLAB to plot x(t) and y(t).
Consider the following three different alternatives (use symbolic MATLAB to implement the
filtering and use any method to design the filters):
• Design a band-eliminating filter to get rid of the 60 Hz hum in the signal. Plot the output
of the band-eliminating filter.
• Design a high-pass filter to get the hum signal and then subtract it from y(t). Plot the
output of the high-pass filter.
• Design a band-pass filter to get rid of the hum. Plot the output of the band-pass filter.
Is any of these alternatives better than the others? Explain.
Solution
The following script is used to compute the different signals, filters, and to plot them.
% Pr 7_4
clear all; clf
syms t w
x=cos(100*pi*t)-2*cos(50*pi*t);
y=x+cos(120*pi*t);
figure(1)
subplot(211)
ezplot(x,[0:1]) % desired signal
subplot(212)
ezplot(y,[0:1]) % signal with hum
X=fourier(x);
% butterworth filter
N=5;
%wn=[2*pi*59.9 2*pi*60.1]; [b,a]=butter(N,wn,’stop’,’s’);C=1 % band-eliminating
% [b,a]=butter(2*N,110*pi,’high’,’s’);C=2 % high-pass
[b,a]=butter(N,[20*pi 90*pi],’s’);C=2 % band-pass
% frequency responses
W=0:1:150*pi;
Hm=abs(freqs(b,a,W));
figure(2)
subplot(211)
plot(W/(2*pi),Hm); axis([0 75 0 1.1]);xlabel(’f (Hz)’); ylabel(’|H|’)
M=2*N
% generation of frequency response from coefficients
n=M:-1:0;
U=(j*w).ˆn;
num=b*conj(U’); den=a*conj(U’);
H=num/den; % Butterworth LPF
% output of filter
Y1=X*H;
y1=real(ifourier(Y1,t));
subplot(212)
if C==1,
ezplot(y1,[0:1])
else
ezplot(x-y1,[0:1])
end
title(’denoised signal’)
The most natural of the three approaches is the one using the stopband filter, the other depend
on information on the desired signal that might or might not be available. The only needed
information for the stopband filter is the frequency of the hum.
cos(100 π t)−2 cos(50 π t)
−1
4
3
2
1
0
−1
−2
1 1 1
|H|
|H|
0 0 0
0 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70
f (Hz) f (Hz) f (Hz)
denoised signal denoised signal denoised signal
3
3 1
2 2
0
1 1
0 0 −1
−1 −1
−2
−2
−2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t t t
Figure 7.1: Problem 4: top: original and original with hum signals; bottom left to right: notch filter
and filtered signal, high-pass filter and denoised signal, an bandpass filter and denoised signal.
Sampling Theory
Solution
(a) To find X(Ω), we use duality or find the inverse Fourier transform of a pulse of amplitude
A and bandwidth Ω0 , that is
X(Ω) = A[u(Ω + Ω0 ) − u(Ω − Ω0 )]
so that
Z Ω0
1 A jΩt Ω0
x(t) = AejΩt dΩ = e −Ω0
2π −Ω0 2πjt
A
= sin(Ω0 t)
πt
which when compared with the given x(t) = sin(t)/t gives that A = π and Ω0 = 1 or
X(Ω) = π[u(Ω + 1) − u(Ω − 1)]
1
Chaparro-Akan — Signals and Systems using MATLAB 8.2
1 Ωmax
fs = ≥2
Ts 2π
which gives a sampling period
π
Ts ≤ = π sec/sample
Ωmax
1
Y (Ω) = (X(Ω) ∗ X(Ω))
2π
which would have a maximum frequency Ωmax = 2, giving a sampling frequency which is
double the one for x(t). The sampling period for y(t) should be
π
Ts ≤ .
2
(d) The signal x(t) = sin(t)/t is zero whenever t = ±kπ, for k = 1, 2, · · · so that choosing Ts = π
(the Nyquist sampling period) we obtain the desired signal xs (0) = 1 and x(nTs ) = 0.
Solution
(a) X(Ω) = ejΩ + e−jΩ = 2 cos(Ω) which has no maximum frequency, so x(t) is not bandlim-
ited.
(b) Y (Ω) = X(Ω)H(jΩ) = 2 cos(Ω)[u(Ω + 1) − u(Ω − 1)] with maximum frequency of 1 and
so y(t) is band-limited.
(c) It is not time limited, as
Z 1 Z 1
1 1
y(t) = 2 cos(Ω)ejΩt dΩ = [ejΩ(t+1) + ejΩ(t−1) ]dΩ
2π −1 2π −1
sin(t + 1) sin(t − 1)
= +
π(t + 1) π(t − 1)
8.3 Suppose you wish to sample an amplitude modulated signal x(t) = m(t) cos(Ωc t) where m(t)
is the message signal and Ωc = 2π104 (rad/sec) is the carrier frequency.
(a) If the message is an acoustic signal with frequencies in a band of [0, 22] kHz, what would
be maximum frequency present in x(t)?
(b) Determine the range of possible values of the sampling period Ts that would allow to
sample x(t) satisfying the Nyquist sampling rate condition.
(c) Given that x(t) is a bandpass signal, compare the above sampling period with the one that
can be used to sample bandpass signals.
Solution
(a) The Fourier transform of x(t) is
where M (Ω) is the Fourier transform of m(t). The maximum frequency present in x(t) is
2π
Ωs = ≥ 2 × 64π × 103 = 128π × 103
Ts
so that
1
Ts ≤ 10−3 sec/sample
64
(c) The bandwidth of the message is B = 2π×22×103 = 44π×103 rad/sec, using this frequency
x(t) can be sampled at 2B = 88π×103 rad/sec which is much smaller than the one found above.
8.4 The input/output relation of a non-linear system is y(t) = x2 (t), where x(t) is the input and
y(t) is the output.
(a) The signal x(t) is band-limited with a maximum frequency ΩM = 2000π (rad/sec), deter-
mine if y(t) is also band-limited and if so what is its maximum frequency Ωmax .
(b) Suppose that the signal y(t) is low-pass filtered. The magnitude of the low-pass filter
is unity and the cut-off frequency is Ωc = 5000π (rad/sec), determine the value of the
sampling period Ts according to the given information.
(c) Is there a different value for Ts that would satisfy the Nyquist sampling rate condition for
sampling both x(t) and y(t) and that is larger than the one obtained above? Explain.
Solution
(a) If the signal x(t) is band-limited y(t) = x2 (t) has a Fourier transform Y (Ω) = (1/2π)(X(Ω) ∗
X(Ω)) having a bandwidth double that of x(t), or
so y(t) is band-limited.
(b) Filtering with a low-pass filter of cut-off frequency 5000π would not change the maximum
frequency of Y (Ω) so that
π 1
Ts ≤ = = 0.25 × 10−3
4000π 4000
(c) No. For x(t),
π
Ts1 ≤ = 0.5 × 10−3
2000π
and if Ts = 0.25 × 10−3 then Ts1 ≤ 2Ts , so that we need to use Ts to sample both x(t) and y(t).
8.5 A sinusoid x(t) = cos(t) is a band-limited signal with maximum frequency Ωmax = 1
(a) Using Fourier transform properties determine the maximum frequency of x2 (t). What
sampling period Ts can be used to sample x2 (t) without aliasing. Verify your results using
trigonometric identities.
(b) Can you generalize the above results to xN (t), for integer N > 1? Is xN (t) band-limited
for all the given values of N ? Consider the case N = 3 and use trigonometric identities to
verify the results.
Answers: (a) Ts ≤ π/2.
Solution
(a) If F(cos(t)) = X(Ω), then F(x2 (t)) = (X ∗ X)(Ω)/(2π), i.e., the convolution of X(Ω)
with itself, having twice the bandwidth of X(Ω), so the maximum frequency of x(t) is
Ωmax = 2, and Ts ≤ π/2. Thus x(t) is band-limited with twice the maximum frequency of
cos(t).
Using cos2 (t) = 0.5(1 + cos(2t)) we see that its maximum frequency is 2 and it is band-
limited and Ts ≤ π/2.
(b) In general F(xN (t)) would be N convolutions of X(Ω) with itself, and the bandwidth
being N ×bandwith of cos(t), i.e., N , so the maximum frequency is Ωmax = N and Ts ≤
π/N .
If N = 3 , cos3 (t) = 0.25(3 cos(t) + cos(3t)) with Ωmax = 3, and Ts ≤ π/3.
8.6 Let x(t) = cos(πt)[u(t) − u(t − 2)] be the input of a zero-order hold sampler. The sampler
samples every Ts sec. starting at t = 0, and the impulse response of the zero-order hold is
h(t) = u(t) − u(t − Ts ). The quantizer and coder shown in Fig. 7.12 in the text is used to
generate a digital signal.
(a) Let Ts = 0.5 sec/sample, find the values of the sampled signal using the ideal sampler.
Determine the corresponding digital signal.
(b) Suppose now that we let Ts = 0.25 sec/sample. Carefully plot the sampled signal using
the ideal sampler and then plot the output y(t) of the zero-order hold system. Determine
the corresponding digital signal.
(c) For the above two cases find and plot the quantization error ε(nTs ) = x(nTs ) − x̂(nTs ),
where x̂(nTs ) is the output of the quantizer.
Answer: The binary code is {01 00 10 00 01 00 · · · }.
Solution
(a) The sample values for Ts = 0.5 are
nTs x(nTs )
0 cos(0) = 1
0.5 cos(π/2) = 0
1.0 cos(π) = −1
1.5 cos(3π/2) = 0
2.0 cos(2π) = 1
x(nTs ) = 0 nTs > 2.
The binary code for nTs = 0, 0.5, 1, 1.5, 2, 2.5 · · · is
01 00 10 00 01 00 · · ·
corresponding to the quantized signal
x̂(nTs ) = ∆δ(nTs ) + 0δ((n − 1)Ts ) − 2∆δ((n − 2)Ts ) + 0δ((n − 3)Ts ) + ∆δ((n − 4)Ts )
where ∆ = 2/4 = 0.5, so that the quantization error is
ε(nTs ) = x(nTs ) − x̂(nTs )
= 0.5δ(nTs ) + 0δ((n − 1)Ts ) + 0δ((n − 2)Ts ) + 0δ((n − 3)Ts ) + 0.5δ((n − 4)Ts )
01 01 00 11 10 11 00 01 01 00, · · ·
Solution
(a) For x(t) to be sampled without aliasing we need Ωs > 2Ωmax = 2.
(b) If y(t) = x2 (t) then Y (Ω) = [X ∗ X](Ω)/(2π) which gives a triangular spectrum in [−2 2]
frequency band, i.e., it is band-limited.
(c) The maximum frequency of z(t) is Ωmax = 2 (rad/sec) so Ωs ≥ 4 (rad/sec) and Ts ≤ π/2
sec.
2∆(t)∆(Ω) ≥ 1
measures the frequency support for which the Fourier representation is significant. The
energy of the signal is represented by Ex . Use MATLAB to compute ∆(t) and ∆(Ω) for
the signals x(t) and x1 (t) and verify that the uncertainty principle is satisfied. Comment
on the difference in results for the two signals; pay attention to the smoothness and the
frequency content of the signals.
Solution
(a) For the Gaussian signal we notice that the signal and its spectrum look very similar, and
that the spectrum is narrower that the signal. For x1 (t), we notice that it has discontinuities
and its spectrum has high frequencies because of them. See Figs. 8.1 and 8.2.
(b) Computing ∆(t) and ∆(Ω) we verify the uncertainty bound. It is 1 for the Gaussian signal
(the only signal that will give this value) and 1.2197 for the square pulse x1 (t). You need to
make some simple changes to the script for x1 (t), the script shown is for the Gaussian signal.
% Pr. 8.8
clear all; clf
syms x t X X1 w z
sigma=1
%x=heaviside(t+sigma)-heaviside(t-sigma);
x=exp(-tˆ2/(2*sigmaˆ2));
X=fourier(x,w);
XX=X*conj(X);
figure(1)
subplot(211)
ezplot(x,[-4*sigma,4*sigma]);grid;
subplot(212)
ezplot(XX,[-4*sigma,4*sigma]);grid;
Ex=int(xˆ2,t,-4*sigma,4*sigma)
figure(2)
subplot(211)
ezplot((xˆ2)*(tˆ2)/Ex,[-4*sigma,4*sigma]);grid;
subplot(212)
ezplot(XX*(wˆ2)/Ex,[-4*sigma,4*sigma]);grid;
z=int((xˆ2)*(tˆ2),t, -4*sigma, 4*sigma)/Ex
zz=subs(sqrt(z)) % Delta(t)
v=int((Xˆ2)*(wˆ2)/(2*pi),w, -4*sigma,4*sigma)/Ex
vv=subs(sqrt(v)) % Delta(W)
2*zz*vv
0.2 0.05
0 0
−4 −3 −2 −1 0 1 2 3 4 −4 −3 −2 −1 0 1 2 3 4
t t
6
5 1
4
3
0.5
2
1
0 0
−4 −3 −2 −1 0 1 2 3 4 −4 −3 −2 −1 0 1 2 3 4
w w
Figure 8.1: Problem 8: Gaussian signal x(t) and its spectrum (both are real) (left); integrands used
to calculate ∆(t) and ∆(Ω) (right)
heaviside(t + 1) − heaviside(t − 1) (t2 (heaviside(t − 1) − heaviside(t + 1))2)/2
0.8
1
0.8 0.6
0.6 0.4
0.4
0.2
0.2
0 0
−4 −3 −2 −1 0 1 2 3 4 −4 −3 −2 −1 0 1 2 3 4
t t
((cos(w) i − sin(w))/w − (cos(w) i + sin(w))/w) ((cos(conj(w)) i − sin(conj(w)))/conj(w) − (cos(conj(w)) i + sin(conj(w)))/conj (w2 ((cos(w) i − sin(w))/w − (cos(w) i + sin(w))/w) ((cos(conj(w)) i − sin(conj(w)))/conj(w) − (cos(conj(w)) i + sin(conj(w)))/conj
4 2
3 1.5
2 1
1 0.5
0 0
−4 −3 −2 −1 0 1 2 3 4 −4 −3 −2 −1 0 1 2 3 4
w w
Figure 8.2: Problem 8: Rectangular pulse x1 (t) and its spectrum (left); integrands used to calculate
∆(t) and ∆(Ω) (right
9.1 For a finite-support signal x[n] = r[n](u[n] − u[n − 11]) where r[n] is the discrete-time ramp
function,
Answers: εx = 385.
Solution
1
Chaparro-Akan — Signals and Systems using MATLAB 9.2
The energy
∞
X
εx = (xe [n] + xo [n])2
n=−∞
X∞
= x2e [n] + x2o [n] + 2xe [n]xo [n]
n=−∞
X∞ ∞
X
= x2e [n] + x2o [n]
n=−∞ n=−∞
since xe [n]xo [n] is odd its sum is zero. Also from plots of xe [n] and xo [n] we have that
X10 X10
n2
εx = 4 = n2 = 385
n=0
4 n=0
Solution
(a) Yes, the quantizer is time-invariant. Indeed, x(nTs ) and x(nTs − M Ts ) will give the same
values just shifted in time. In this case the M must be an integer as the sampled signal is dis-
crete in time.
(b) The samples are given by x(nTs ) = nTs = 0.1n. For ∆ = 0.25 the sample x[1] = 0.1 when
quantized gives x̂[1] = 0 since 0 < 0.1 < ∆. If we multiply the signal by 3 so that the value
3x[1] = 0.3 would give x̂[1] = 1 since ∆ < 0.3 < 2∆. Since the second output is not the first
output multiplied by 3, the quantizer is not linear.
We saw before that the sampler is linear but time-varying, while the quantizer is time-invariant
but non-linear, so the A/D converter which is composed of these systems is not LTI.
dy(t) dx(t)
+ y(t) = 2x(t) +
dt dt
This equation is discretized by approximating the derivatives for a signal ρ(t) as
Solution
Discretization of differential equation:
x[n] x[n − 1]
D
y[n]
+
9.4 An LTI discrete-time system has the impulse response h[n] = (−1)n u[n]. Use the convolution
sum to compute the output response y[n], n ≥ 0 when the input is x[n] = u[n] − u[n − 3] and the
initial conditions are zero. In particular, find the following values of the output y[n], 0 ≤ n ≤ 4.
Answers: y[0] = 1, y[1] = 0, y[2] = 1, y[3] = −1, y[4] = 1.
Solution
The convolution sum is
n
X
y[n] = x[m]h[n − m]
m=0
y[0] = x[0]h[0] = 1
y[1] = x[0]h[1] + x[1]h[0] = −1 + 1 = 0
y[2] = x[0]h[2] + x[1]h[1] + x[2]h[0] = 1 − 1 + 1 = 1
y[3] = x[0]h[3] + x[1]h[2] + x[2]h[1] + x[3]h[0] = −1 + 1 − 1 + 0 = −1
y[4] = x[0]h[4] + x[1]h[3] + x[2]h[2] + x[3]h[1] + x[4]h[0] = 1 − 1 + 1 + 0 + 0 = 1
y[n] = x[0]h[n] + x[1]h[n − 1] + x[2]h[n − 2] = (−1)n + (−1)n−1 + (−1)n−2 = (−1)n−2 n ≥ 5
9.5 An LTI system represented by the difference equation y[n] = −y[n − 1] + x[n], n ≥ 0, is initially
at rest. The input of the system is x[n] and the output y[n].
(a) Using the difference equation find the system response y[n] for n ≥ 0 when the input is
1 n = 0, 1, 2
x[n] =
0 otherwise
Answers: y[n] = 0 for n < 0, y[0] = 1, y[1] = 0, y[n] = (−1)n for n ≥ 2; y1 [n] = y[n] + y[n − 3].
Solution
(a) According to the difference equation the output is
or
y[n] = 0 n<0
y[0] = 1
y[1] = 0
y[n] = (−1)n n≥2
9.6 The output of a discrete-time system is y[n] = w[n]x[n] where x[n] is the input, and w[n] =
u[n] − u[n − 5] is a rectangular window.
(a) The input is x[n] = 4 sin(πn/2), −∞ < n < ∞, determine if x[n] is periodic, and if so
indicate its fundamental period N0 .
(b) Let z[n] = w[n]x[n − N ] for what values of N is y[n] = z[n]? According to this, is the
system time-invariant? Explain.
(c) Suppose that we let x1 [n] = x[n]u[n], i.e., x1 [n] is the causal component of x[n]. Determine
the outputs due to x1 [n] and to x1 [n − 4]. Again according to this, is the system time-
invariant? Explain.
Solution
(a) The discrete frequency of x[n] is ω0 = π/2 which can be written as 2πm/N0 with m = 1,
N0 = 4. So x[n] is periodic of fundamental period N0 = 4.
(b) Since x[n − kN0 ] = x[n] for any integer k and the fundamental period N0 = 4, then taking
N = kN0 we will have that
so that delaying the input by a multiple of N0 does not change the output. The output in
this case is just the period of x[n] between 0 and 4. However, that is not so if the delay is
different from a multiple of the fundamental period, thus the system is time-varying.
(c) The output corresponding to x1 [n] = x[n]u[n] is
which is not the previous result shifted, so that the system is time-varying.
9.7 Express the signal u12 [m, n] with support in the first and second quadrant in terms of one-
dimensional unit-step functions and determine if it is separable.
Solution
We have
1 −∞ < m < ∞, n ≥ 0
u12 [m, n] =
0 otherwise
= u1 [m, n] + u2 [m, n] − u1 [0, n] = u[m]u[n] + u[−m]u[n] − u[n]
= (u[m] + u[−m] − 1)u[n]
i.e., it is separable.
9.8 Although few signals encountered in practice are separable, any signal x[m, n], 0 ≤ m ≤ M −
1, 0 ≤ n ≤ N − 1 and zero otherwise, can be written as as the sum of finite number of separable
sequences
NX−1
x[m, n] = xi1 [m]xi2 [n].
i=0
Choose xi1 [m] and xi2 [n] to obtain this representation and show it can be expressed in terms of
δ[m, n].
Solution
A simple way to obtain such a representation is by letting the signal be expressed as the sum of
isolated columns. This is done by choosing for 0 ≤ i ≤ N − 1
so that
N
X −1
x[m, n] = x[m, i]δ[n − i]
i=0
and thus
N
X −1 N
X −1 M
X −1
x[m, n] = x[m, i]δ[n − i] = x[k, i] δ[m − k]δ[n − i]
| {z }
i=1 i=0 k=0
δ[m−k,n−i]
9.9 Suppose both the input x[m, n] and the impulse response h[m, n] of a LSI system are separable.
(a) Obtain a simpler expression for the output y[m, n] of the system.
(b) Suppose x[m, n] = (u[m] − u[m − 2])(u[n] − u[n − 2]) and h[m, n] = (u[m] − u[m − 2])(u[n] −
u[n − 2]) where u[n] is the one-dimensional unit-step signal. Find y[m, n]
Solution
(a) Letting x[m, n] = x1 [m]x2 [n] and h[m, n] = h1 [m]h2 [n] we have that
∞
X ∞
X
y[m, n] = x1 [k]x2 [`]h1 [m − k]h2 [n − `]
k=−∞ `=−∞
X∞ ∞
X
= x1 [k]h1 [m − k] x2 [`]h2 [n − `]
k=−∞ `=−∞
= (x1 ∗ h1 )[m] (x2 ∗ h2 )[n]
Find the rest of the sequence for 0 ≤ n ≤ 5 and plot it using function stem.
Answers: x[3] = 2, x[4] = 3, x[5] = 5.
Solution
The equation
with initial conditions x[0] = 0, x[1] = 1 and x[2] = 2 can be solved recursively for n ≥ 3 as
follows
% Pr 9_10
clear all; clf
x(1)=0; x(2)=1; x(3)=2; % initial conditions
for m=4:50,
x(m)=x(m-1)+x(m-3);
end
n=0:49;
figure(1)
stem(n,x); title(’x[n]’); xlabel(’n’)
grid
0
0 5 10 15 20 25 30 35 40 45 50
n
9.11 Generation of periodic discrete-time signal — Periodic signals can be generated by obtaining
a period and adding shifted versions of this period. Suppose we wish to generate a train of
triangular pulses. A period of the signal is x[n] = 0.5(r[n] − 2r[n − 2] + r[n − 4]) where r[n]
is the discrete-time ramp signal. Obtain the values of x[n] for 0 ≤ n ≤ 4 and write a script to
generate x[n]and the periodic signal y[n]
∞
X
y[n] = x[n − 4k]
k=−∞
Solution
(a) A period of a train of pulses is
0 n<0
r[n] = 0.5n 0≤n≤2
x[n] = 0.5(r[n] − 2r[n − 2] + r[n − 4]) =
r[n] − 2r[n − 2] = 0.5(4 − n) 2<n≤4
r[n] − 2r[n − 2] + r[n − 4] = 0 n>4
The following MATLAB script will generate x[n] for n ≥ 0 (it is zero before) and the periodic
signal y[n].
% Pr 9_11
x=zeros(1,5);
for m=1:3,
x(m)=m-1;
end
for m=3:5,
x(m)=4-(m-1)
end
x=x/2;
n=0:4;
figure(1)
subplot(211)
stem(n,x); axis([0 4 -0.1 1.1])
hold on; plot(n,x,’:r’); hold off; title(’x[n]’)
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3 3.5 4
y[n]
1
0.8
0.6
0.4
0.2
0
0 2 4 6 8 10 12 14 16 18
n
Figure 9.3: Problem 11: One period, 4 periods of y[n].
9.12 Envelope modulation — In the generation of music by computer, the process of modulation is
extremely important. When playing an instrument, the player typically does it in three stages:
(1) rise time, (2) sustained time and (3) decay time. Suppose we model these three stages as an
envelope continuous-time signal given by
1 1
e(t) = [r(t) − r(t − 0.025)] − [r(t − 200) + r(t − 300)]
0.025 0.1
(a) For a simple tone x(t) = cos(2π/T0 t), the modulated signal is y(t) = x(t)e(t). Find the
period T0 so that 100 cycles of the sinusoid occur for the duration of the envelope signal.
(b) Simulate in the modulated signal using the value of T0 and a simulation sampling time
of Ts = 0.1. Plot y(t) and e(t) (discretized with the sampling period Ts ) and listen to the
modulated signal using sound.
Solution
(a) If the duration of e(t) is Te = 300Ts and for 100 cycles of the sinusoid 100T0 then T0 =
Te /100 = 3Ts .
(b) If we let t = nTs , then the sampled envelope is
The following script generates the window e[n] and the modulated signal y[n] = e[n]x[n]
% Pr 9_12
clear all;clf
Ts=0.1;
t1=0:Ts:3;
t2=20:Ts:30;
e=[t1/3 ones(1,200-30) 3-t2/10]
t=0:Ts:30+Ts;
figure(1)
plot(t,e,’r’);grid
n=0:length(e)-1;
x=cos(2*pi*n/3);
hold on
plot(t,x.*e,’b’)
hold off; axis([0 30 -1.2 1.2]);
xlabel(’nT_s’);ylabel(’x(nT_s)e(nT_s)’)
0.8
0.6
0.4
0.2
x(nTs)e(nTs)
0
−0.2
−0.4
−0.6
−0.8
−1
0 5 10 15 20 25 30
nTs
Figure 9.4: Problem 12: window e(nTs ) and modulated cosine y(n) = x(nTs )e(nTs ) for Ts = 0.1.
9.13 LTI and convolution sum — The impulse response of a discrete-time system is
h[n] = (−0.5)n u[n].
(a) If the input of the system is x[n] = δ[n] + δ[n − 1] + δ[n − 2], use the linearity and time-
invariance of the system to find the corresponding output y[n].
(b) Find the convolution sum corresponding to the above input, and show that your solution
coincides with the output y[n] obtained above.
(c) Use the function conv to find the output y[n] due to the given input x[n]. Plot x[n], h[n]
and y[n] using .
Solution
(a) Using the superposition property of LTI systems and that h[n] is the response to δ[n] gives
y[n] = h[n] + h[n − 1] + h[n − 2] = (−0.5)n u[n] + (−0.5)n−1 u[n − 1] + (−0.5)n−2 u[n − 2]
1 n=0
0.5 n=1
=
0.75 n=2
3(−0.5)n n ≥ 3
since the multiplication by the delta functions gives zero except when their arguments are zero,
or k = 0, k = 1 and k = 2. This result coincides with the one in part (a).
(c) The following script is used to compute the convolution sum of the given x[n] and h[n]:
% Pr 9_13
clear all; clf
x=[1 1 1 zeros(1,98)];
n=0:100;
h=(-0.5).ˆn;
y=conv(h,x);
n1=0:20;
figure(1)
subplot(311)
stem(n1,x(1:21)); ylabel(’x[n]’);grid
subplot(312)
stem(n1,h(1:21)); ylabel(’h[n]’);grid
subplot(313)
stem(n1,y(1:21)); ylabel(’y[n]’);grid;xlabel(’n’)
x[n]
0.5
0
0 2 4 6 8 10 12 14 16 18 20
0.5
h[n] 0
−0.5
0 2 4 6 8 10 12 14 16 18 20
0.5
y[n]
−0.5
0 2 4 6 8 10 12 14 16 18 20
n
Figure 9.5: Problem 13: Input x[n], impulse response h[n] and convolution sum y[n].
(a) Generate a signal x[m, n] with period p[n, m] and rectangular periods N1 = 3 and N2 = 2,
i.e., x[m, n] = x[m + 3, n] = x[m, n + 2]. Use the functions imshow and stem3 to display the
signals.
(b) Use p[m, n] to create an array y[m, n] with periodic matrix
3 1
N=
0 2
i.e.,
y[m, n] = y[m + 3, n + 1] = y[m, n + 2],
and plot several periods as an image. Is this signal rectangular periodic? If so indicate the
periods.
Solution
(a) The array will be a stacking of p[m, n] horizontally and vertically to generate x[m, n] which
can be visualized for m ≥ 0, n ≥ 0 for 9 periods in first quadrant
(b) In this case the [0, 0] value appears in [3, 1] and in [0, 2] as indicated by the periodicity
matrix. The overall array can be seen as in the first quadrant
The generated y[m, n] is also periodic of periods N1 = 6 and N2 = 2. The period for these
values contains 12 points instead of the 9 for the given p[m, n].
% Pr 9_14
clear all;clf
p=[0 0 0;1 0 0]
X=[p;p;p];
X=[X X X;X X X]
Y=zeros(4,9);
Y(1,4)=1; Y(2,1)=1;Y(2,7)=1;
Y(3,4)=1;Y(4,1)=1;Y(4,7)=1;
Y=[Y;Y;Y]
figure(1)
subplot(121)
imshow(X)
hold on
stem3(X,’filled’);title(’x[m,n]’)
hold off
subplot(122)
imshow(Y)
hold on
stem3(Y,’filled’);title(’y[m,n]’)
hold off
x[m,n] y[m,n]
Figure 9.6: Rectangular periodic x[m, n] and block periodic y[m, n].
(x[m, n] − M )
y[m, n] = α + β
(N − M )
where the output is y[m, n] and the input is x[m, n] and α, β, M and N are parameters to be
determined for a certain image. Let the input x[m, n be the image clown in a gray scale.
(a) Suppose that we wish to change the pixel values of the output y[m, n] so that they range
from 6 to 256. Determine the values of M and N from the input image, and the parameters
α and β so that the desired scaling is possible. Use the functions colormap and imagesc to
plot the output in a gray scale with a range [0, 255].
(b) Let α = 9 and β = 246, and the input signal x[m, n] be the test image clown. Is this scaling
system linear? Explain why or why not. Plot the resulting image with a range of [0, 255]
for the pixel values.
(c) Compute and plot the histogram of the two images and verify the range of values of the
pixels of the two images.
Solution
(a) The values needed to obtain the desired image are M = 1 and N = 81 corresponding to
the minimum and maximum pixel values of clown and α = 6 and β = 200.
(b) The equation relating the input and the output for the given values of α and β is
(x[m, n] − 1)
y[m, n] = 10 + 246
80
This system is non-linear because of the term 10 which is a bias that does not change with
changes in the input.
% Pr 9_15
clear all;clf
load clown; x=X;
M1=min(x(:))
M2=max(x(:))
for m=1:200.
for n=1:320,
y(m,n)=9+246*(x(m,n)-M1)/(M2-M1);
end
end
M3=min(y(:))
M4y=max(y(:))
figure(1)
colormap(gray)
subplot(221)
imagesc(x,[0 256])
subplot(222)
hist(x,100)
subplot(223)
imagesc(y,[0 256])
subplot(224)
hist(y,100)
The Z–transform
10.1 The sign function extracts the sign of a real valued signal, i.e.,
−1
x[n] < 0
s[n] = sign(x[n]) = 0 x[n] = 0
1 x[n] > 0
(a) Let s[n] = s1 [n]+s2 [n], where s1 [n] is causal and s2 [n] anti–causal; find their Z-transforms
and indicate the corresponding ROCs.
Solution
(a) We have s[n] = u[n] − u[−n], the Z-transform of s1 [n] = u[n] is
∞
X 1
S1 (z) = z −n =
n=0
1 − z −1
provided that |z −1 | < 1, thus |z| > 1 is the region of convergence for S1 (z). The Z-transform
of s2 [n] = −u[−n] is given by
0
X ∞
X −1
S2 (z) = − z −n = − zm =
n=−∞ m=0
1−z
1
Chaparro-Akan — Signals and Systems using MATLAB 10.2
1 1 1 + z −1
S(z) = S1 (z) + S2 (z) = −1
− =
1−z 1−z 1 − z −1
to converge is that |z| > 1 and that |z| < 1 simultaneously, which is not possible. Since there
is no region of convergence for S(z), the Z-transform of s[n] does not exist.
10.2 An analog pulse x(t) = u(t) − u(t − 1) is sampled using a sampling period Ts = 0.1.
(a) Obtain the discrete-time signal x(nTs ) = x(t)|t=nTs and plot it as a function of nTs
(b) If the sampled signal is represented as an analog signal as
N
X −1
xs (t) = x(nTs )δ(t − nT s)
n=0
determine the value of N in the above equation. Compute the Laplace transform of the
sampled signal. i.e., Xs (s) = L[xs (t]
(c) Determine the Z-transform of x(nTs ) or X(z). Indicate how to transform Xs (s) into X(z).
Solution
(a) If Ts = 0.1 the discrete-time signal is
(
1 0 ≤ 0.1n ≤ 1 or 0 ≤ n ≤ 10
x(0.1n) = [u(t) − u(t − 1)] |t=0.1n =
0 otherwise
10.3 When finding the inverse Z-transform of a function with z −1 terms in the numerator, z −1 can
be thought of a delay operator to simplify the calculation. For
1 − z −10
X(z) =
1 − z −1
(a) Use the Z-transform of u[n] and the properties of the Z-transform to find x[n].
(b) Use your result for x[n] above to express X(z) as a polynomial in negative powers of z.
(c) Find the poles and zeros of X(z) and plot them on the Z–plane. Is there a pole or zero at
z = 1? Explain.
Solution
(a) Writing X(z) as
1 z −10
X(z) = −
1 − z −1 1 − z −1
since the inverse Z-transform of the first term is u[n], then the inverse of the second is −u[n −
10] given that z −10 indicates a delay of 10 samples. Thus,
(
1 0≤n≤9
x[n] = u[n] − u[n − 10] =
0 otherwise
(b) Although X(z) has been shown as a ratio of two polynomials, using the above representa-
tion of x[n] its Z-transform is
X(z) = 1 + z −1 + · · · + z −9
which is obtained by finding that the zeros of X(z) are values zk10 = 1 or zk = ej2πk/10 for
k = 0, · · · , 9. For k = 0 the zero is z0 = 1, which cancels the pole at 1.
1 − αM z −M
H(z) = for any value of α?
1 − αz −1
ii. Let M = 3, 0 ≤ α < 1 write H(z) as a polynomial, and then show that it equals
1 − α3 z −3
. Determine the region of convergence of H(z).
1 − αz −1
Solution
1 − α3 z −3 z 3 − α3
H(z) = =
1 − αz −1 z 2 (z − α)
so the zeros of H(z) are zk = αej2πk/3 for k = 0, 1, 2, and the poles z = 0 of order
2, and z = α. The pole z = α is cancelled by the zero at k = 0 so ROC is the whole
Z–plane, except z = 0.
(b) i. If we let
0.5 n=0
h1 [n] = 0.5n 1≤n≤N −1
0 otherwise
noticing that the initial expression for H1 (z) is a polynomial in z −1 then the ROC
is |z| > 0. The final expression being rational indicates that if z −1 = 2 or z = 1/2
is a zero and a pole and they cancel so that only poles at the origin remain. The
Z-transform
3
X 0.5 + 0.25z − 0.54 z 4
Z[h1 [−n]] = H1 (1/z) = 0.5n z n − 0.5 =
n=0
1 − 0.5z
with ROC the whole Z–plane as it is a polynomial in z. Verify that the pole at z = 2
of the final expression for Z[h1 [−n]] is cancelled by a zero at the same place. Thus
10.5 Consider a discrete-time LTI system represented by the difference equation with the given
initial condition
Solution
We have
Y (z) 2(1 − z −1 )
H(z) = =
X(z) 1 + 0.5z −1
Y (z) = H(z)(Z[u[n]] + Z[0.5n cos(πn)u[n]])
only the first of the above terms would give steady state due to its pole at z = 1, so
1 2(1 − z −1 ) 2
H(z) = =
1 − z −1 (1 − z −1 )(1 + 0.5z −1 ) 1 + 0.5z −1
which gives a transient, so yss [n] = 0. Since system is stable, changing the initial condition
will not change the steady state–response.
10.6 Determine the impulse response h[n] of the feedback system shown in Fig. 10.1. Determine if
the system is BIBO stable.
Answer: System is not BIBO stable.
e[n]
x[n] Delay y[n]
+
−
Moreover, the pole z = 1 is on the unit disc and as such the system is not BIBO stable.
(a) The input and the output of a LTI discrete-time system are
(
1 n = 0, 1 1 n = 0, 3
Input x[n] = Output y[n] = 2 n = 1, 2
0 otherwise
0 otherwise
1 − 0.5z −1
H(z) = |z| > 0.5
1 − 0.25z −2
i. Is this system causal? Explain.
ii. Determine the impulse response of the system.
(c) The transfer function of an LTI system is
z+1
H(z) = |z| > 1
z(z − 1)
What are the values of h[0], h[1] and h[1000] of the impulse response.
Solution
(a) We have
Y (z) = 1 + 2z −1 + 2z −2 + z −3
X(z) = 1 + z −1
Y (z)
H(z) = = 1 + z −1 + z −2
X(z)
after division.
(b) i. Yes, because of the ROC.
ii. There is a pole-zero cancellation
1 − 0.5z −1
H(z) =
(1 − 0.5z −1 )(1 + 0.5z −1 )
H(z) = z −1 + 2z −2 + 2z −3 + · · · + 2z −10000 + · · ·
10.8 The impulse response of a LTI discrete-time system is h[n] = u[n] − u[n − 3]. If the input to
this system is
0 n<0
x[n] = n n = 0, 1, 2
1 n≥3
(a) Find the output y[n] of this system by calculating graphically the convolution sum.
(b) Express x[n] in terms of basic functions and determine X(z).
(c) Use the Z-transform to find the output y[n].
Answers: x[n] = δ[n − 1] + 2δ[n − 2] + u[n − 3]; y[n] = x[n] + x[n − 1] + x[n − 2]
Solution
y[0] = 0 y[1] = 1
y[2] = 3 y[3] = 4
y[4] = 4 y[n] = 3 n ≥ 5
thus
so that
10.9 The transfer function of an RLC circuit is H(s) = Y (s)/X(s) = 2s/(s2 + 2s + 1).
(a) Obtain the ordinary differential equation with input x(t) and output y(t). Approximating
the derivatives by differences (let T = 1) obtain the difference equation that approximates
the differential equation.
(b) Let the input be a constant source, so that x[n] = u[n], and let the initial conditions for
the difference equation be zero. Solve the difference equation using the Z-transform.
Solution
Y (s) 2s
H(s) = = 2
X(s) s + 2s + 1
(y[n] − 2y[n − 1] + y[n − 2]) + 2(y[n] − y[n − 1]) + y[n] = 2(x[n] − x[n − 1])
y[n] − y[n − 1] + 0.25y[n − 2] = 0.5(x[n] − x[n − 1])
Using Z-transform
10.10 Suppose we cascade a “differentiator” and a “smoother” systems characterized by the follow-
ing input/output equations
(a) If x[n] = u[n] and the initial conditions for the smoother are zero, find the output of the
overall system y[n].
(b) If x[n] = (−1)n , −∞ < n < ∞, find the steady-state response yss [n] of the overall system
yss [n].
Answers: The overall transfer function is H(z) = (1/3)(1−z −1 )/(1−z −1 /3); yss [n] = 0.5(−1)n .
Solution
Let Hs (z) and H(z) be the transfer functions of the smoother and of the overall system.
(a) If x[n] = u[n] then w[n] = δ[n] and y[n] = hs [n]. Thus, from
Y (z) 1/3
Hs (z) = = ⇒ y[n] = hs [n] = (1/3)n+1 u[n]
W (z) 1 − z −1 /3
10.11 Suppose we are given a finite-length sequence h[n] (it could be part of an infinite-length im-
pulse response from a discrete system that has been windowed) and would like to obtain a
rational approximation for it. This means that if H(z) = Z[h[n]], a rational approximation of
it would be H(z) = B(z)/A(z), from which we get
H(z)A(z) = B(z)
Letting
M
X −1 N
X −1
B(z) = bk z −k , A(z) = 1 + ak z −k
k=0 k=1
for some choice of M and N , equations from H(z)A(z) = B(z) should allow us to find the
M + N − 1 coefficients {ak bk }.
(a) Find a matrix equation that would allow us to find the coefficients of B(z) and A(z).
(b) Let h[n] = 0.5n (u[n] − u[n − 101]) be the sequence we wish to obtain a rational approx-
imation and let B(z) = b0 while A(z) = a0 + a1 z −1 , find the equations to solve for the
coefficients {b0 , a0 , a1 }. Use MATLAB to solve the equations.
PN −1
Answers: h[m] + k=1 ak h[m − k] = bm where h[m] = 0 for m < 0 and bm = 0 for m < 0 and
m > M − 1.
Solution
(a) The equation H(z)A(z) = B(z) using convolution and that a0 = 1
N
X −1
h[m] + ak h[m − k] = bm
k=1
Solving first for the denominator coefficients (the bottom equations above), we get a matrix
equation
h[M ] h[M − 1] ··· h[M − N + 1] 1 0
h[M + 1] h[M ] ··· h[M − N + 2] a1 0
.. .. .. .. .. = ..
. . . . . .
h[M + N − 2] h[M + N − 3] · · · h[M − 1] aN −1 0
To solve these equations we move the first column of the matrix on the left to replace the zero
vector on the right, using the fact that a0 = 1 is known.
Once the denominator coefficients are found, the following set of equations are solved for the
numerator coefficients:
h[0] 0 ··· 0 1 b0
h[1] h[0] ··· 0 a1 b1
. .. .. .. .. = ..
.
. . . . . .
h[M − 1] h[M − 2] · · · h[M − N ] aN −1 bM −1
h[0]a1 = −h[1]
b0 1
Ĥ(z) = =
1 + a1 z −1 1 − 0.5z −1
with impulse response ĥ[n] = 0.5n u[n] which approximates very well the given impulse re-
sponse.
Solution
(a) The support of x[m, n] coincides with that of u12 [m, n] which is the first and second-
quadrant in the (m, n) plane.
(b) We have
∞
" ∞
# " ∞ −1
#
X X 1 X X
X(z1 , z2 ) = 0.5m z1−m 0.5n z2−n = 0.5m z1−m + 0.5m z1−m
m=−∞ n=0
1 − 0.5z −1 m=0 m=−∞
| {z }
1/(1−0.5z −1
1 1 1
= + −1
1 − 0.5z −1 1 − 0.5z1−1 1 − 2z1−1
ROC: |z2 | > 0.5, (|z1 | > 0.5, and |z1 | < 0.5)
and because the last two are contradictory the Z-transform does not exist.
10.13 The region of convergence of H(z1 , z2 ) = 1/(1 − (z1−1 + z2−1 ) is |z1−1 | + |z2−1 | < 1. Determine if
the unit bidisc is included in this ROC and sketch a possible ROC in a |z1 | × |z2 | plane.
Solution
If |z1 | = |z2 | = 1, the bidisc, their sum is 2 which is bigger than 1, so the bidisc is not in the
ROC. The ROC can be expressed as
|z1 |
|z2 | >
|z1 | − 1
10.14 For what values of α is the filter with the following transfer function not BIBO stable?
1
H(z1 , z2 ) =
1 + 0.5z −1 + αz2−1
Solution
The condition for stability
A(z1 , 1) = 1 + 0.5z1−1 + α
has as root z1 = −0.5/(1 + α) and if |z1 | = |0.5/(1 + α)| ≥ 1, i.e., on or outside the unit
circle, the filter is unstable. So for any value of α that makes |1 + α| ≤ 0.5, or any α such that
−1.5 ≤ α ≤ −0.5
10.15 Fibonacci sequence generation — Consider the Fibonacci sequence generated by the differ-
ence equation f [n] = f [n − 1] + f [n − 2], n ≥ 0 with initial conditions f [−1] = 1, f [−2] = −1.
(a) Find the Z-transform of f [n], or F (z). Find the poles φ1 , and φ2 and the zeros of F (z).
How are the poles connected? How are they related to the ”golden-ratio”?
(b) The Fibonacci difference equation has zero input, but its response is a sequence of ever-
increasing integers. Obtain a partial fraction expansion of F (z) and find f [n] in terms
of the poles φ1 and φ2 , and show that the result is always integer. Use MATLAB to
implement the inverse in term of the poles.
Solution
(a) Using the difference equation and the given initial conditions we get the values of the
Fibonacci series:
or the sequence {0, 1, 1, 2, 3, 5, 8, 13, 21, 34, · · · } for n ≥ 0. The Z-transform of f [n] is then
which is not in a closed-form. To get that, we use the Z-transform shift property in the given
difference equation:
z −1
F (z) =
1− z −1
− z −2
We have obtained the Z-transform of the Fibonacci sequence even though we know that the
sequence blows up, and that the region of convergence for this sequence is not obvious from
the development.
The value φ1 = 1.6183 is the golden ratio which in ancient Greece was considered the most
pleasing ratio for many designs. The letter φ has been used in honor of the Greek sculptor
Phidias who is said to have used it in his work.
(b) Notice that φ1 + φ2 = 1. The partial fraction expansion is
A B
F (z) = −1
+
1 − φ1 z 1 − φ2 z −1
1
f [n] = √ (φn1 − φn2 )u[n]
5
This formula was first presented by Leonhard Euler in 1765 as a formula to obtain the Fi-
bonacci numbes. Let us verify that this equation gives the Fibonacci numbers, f [0] = 0,
f [1] = √1 (φ1 − φ2 ) = 1, and so on. For the region of convergence |z| > φ1 = 1.6183, the
5
system represented by the Fibonacci difference equation is not BIBO stable as one of its poles
is outside the unit circle.
Using the poles φi we compute the Fibonacci sequence. The obtained values do not coincide
exactly with the sequence, only if we find integers that smaller using floor we obtain it.
% Pr. 10.15
clear all; clf
phi1=1.6183; phi2=-0.6183;
n=0:100;
f=(phi1.ˆn-phi2.ˆn)/sqrt(5);
figure(1)
stem(n(1:10),f(1:10));grid;title(’Fibonacci sequence’);xlabel(’n’)
floor(f(1:10))
0 1 1 2 3 5 8 13 21 34
30
25
20
15
10
0
0 1 2 3 4 5 6 7 8 9
n
10.16 Generation of discrete-time sinusoid — Given that the Z-transform of a discrete-time cosine
A cos(ω0 n)u[n] is
A(1 − cos(ω0 )z −1 )
1 − 2 cos(ω0 )z −1 + z −2
(a) Use the given Z-transform to find a difference equation whose output y[n] is a discrete-
time cosine A cos(ω0 n) and input x[n] = δ[n]. What should you use as initial conditions?
(c) Indicate how to change your previous algorithm to generate a sine function y[n] =
2 sin(πn/2)u[n]. Use MATLAB to find y[n], and to plot it.
Solution
(a) Thinking of the given Z-transform as a transfer function with x[n] as the input and y[n] as
the output we have that
with zero initial conditions and input x[n] = δ[n] that will generate the cosine function.
(b) We have that ω0 = π/2 so that cos(ω0 ) = 0 giving the difference equation
y[0] = 2
y[1] = 0
y[2] = −2
y[3] = 0
..
.
or y[n] = 2 cos(πn/2)u[n].
y1 [0] = 0
y1 [1] = 2
y1 [2] = 0
y1 [3] = −2
..
.
i.e., y1 [n] = y[n − 1], so that the above difference equation shifted by 1 gives
which gives
y1 [0] = −y[−2] + 2 = 0
y1 [1] = −y[−1] + 0 = 2
y1 [2] = y1 [0] = 0
y1 [3] = −y[1] = 2
..
.
2z −1
Y (z) =
1 + z −2
which is the Z-transform of the desired signal.
The following script computes and plots the two sinusoids. To generate the sine, we compute
the first two values using the initial conditions and the input, and then use the difference
equation to obtain the rest using a for loop.
% Pr. 10.16
clear all; clf
% generation of 2cos(pi n/2)
b=2;
a=[1 0 1];
x=[1 zeros(1,100)];
y=filter(b,a,x);
n=0:length(y)-1;
figure(1)
subplot(211)
stem(n,y); axis([0 100 -2.2 2.2]); grid; ylabel(’y[n]’)
hold on; plot(n,y,’:r’); hold off
% generation of 2sin(pi n/2)
y1(1)=0;y1(2)=2;
for m=3:101,
y1(m)=-y1(m-2);
end
subplot(212)
stem(n,y1); axis([0 100 -2.2 2.2]);grid; ylabel(’y_1[n]’);xlabel(’n’)
hold on; plot(n,y1,’:r’); hold off
1
y[n]
−1
−2
0 10 20 30 40 50 60 70 80 90 100
1
y1[n]
−1
−2
0 10 20 30 40 50 60 70 80 90 100
n
Figure 10.3: Problem 16: Generation of 2 cos(πn/2) (top) and 2 sin(πn/2) (bottom).
10.17 Inverse Z-transform — We are interested in the unit-step solution of a system represented by
the following difference equation y[n] = y[n − 1] − 0.5y[n − 2] + x[n] + x[n − 1]
Solution
(a) To find the unit–step response, the input is x[n] = u[n] and the initial conditions are zero.
We then have
1
Y (z)[1 − z −1 + 0.5z −2 ] = (1 + z −1 )
1 − z −1
1 + z −1
Y (z) =
(1 − z −1 )(1 − 0.5(1 + j)z −1 )(1 − 0.5(1 − j)z −1 )
A B B∗
Y (z) = + +
1 − z −1 1 − 0.5(1 + j)z −1 1 − 0.5(1 − j)z −1
A = Y (z)(1 − z −1 )|z−1 =1 = 4
o
B = Y (z)(1 − 0.5(1 + j)z −1 )|z−1 =0.5(1+j) = 1.58e−j161.5
The following script is used to verify the above result using MATLAB
%% Pr. 10.17
clear all; clf
N=[1 1 0 0];
D=conv([1 -1],[1 -1 0.5])
[r,p,k]=residuez(N,D)
abs(r(2))
angle(r(2))*180/pi
n=0:49;
y=r(1)*(p(1).ˆn)+r(2)*(p(2).ˆn)+r(3)*(p(3).ˆn);
figure(1)
subplot(311)
zplane(N,D)
subplot(312)
stem(n,y); ylabel(’y[n]’)
% verification
y1=4+3.16*(0.707).ˆn.*cos(pi.*n/4-161.5*pi/180);
subplot(313)
stem(n,y1); ylabel(’y_1[n]’); xlabel(’n’)
1
Imaginary Part
0 2
−1
−5 −4 −3 −2 −1 0 1 2 3 4 5
Real Part
6
4
y[n]
0
0 5 10 15 20 25 30 35 40 45 50
4
y1[n]
0
0 5 10 15 20 25 30 35 40 45 50
n
Figure 10.4: Problem 17: Verification of unit-step response using residues and analytic solution. Poles and
zeros of Y (z); unit–step response found using the partial fraction expansion with MATLAB (middle) and
analytically (bottom).
10.18 Prony’s Rational Approximation — The Pade approximant provides an exact matching of
M + N − 1 values of h[n], where M and N are the orders of the numerator and denominator
of the rational approximation. But there is no method for choosing the numerator and de-
nominator orders, M and N . Also, there is no guarantee on how well the rest of the signal is
matched. Prony’s rational approximation considers how well the rest of the signal is approxi-
mated. Let h[n] = 0.9n u[n], be the impulse response we wish to find a rational approximation.
Take the first 100 values of this signal as the impulse response.
(a) Assume the order of the numerator and the denominators are equal M = N = 1, use the
MATLAB function prony to obtain the rational approximation, and then use filter to verify
that the impulse response of the rational approximation is close to the given 100 values.
Plot the error between h[n] and the impulse response of the rational approximation for
the first 100 samples. Plot the poles and zeros of the rational approximation and compare
them to the poles and zeros of H(z) = Z(h[n].
(b) Suppose that h[n] = (h1 ∗ h2 )[n], i.e., the convolution of h1 [n] = 0.9n u[n] and h2 [n] =
0.8n u[n]. Use again Prony to find the rational approximation when the first 100 values
of h[n] are available. Use conv from MATLAB to compute h[n]. Compare the impulse
response of the rational approximation to h[n]. Plot the poles and zeros of H(z) = Z(h[n])
and of the rational approximation.
(c) Consider the h[n] given above, and perform the Prony approximation using orders M =
N = 3, explain your results. Plot poles and zeros.
Solution
(a) The following script obtains the Prony rational approximation of the first 100 value of
h[n] = 0.9n u[n].
% Pr. 10.18
clear all; clf
n=[0:100];
h=0.9.ˆn;
[b,a]=prony(h,1,1); % rational approximation of order (1,1) of h
delta=[1 zeros(1,100)];
hd=filter(b,a,delta); % approximate sequence
figure(1)
subplot(211)
plot(n,h,’o’,n,hd,’+’)
xlabel(’n’)
ylabel(’h[n],hd[n]’)
grid
p=roots(a) % poles
z=roots(b) % zeros
subplot(223)
stem(n,hd-h); ylabel(’error’)
subplot(224)
zplane(b,a)
h[n],hd[n] 0.8
0.6
0.4
0.2
0
0 10 20 30 40 50 60 70 80 90 100
n
−15
x 10
2.5 1
2
0.5
Imaginary Part
1.5
error
0
1
−0.5
0.5
−1
0
0 20 40 60 80 100 −1 −0.5 0 0.5 1
Real Part
Figure 10.5: Problem 18: Prony approximation, of order (1, 1), of finite segment of h[n] = 0.9n u[n] .
The poles and zeros corresponding to the actual and the approximate are very close.
(b) The convolution of h1 [n] = 0.9n u[n] and h2 [n] = 0.8n u[n] is done using the MATLAB func-
tion conv. The Prony approximation of h1 [n] ∗ h2 [n] gives an approximation to H1 (z)H2 (z) for
the correct orders. Using prony of order (2, 2) (in general this has to be guessed, but here we
have the advantage of knowing the sequences, a big advantage!).
p=roots(a) % poles
z=roots(b) % zeros
subplot(212)
zplane(b,a)
The following are the results obtained. Notice that the poles/zeros are those of H1 (z) and
H2 (z) or the product H(z) = H1 (z)H2 (z).
2.5
2
h[n],hd[n]
1.5
0.5
0
0 10 20 30 40 50 60 70 80 90 100
n
0.5
Imaginary Part
−0.5
−1
−3 −2 −1 0 1 2 3
Real Part
Figure 10.6: Problem 18: Prony approximation of convolution h1 [n] and h2 [n].
(c) When the chosen order is higher than the actual one, it is possible to have a pole-zero can-
celations. For the h[n] = [h1 ∗h2 ][n] when using prony instead of a second-order approximation
we choose a third-order, the results has a pole/zero cancellation.
[b1,a1]=prony(h,3,3)
p=roots(a1)
z=roots(b1)
subplot(211)
zplane(z,p)
delta=[1 zeros(1,120)];
hd=filter(b1,a1,delta);
M=length(h);
n=[0:ceil(M/2)-1];
subplot(212)
plot(n,h(1:ceil(M/2)),’o’,n,hd(1:ceil(M/2)),’+’)
b1 =
1.0000 0.0000 -0.0000 -0.0000
a1 =
1.0000 -1.7000 0.7200 0.0000
p =
0.9000
0.8000
-0.0000
z =
1.0e-03 *
0.9266
-0.4633 + 0.8017i
-0.4633 - 0.8017i
showing the same poles as before, which means a third one has been cancelled.
(a) Find the impulse response h[n] of this filter. Is this filter causal?
(b) Suppose that the input of this filter x[n] has a DTFT X(ejω ) = H(ejω ), what is the output
of the filter y[n].
(c) According to the above result, does that mean that (h ∗ h)[n] = h[n]? Explain.
Answers: h[n] is non–causal; yes, h[n] = (h ∗ h)[n].
Solution
(a) Impulse response
Z π/2
1 jωn 0.5 n=0
h[n] = e dω =
2π −π/2 6 0
sin(πn/2)/(πn) n =
1
Chaparro-Akan — Signals and Systems using MATLAB 11.2
Solution
(a) (a) (b) DTFT
X(ejω ) = 1 − e−j2ω = 2je−jω sin(ω)
Yes, X(ejω ) is periodic of period 2π since ω + 2kπ = ω. |X(ejω )| = 2| sin(ω)| is periodic
|X(ejω )|
ω
−π π
∠X(ejω )
ω
−π π
and because of the odd symmetry of the phase it is zero at ω = 0. See Fig. 11.1.
11.3 The transfer function of an FIR filter is H(z) = z −2 0.5z + 1.2 + 0.5z −1 .
(a) Find the frequency response H(ejω ) of this filter. Is the phase response of this filter linear?
(b) Find the impulse response h[n] of this filter. Is h[n] symmetric with respect to some n?
How does this relate to the phase?
Solution
(a) Frequency response
0.5 z −1 + |{z}
H(z) = |{z} 1.2 z −2 + |{z}
0.5 z −3
h[1] h[2] h[3]
11.4 A periodic discrete-time signal x[n] with a fundamental period N = 3 is passed through a
filter with impulse response h[n] = (1/3)(u[n] − u[n − 3]). Let y[n] be the filter output. We
begin the filtering at n = 0, i.e., we are interested in y[n], n ≥ 0 and could assume y[n] = 0,
n < 0. For a period, x[n] is x[0] = 1, x[1] = −2, x[2] = 1.
(a) Determine the values of y[0], y[1] and y[2]. Because x[n] is periodic you can assume its
values for n < 0 are given.
(b) Use the convolution sum to calculate y[n], n ≥ 0. What is the steady state output? When
is the steady-state of y[n] attained?
(c) Compute the output of the filter using the DTFT.
Answers: y[0] = y[1] = y[2] = 0; H(ej0 ) = 1, H(ej2π/3 ) = 0.
Solution
(a) The impulse response is h[n] = (1/3)[δ[n] + δ[n − 1] + δ[n − 2] so that y[n] = (1/3)(x[n] +
x[n − 1] + x[n − 2]) for n ≥ 0, so
y[0] = (1/3)(x[0] + x[−1] + x[−2]) = 0
y[1] = (1/3)(x[1] + x[0] + x[−1]) = 0
y[2] = (1/3)(x[2] + x[1] + x[0]) = 0
e−jω jω e−jω
H(ejω ) = (e + 1 + e−jω ) = (1 + 2 cos(ω))
3 3
Solution
which gives X[0] = 2 and X[1] = 0. Since det(F) = −2 6= 0 the matrix is invertible with
inverse 0.5FT = 0.5F and
−1 X[0] x[0]
F =
X[1] x[1]
(b) We have
FT
F=I
|det(F)|
| {z }
F−1
so α = 2.
11.6 The input of a discrete-time system is x[n] = u[n] − u[n − 4] and the impulse of the system is
h[n] = δ[n] + δ[n − 1] + δ[n − 2]
(a) Calculate the DFTs of h[n], x[n] of length N = 7. Call them H(k) and X(k). Calculate
X(k)H(k) = Ŷ (k)
(b) If y(n) = (x ∗ h)[n], would you expect ŷ[n] = y[n]? Explain. If you had chosen the lengths
of the DFTs to be N = 4 would you have had the circular and the linear convolutions
coincide? Explain.
P3 P2
Answers: (a) X[k] = n=0 e−j2πnk/7 , H[k] = n=0 e−j2πnk/7 ; (b) yes, they coincide.
Solution
(a) The DFTs of length 7 are
2
X
H[k] = e−j2πnk/7
n=0
3
X
X[k] = e−j2πnk/7
n=0
(b) Yes, the circular convolution coincides with the linear convolution if the DFTs are of
length 7, which is larger than the length of the linear convolution (4 + 3 − 1 = 6). If N < 6
the two convolutions will not coincide.
11.7 Frequency shift of FIR filters — Consider a moving average FIR filter with an impulse re-
sponse h[n] = 13 (δ[n] + δ[n − 1] + δ[n − 2]). Let H(z) be the Z-transform of h[n].
Solution
(a) The Z-transform of h[n] is
1
H(z) = (1 + z −1 + z −2 )
3
with ROC the whole Z–plane, except the origin, so we can let z = ejω to get
1 −jω jω 1
H(ejω ) = e (e + 1 + e−jω ) = e−jω (1 + 2 cos(ω))
3 3
(b)
X 1 −j(ω−π)
H1 (ejω ) = h[n]e−j(ω−π)n = H(ej(ω−π) ) = e (1 + 2 cos(ω − π))
n
3
1 −jω
= e (2 cos(ω) − 1)
3
(c) The frequency response of the two filters is computed and plotted using the following
script
% Pr. 11_7
clear all;clf
b1=(1/3)*[1 1 1];
b2=(1/3)*[1 -1 1];
[H1,w]=freqz(b1,1);
[H2,w]=freqz(b2,1);
figure(1)
subplot(221)
plot(w/pi,abs(H1));ylabel(’|H_1|’)
subplot(222)
plot(w/pi,angle(H1));ylabel(’<H_1’)
subplot(223)
plot(w/pi,abs(H2));ylabel(’|H_2|’);xlabel(’\omega/\pi’)
subplot(224)
plot(w/pi,angle(H2));ylabel(’<H_2’);xlabel(’\omega/\pi’)
0.8 1
0.6 0
|H1|
<H1
0.4 −1
0.2 −2
0 −3
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
1 3
0.8 2
0.6 1
|H2|
<H2
0.4 0
0.2 −1
0 −2
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
ω/π ω/π
Figure 11.2: Problem 7: Magnitude and phase of low–pass H1 (z) (top) and high–pass H2 (z) filters
11.8 Computations from DTFT definition — For simple signals it is possible to obtain some infor-
mation on their DTFTs without computing it. Let
(a) Find X(ej0 ) and X(ejπ ) without computing the DTFT X(ejω ).
(b) Find Z π
|X(ejω )|2 dω
−π
jω
(c) Find the phase of X(e ). Is it linear?
Use MATLAB to verify your results.
Answers: X(ej0 ) = 9, X(ejπ ) = 1, phase is linear.
Solution
(a) Using
∞
X
X(ejω ) = x[n]e−jnω
n=−∞
we have
∞
X
X(ej0 ) = x[n] = 1 + 2 + 3 + 2 + 1 = 9
n=−∞
∞
X
X(ejπ ) = (−1)n x[n] = 1 − 2 + 3 − 2 + 1 = 1
n=−∞
we get that Z π
|X(ejω )|2 dω = 2π(1 + 4 + 9 + 4 + 1) = 38π
−π
X(z) = 1 + 2z −1 + 3z −2 + 2z −3 + z −4 = z −2 [z 2 + 2z 1 + 3 + 2z −1 + z −2 ]
% Pr. 11_8
clear all;clf
x=[1 2 3 2 1];
N=256; X=fft(x,N);
w=[0:N-1]*2*pi/N;
figure(1)
subplot(211)
plot(w,abs(X));axis([min(w) max(w) -0.1 10]);grid
ylabel(’|X|’)
subplot(212)
plot(w,unwrap(angle(X)));axis([min(w) max(w) -15 1]);grid
ylabel(’<X’);xlabel(’\omega’)
10
6
|X|
0
0 1 2 3 4 5 6
−5
<X
−10
−15
0 1 2 3 4 5 6
ω
11.9 Linear equations and Fourier series — The Fourier series of a signal x[n] and its coefficients
Xk are both periodic of some value N and as such can be written
N
X −1
(i) x[n] = Xk ej2πnk/N 0≤n≤N −1
k=0
N −1
1 X
(ii) Xk = x[n]e−j2πnk/N 0≤k ≤N −1
N n=0
X = Sx
% Pr. 11_9
clear all; clf
x1=[0 1 2 0];
N=4;n=0:N-1;k=n;
S=exp(-j*2*pi.*(n’*k)/N)
Xk=(1/N)*S*x1’
x=inv(S)*N*Xk
S =1.0000 1.0000 1.0000 1.0000
1.0000 0.0000 - 1.0000i -1.0000 - 0.0000i -0.0000 + 1.0000i
where we have simplified some of the exponentials. Using that e−j2πk/2 = (−1)k and e−j2πk/4 =
(−j)k we finally obtain
Xk = 0.25 x[0] + x[2](−1)k + (−j)k x[1] + x[3](−1)k 0≤k≤3
1 T 1 T
X= S x[0] x[1] x[2] x[3] = SP x[0] x[2] x[1] x[3]
4 4
where 41 SP = S1 and P is a permutation matrix
1 0 0 0 1 1 1 1
1 0 0 1 0 = 1
1 −1 −j −j
S1 = S
4 0 1 0 0 4 1 1 −1 −1
0 0 0 1 1 −1 j −j
| {z }
P
so that
X0 x[0]
X1 x[2]
X2 = S1 x[1]
X3 x[3]
11.10 DFT and IIR filters — A definite advantage of the FFT is that it reduces considerably the
computation in the convolution sum. Thus if x[n], 0 ≤ n ≤ N − 1, is the input of an FIR
filter with impulse response h[n], 0 ≤ n ≤ M − 1, their convolution sum y[n] = (x ∗ h)[n] will
be of length M + N − 1. Now if X[k] and H[k] are the DFTs (computed by the FFT) of x[n]
and h[n], then Y [k] = X[k]H[k] is the DFT of the convolution sum, of length M + N − 1. To
multiply the FFTs X[k] and H[k] they both should be of the same length as Y [k], i.e., M +N −1.
Consider what happen when the filter is IIR which has possibly an impulse response of very
large length. Let
y[n] − 1.755y[n − 1] + 0.81y[n − 2] = x[n] + 0.5x[n]
be the difference equation representing an IIR filter with input x[n] and output y[n]. Assume
the initial conditions are zero, and the input is x[n] = u[n] − u[n − 50]. Use MATLAB to obtain
the filter output.
(a) Compute using filter the first 40 values of the impulse response h[n], and call them ĥ[n],
an approximate of h[n]. Compute the filter output ŷ[n] using the FFT as indicated above.
In this case, we are approximating the IIR filter by and FIR filter of length 40. Plot the
input and the output. Use FFTs of length 128.
(b) Suppose now that we do not want to approximate h[n], so consider the following proce-
dure. Find the transfer function of the IIR filter, say H(z) = B(z)/A(z), and if X(z) is the
Z-transform of the input then
B(z)X(z)
Y (z) =
A(z)
Compute as before the FFT for x[n], of length 128, call it X[k], and compute the 128-length
of the coefficients of B(z) and A(z) to obtain DFTs B[k] and A[k]. Multiplying X[k] by
B[k] and dividing by A[k], all of length 128, results in a sequence of length 128 that should
correspond to Y [k], the DFT of y[n]. Compute the inverse FFT to get y[n] and plot it.
(c) Use filter to solve the difference equation and obtain y[n] for x[n] = u[n] − u[n − 50].
Considering this the exact solution, calculate the error with respect to the other responses
in (a) and (b). Comment on your results.
Answers: h[n] decreases by n = 40.
Solution
(a) The impulse response of the IIR filter appears to taper down around n = 30 so the approx-
imation of the impulse response seems good.
% Pr. 11_10
clear all;clf
x=ones(1,50);
b=[1 0.5 0]; a=[1 -1.755 0.81];
delta=[1 zeros(1,30)];
h=filter(b,a,delta);
N=128;
X=fft(x,N);H=fft(h,N);
Y=X.*H; y=real(ifft(Y));
figure(1)
subplot(311)
stem([x zeros(1,10)])
subplot(312)
stem([h zeros(1,10)])
subplot(313)
stem(y)
(b)(c) In this case we do not truncate the impulse response, and the results y1 [n] are better than
before y[n], (the control is the solution obtained using filter, yc [n]). See right figure in Fig. 11.4.
The error is clearly seen when compared to the first result, but there is no visible error when
compared to the second.
% part 2
X=fft(x,N);A=fft(a,N);B=fft(b,N);
Y1=X.*B./A;
y1=ifft(Y1);
yc=filter(b,a,x);
figure(2)
subplot(211)
plot(y)
hold on
plot(yc,’r’);grid
hold off
subplot(212)
plot(y1)
hold on
plot(yc, ’r’); grid
hold off
1 40
30
x[n]
0.5
y1[n], yc[n]
20
0 10
0 10 20 30 40 50 60
0
4
3 −10
0 20 40 60 80 100 120 140
2
h[n]
1
0 40
0 5 10 15 20 25 30 35 40 30
y1[n],yc[n]
20
30
20 10
y[n]
10
0
0
−10
0 20 40 60 80 100 120 0 20 40 60 80 100 120 140
n n
Figure 11.4: Problem 10: Input x[n], impulse response ĥ[n] and output y[n] (left); right: exact re-
sponse compared to solution when h[n] is truncated (top) and when it is not (bottom)
Since the symmetry of the window is circular, express the above equation in terms of polar
coordinates, i.e., let x + jy = rejθ and ω1 + jω2 = qejφ .
(a) Using polar coordinates show that the 2D-Fourier transform becomes
Z ∞ Z 2π
W (q) = wc (r) exp (−2πjrq(cos(θ − φ))rdrdθ
0 0
(b) Lettingρ = θ − φ, show that the above equation for W (q) can be written
Z ∞
W (q) = rwc (r)J0 (2πqr)dr
0
where Z 2π
J0 (2πqr) = e−2πrq cos(ρ) dρ
0
is the zero-order Bessel function of the first kind.
(c) Replacing the circular window we finally obtain
Z 1 Z 2πq
1 J1 (2πq)
W (q) = rJ0 (2πqr)dr = r0 J0 (r0 )dr0 =
0 (2πq)2 0 2πq
where J1 (.) is the first order Bessel function of the first king.
(d) Use MATLAB to generate a circular discrete window wc [m, n] and compute its 2D-FFT
to get Wc (k, `). Determine what type of symmetry Wc (k, `) has. What does it look like?
Solution
%Pr. 11_11
clear all; clf
N=125; [m,n] = freqspace(N,’meshgrid’);
wc = ones(N);
r = sqrt(m.ˆ2 + n.ˆ2);
wc(r>0.05)=0;
figure(1)
subplot(311)
colormap(’gray’)
20 40 60 80 100 120
Magnitude response
40
20
0
150 150
100 100
50 50
0 0
20 40 60 80 100 120
Figure 11.5: Problem 11: Magnitude response of circular windows. Contour plot of circular win-
dow (top), magnitude response of window (middle) and its contour.
11.12 Separable and non-separable cosines —Consider the separable and non-separable cosine sig-
nals
(a) Use the MATLAB function fft2 to compute the FFTs of these two functions, and plot their
magnitude functions |Xi (ω1 , ω2 )|, i = 1, 2 using the function mesh (let the frequencies
ω1 , ω2 in the plots be in the domain [−π, π] × [−π, π]) and comment on their difference
(use the function contour to help you).
(b) Given a signal
1 0 ≤ m, n ≤ 99
x[m, n] =
0 100 ≤ m, n ≤ 123
Generate modulated signals y1 [m, n] = x[m, n]x1 [m, n] and y2 [m, n] = x[m, n]x2 [m, n],
and using fft2 determine the magnitude functions |Yi (ω1 , ω2 )|, i = 1, 2 and comment on
their differences.
Solution
%Pr. 11_12
% 2D-FFT of separable, non-separable and modulated cosines
clear all; clf
% (a) separable and non-separable cosines
for m=1:124,
for n=1:124,
x1(m,n)=cos(0.2*pi*(n-1))*cos(0.3*pi*(m-1)); % separable
x2(m,n)=cos(0.25*pi*(n-1)+0.5*pi*(m-1)); % non-separable
end
end
X1=fft2(x1); X2=fft2(x2);
X1a=abs(fftshift(X1)); X2a=abs(fftshift(X2));
w1=-pi:2*pi/124:pi-2*pi/124;w1=w1/pi; w2=w1;
figure(1)
subplot(221)
mesh(w1,w2,X1a);
axis([min(w1) max(w1) min(w2) max(w2) 0 max(max(X1a))])
xlabel(’\omega_1/\pi ’), ylabel(’\omega_2 /\pi’)
title(’Magnitude spectrum of x1[m,n]’)
subplot(222)
contour(w1,w2,X1a);grid
xlabel(’\omega_1/\pi’), ylabel(’\omega_2 /\pi’)
title(’Contour plot of the magnitude spectrum’)
subplot(223)
mesh(w1,w2,X2a); view(25,20);
M=max(max(X2a));
axis([min(w1) max(w1) min(w2) max(w2) 0 M])
xlabel(’\omega_1 /\pi’), ylabel(’\omega_2 /\pi’)
title(’Magnitude spectrum of x2[m,n]’)
subplot(224)
contour(w1,w2,X2a);grid
xlabel(’\omega_1 /\pi’), ylabel(’\omega_2 /\pi’)
title(’Contour plot of the magnitude spectrum’)
% (b) modulated cosines
x=zeros(124,124);
for m=1:100,
for n=1:100,
x(m,n)=1;
end
end
y1=x*x1; y2=x*x2;
Y1=fft2(y1); Y1a=abs(fftshift(Y1));
Y2=fft2(y2); Y2a=abs(fftshift(Y2));
figure(2)
subplot(221)
mesh(w1,w2,Y1a);
axis([min(w1) max(w1) min(w2) max(w2) 0 max(max(Y1a))])
xlabel(’\omega_1 /\pi’), ylabel(’\omega_2 /\pi’)
title(’Magnitude spectrum of y1[m,n]’)
subplot(222)
contour(w1,w2,Y1a);grid
xlabel(’\omega_1 /\pi’), ylabel(’\omega_2 /\pi’)
title(’Contour plot of the |Y_1(k,l)|’)
subplot(223)
mesh(w1,w2,Y2a);
axis([min(w1) max(w1) min(w2) max(w2) 0 max(max(Y2a))])
xlabel(’\omega_1 /\pi’), ylabel(’\omega_2 /\pi’)
title(’Magnitude spectrum of y2[m,n]’)
subplot(224)
contour(w1,w2,Y2a);grid
xlabel(’\omega_1 /\pi’), ylabel(’\omega_2 /\pi’)
title(’Contour plot of the |Y_2(k,l)|’)
Magnitude spectrum of x1[m,n] Contour plot of the magnitude spectrum Magnitude spectrum of y1[m,n] Contour plot of the |Y 1 (k,l)|
10 -12
2000 0.5 0.5
6
4
1000
/
0 0
2
2
0 -0.5 0 -0.5
0.5 0.5
0 0.5 0 0.5
-0.5 0 -1 -0.5 0 -1
-0.5 -0.5
-1 -1 -1 -0.5 0 0.5 -1 -1 -1 -0.5 0 0.5
2
/ / 2
/ /
1 1
1
/ 1
/
Magnitude spectrum of x2[m,n] Contour plot of the magnitude spectrum Magnitude spectrum of y2[m,n] Contour plot of the |Y 2 (k,l)|
10 -10
0.5 0.5
1.5
4000
1
3000
/
/
0 0
2
2
2000 0.5
0 0.5 0.5
0 0 0.5
-1 -0.5 -0.5 -1 -0.5 0 -1
0 0.5 -0.5
-1 -1 -0.5 0 0.5 -1 -1 -1 -0.5 0 0.5
/ / /
/ 2 2 1
1
1
/ 1
/
Figure 11.6: Problem 12: Results of parts (a) separable and non-separable and (b) modulated.
11.13 DCT vs DFT —To see the advantage of the DCT over the DFT consider a one-dimensional
edge having the values
x = [8 16 24 32 40 48 4 8]
Compute the FFT and DCT of length 8 and choose the first 4 values of the FFT and DCT
sequences and set the other 4 to zero . Compute the IFFT and IDCT of the modified FFT and
DCT sequences to obtain sequences that approximates the edge sequence. Plot the original
sequence, and the approximate sequences obtained from the IFFT and IDCT and plot the
approximation error. Find the mean square error corresponding to the two approximating
sequences and determine which gives the smallest error. Why could you say that the FFT-
based procedure has an ’unfair’ advantage over the DCT-based procedure.
Solution
% Pro 11.13
% comparison of DCt and FFT
%
x=[8 16 24 32 40 48 4 8];
X=fft(x)
Cx=dct(x)
X1=zeros(1,8); X1(1:4)=X(1:4);
x1=floor(real(ifft(X1)))
Cx1=zeros(1,8); Cx1(1:4)=Cx(1:4);
xx=floor(idct(Cx1))
figure(1)
subplot(211)
stem(x,’filled’); hold on; stem(x1,’red’,’filled’); hold on; plot(x);grid;hold on;
plot(x1,’r’); hold off; title(’DFT approximation’)
subplot(212)
stem(x,’filled’); hold on; stem(xx,’red’,’filled’);hold on; plot(x);grid; hold on
plot(xx,’r’); hold off;title(’DCT approximation’)
figure(2)
subplot(211)
stem(x-x1,’filled’);hold on; plot(x-x1);hold off;title(’DFT error’)
subplot(212)
stem(x-xx,’filled’);hold on; plot(x-xx); hold off; title(’DCT error’)
edft=sum((x-x1).*(x-x1))/8
edct=sum((x-xx).*(x-xx))/8
The FFT-based procedure is using complex values and the DCT-based procedure real values,
so the first is using twice as many values as the second one.
40 10
5
30
0
20
-5
10 -10
0 -15
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
40 10
5
30
0
20
-5
10 -10
0 -15
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
Answers: (a) Yes to both; (b) h[n] = 0.5δ[n + 2] + 1.5 + 0.5δ[n − 2].
Solution
1
Chaparro-Akan — Signals and Systems using MATLAB 12.2
which is different from zero for n = 0 and n = ±2, and h[0] = 1.5 and h[−2] = h[2] = 0.5.
Indeed, then
12.2 An FIR filter has a transfer function H(z) = z −2 (z − ejπ/2 )(z − e−jπ/2 ).
(a) Find and plot the poles and zeros of this filter.
(b) Expressing the frequency response at some frequency ω0 as
Z~1 (ω0 )Z~2 (ω0 )
H(ejω0 ) =
P~1 (ω0 )P~2 (ω0 )
carefully draw the vectors on the pole-zero plot.
(c) Consider the case when ω0 = 0, find H(ej0 ) analytically and using the vectors.
(d) Repeat the above calculations for ω0 = π and ω0 = ±π/2 (verify your result from the
given Z-transform). Can you classify this filter?
Answers: H(ej0 ) = 2; H(ejπ ) = 2; H(e±jπ/2 ) = 0.
Solution
(a) Poles z = 0 double; zeros z1,2 = e±jπ/2 = ±j.
(b) Vectors Z~1 (ω0 ) and Z~2 (ω0 ) go from the zeros to the frequency where we are finding the
frequency response. Likewise, vectors P~1 (ω0 ) and P~2 (ω0 ) go from the poles to the fre-
quency for which we are finding the frequency response.
(c) Analytically, |H(ej0 )| = (1 − j)(1 + j)/1 = 2. Graphically
√ √
j0 |Z~1 (0)||Z~2 (0)| 2 2
|H(e )| = = =2
~
(|P1 (0)| 2 1
and the phase would be
∠Z~1 (0) + ∠Z~2 (0) − 2∠P~1 (0) = (2π − φ) + φ + 0 = 2π
so H(ej0 ) = 2. Notice that when ω = 0 then H(ej0 ) = H(z)|z=1 .
2
×
Figure 12.1: Problem 2: Poles and zeros and vectors to ω = 0 to find H(ej0 ).
Solution
1 + (−2.5) z −1 + |{z}
H(z) = |{z} 1 z −2
| {z }
h[0] h[1] h[2]
so the phase is
∠H(ejω ) = −ω ± π
or linear. This corresponds to the impulse response being symmetric about n = 1 and
having no zeros on the unit circle.
12.4 Consider the following problems related to the specification of IIR filters
(a) The magnitude specifications for a low-pass filter are
1 − δ ≤ |H(ejω )| ≤ 1 0 ≤ ω ≤ 0.5π
jω
0 < |H(e )| ≤ δ 0.75π ≤ ω ≤ π
i. Find the value of δ that gives the following equivalent loss specifications for this filter
Solution
(a) i. From αmin = 20 = −20 log10 (δ) then δ = 10−1 = 0.1 and it can be verified that
αmax = 0.92 = −20 log10 (1 − δ) = −20 log10 0.9.
ii. ωp = 0.5π and ωst = 0.75π
(b) i. The dc loss α(ej0 ) = 10 dB, and subtracting it we get that
αmax = 0.1 dB
αmin = 50 dB.
ii. The maximum frequency 5 kHz≤ fs /2, so the sampling frequency is fs ≥ 10 kHz.
iii. Using ω = 2πf /fs and defining α̂(ejω ) = α(ejω ) − 10 we get
find the poles and zeros for H1 (z). What type of filter is it?
(c) Is it true, in general, that if every z in a low-pass filter transfer function is changed into
−z you obtain a high-pass filter? Explain.
Answers: H(z) = z/(z − 0.5); H1 (z) is high-pass filter.
with zero at z = 0 and pole at z = −0.5. The filter is a high-pass filter, as the frequency
response of h1 [n] = ejπn h[n] so that H1 (ejω ) = H(ej(ω−π) ) which is the low-pass filter shifted
to π.
(c) Suppose X
H(Z) = h[n]Z −n
n
so that the resulting filter has h[n](−1)n as impulse response, therefore it is high-pass.
√
2(z − 1)(z 2 + 2 z + 1)
H(z) =
(z + 0.5)(z 2 − 0.9z + 0.81)
(a) Develop a cascade realization of H(z) using a first-order and a second-order sections.
Use minimal direct form to realize each of the sections.
(b) Develop a parallel realization of H(z) by considering first and second-order sections,
each realized using minimal direct form.
Answers: Parallel H(z) = −4.94 + 2.16/(1 + 0.5z −1 ) + (4.78 − 1.6z −1 )/(1 − 0.9z −1 + 0.81z −2 )
2
x[n] + + + + y[n]
− + + +
−1
z z −1
√
0.5 −2 0.9 2
z −1
−0.81
(b) To obtain the parallel realization we do a partial fraction expansion of H(z). Notice that
H(z) is not proper rational and so we have to obtain a constant term before performing the
partial fraction expansion
−4.94
z −1
−0.5
4.78
+ +
z −1
0.9 −1.6
−1
z
−0.81
Figure 12.3: Problem 6: Parallel realization of H(z).
where α and β are constants, and x[n] is the input and y[n] is the output of the filter.
(a) Determine the transfer function H(z) = Y (z)/X(z) of the filter and from it find the fre-
quency response H(ejω ) of the filter in terms of α and β.
(b) Let α = 0.5, find then β so that the d.c. gain of the filter is unity, and the filter has a zero
phase. For the given and obtained values of α and β, sketch H(ejω ) and find the poles
and zeros of H(z) and plot them in the Z-plane.
(c) Suppose we let v[n] = y[n−1] be the output of a second filter. Is this filter causal? Find its
transfer function G(z) = V (z)/X(z). Use MATLAB to compute the unwrapped phases
of G(z) and to plot the poles and zeros of G(z) and H(z) and explain the relation between
G(z) and H(z).
Answers: H(ejω ) = β(1 + 2α cos(ω)); β = 0.5
letting z = ejω
H(ejω ) = β(αe−jω + 1 + αejω ) = β(1 + 2α cos(ω))
(b) The dc gain is
H(ej0 ) = β(1 + 2α) = 1
so that
1
β=
1 + 2α
For the phase to be zero H(ejω ) ≥ 0 and so we could let β > 0, and 1 + 2α cos(ω) ≥ 0. The last
equation requires that for all ω the amplitude of the cos(ω) be 2α ≤ 1 or α ≤ 1/2.
If we let α = 0.5 then β = 0.5 and
with double poles at z = 0 and double zero at z = −1. Since H(z) is zero phase, G(z) has
linear phase ∠G(ejω = −ω.
The following script computes and plots the frequency responses of the two filters.
% Pr. 12.7
clear all; clf
N=128; w=[0:N/2-1]*(2*pi/N);
H=0.5*(1+cos(w));
figure(1)
subplot(221)
plot(w, abs(H)); axis([0 pi 0 1]); ylabel(’|H|’)
subplot(222)
plot(w,angle(H));axis([0 pi -0.1 0.1]); ylabel(’<H’)
beta=0.5;alpha=0.5;
b=[0.25 0.5 0.25]
[G,w]=freqz(b,1);
subplot(223)
plot(w,abs(G)); axis([0 pi 0 1]); ylabel(’|G|’)
subplot(224)
plot(w,unwrap(angle(G))); axis([0 pi -pi 0]); ylabel(’<G’); xlabel(’\omega’)
1 0.1
0.8
0.05
0.6
|H|
<H
0
0.4
−0.05
0.2
0 −0.1
0 1 2 3 0 1 2 3
1 0
0.8 −0.5
−1
0.6
−1.5
<G
|G|
0.4 −2
0.2 −2.5
−3
0
0 1 2 3 0 1 2 3
ω ω
Figure 12.4: Problem 7: Frequency responses of H(e ) (top) and G(ejω ) (bottom).
jω
12.8 Effect of phase on filtering — Consider two filters with transfer functions
10
−100 1 − 2z −1
(i) H1 (z) = z , (ii) H2 (z) = 0.5
1 − 0.5z −1
(a) The magnitude response of these two filters is unity, but that they have different phases.
Find analytically the phase of H1 (ejω ) and use MATLAB to find the unwrapped phase of
H2 (ejω ) and to plot it.
(b) Consider the MATLAB signal handel.mat (a short piece of the Messiah by composer George
Handel), use the MATLAB function filter to filter it with the two given filters. Listen to
the outputs, plot them and compare them. What is the difference (look at the first 200
samples of the outputs from the two filters)?
(c) Can you recover the original signal by advancing either of the outputs? Explain.
Answers: Use MATLAB function conv to find coefficients H2 (z).
Solution: (a) Clearly H1 (z) has unit magnitude. H2 (z) is all-pass filter: let H2 (z) = (F (z))10 ,
1 − 2e−jω 1 − 2ejω 1 − 2e−jω −2ejω (1 − 0.5e−jω )
|F (ejω )|2 = 0.5 0.5 = 0.5 0.5 =1
1 − 0.5e−jω 1 − 0.5ejω 1 − 0.5e−jω −0.5ejω (1 − 2e−jω )
so |H2 (ejω )| = |F (ejω )|10 = 1. Since H2 (z) is the product of F (z) ten times, its coefficients can
be computed using the function conv to find the coefficients.
(b) The following script computes the coefficients of the filters and their frequency responses
and then processes handel.
% Pr. 12.8
clear all; clf
b1=[zeros(1,99) 1];
b=[0.5 -1];a=[1 -0.5];
b2=b;a2=a;
for k=2:10; % computing coefficients of H_2(z)
b2=conv(b,b2);
a2=conv(a,a2);
end
[H1,w]=freqz(b1,1);
[H2,w]=freqz(b2,a2);
figure(1)
subplot(221)
plot(w,abs(H1));ylabel(’|H1|’);axis([0 pi 0 1.2])
subplot(222)
plot(w,unwrap(angle(H1)));ylabel(’<H1’);axis([0 pi -320 0])
subplot(223)
plot(w,abs(H2));ylabel(’|H2|’);axis([0 pi 0 1.2])
subplot(224)
plot(w,unwrap(angle(H2)));ylabel(’<H1’);axis([0 pi -32 0])
subplot(311)
plot(x(1:200));ylabel(’x[n]’);axis([0 199 -0.4 0.4])
subplot(312)
plot(y1(1:200));ylabel(’y_1[n]’);axis([0 199 -0.4 0.4])
subplot(313)
plot(y2(1:200));ylabel(’y_2[n]’);axis([0 199 -0.4 0.4])
xlabel(’n’)
The output of H1 (z) is the input signal delayed 100 samples, the output of H2 (z) is not due to
its phase not being exactly linear, although the filter is all-pass.
(c) If the phase is linear advancing the signal 100 samples gives the original signal, while that
is not possible for the non-linear phase.
0 0.4
1 −50 0.2
x[n]
0.8 −100 0
|H1|
−150
<H1
0.6 −0.2
−200 −0.4
0.4 0 20 40 60 80 100 120 140 160 180
−250
0.2
0.4
−300
0 0.2
0 1 2 3 0 1 2 3
y1[n]
0 −0.2
1 −5 −0.4
0 20 40 60 80 100 120 140 160 180
0.8 −10
0.4
−15
|H2|
<H1
0.6 0.2
−20
y2[n]
0.4 0
−25
0.2 −0.2
−30
0 −0.4
0 1 2 3 0 1 2 3 0 20 40 60 80 100 120 140 160 180
n
Figure 12.5: Problem 8: Magnitude and phase of filters (left). Part of handel signal (top right) and outputs
of linear–phase and non–linear phase filters.
12.9 Bilinear transformation and pole location — Find the poles of the discrete filter obtained by
applying the bilinear transformation with K = 1 to and analog second–order Butterworth
low-pass filter. Determine the half-power frequency ωhp of the resulting discrete filter. Use
the MATLAB function bilinear to verify your results. q
√ √
Answers: Double zero at z = −1 and poles at z1,2 = ±j (2 − 2)/(2 + 2).
1
H(s) = √
s2 + 2s + 1
(1 + z −1 )2
H(z) = √
(1 − z −1 )2 + 2(1 − z −2 ) + (1 + z −1 )2
(1 + z −1 )2 (z + 1)2
= √ √ = √ √ √
(2 + 2) + (2 − 2)z −2 (2 + 2)(z 2 + (2 − 2)/(2 + 2))
q √ √
so we have a double zero at z = −1 and poles at z1,2 = ±j (2 − 2)/(2 + 2).
We have that Kb transforms the normalized half-power analog frequency Ωhp = 1 into the
normalized half-power frequency ωhp of the discrete filter, so that
1
Kb = 1 = ⇒ ωhp = π/2
tan(ωhp /2)
% Pr. 12.9
clear all;clf
n=1;d=[1 sqrt(2) 1];
[b,a]=bilinear(n,d,0.5)
[H,w]=freqz(b,a)
figure(1)
subplot(211)
zplane(b,a)
subplot(212)
plot(w/pi,abs(H));grid;ylabel(’|H(eˆ{j\omega})|’); xlabel(’\omega/\pi’)
0.5
Imaginary Part
0 2
−0.5
−1
−3 −2 −1 0 1 2 3
Real Part
0.8
|H(ejω)| 0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ω/π
12.10 Discrete Butterworth filter for analog processing — Design a Butterworth low-pass discrete
filter that satisfies the following specifications
0 ≤ α(ejω ) ≤ 3 dB for 0 ≤ f ≤ 25 Hz
jω
α(e ) ≥ 38 db for 50 ≤ f ≤ F s/2 Hz
and the sampling frequency is Fs = 2000 Hz. Express the transfer function H(z) of the de-
signed filter as a cascade of filters. Use first the design formulas and then use MATLAB to
confirm your results. Show that the designed filter satisfies the specifications, plotting the loss
function of the designed filter.
Q4
Answers: N = 7; H(z) = G i=1 Hi (z); H1 (z) = 0.04(1 + z −1 )/(1 − 0.92z −1 ).
% Pr. 12.10
N = 7; fhp = 25; Fs = 2000;
whp = 2*pi*fhp/Fs;
[b, a] = butter(N, whp/pi);
% analytic results
% coefficients of H1
b1 = [0.04, 0.04]; a1 = [1, -0.92];
% coefficients of H2
b2 = 1.44e-3*[1, 2, 1]; a2 = [1, -1.86, 0.87];
% coefficients of H3
b3 = 1.47e-3*[1, 2, 1]; a3 = [1, -1.9, 0.91];
% coefficients of H4
b4 = 1.51e-3*[1, 2, 1]; a4 = [1, -1.96, 0.91];
% comparison of analytic and MATLAB results
analyticalB = conv(conv(conv(b1, b2), b3), b4)
matlabB = b
analyticalA = conv(conv(conv(a1, a2), a3), a4)
matlabA = a
analyticalB =1.0e-08 *
0.0128 0.0895 0.2685 0.4475 0.4475 0.2685 0.0895 0.0128
matlabB =1.0e-08 *
0.0121 0.0848 0.2544 0.4240 0.4240 0.2544 0.0848 0.0121
12.11 Butterworth filtering of analog signal — We wish to design a discrete Butterworth filter that can be
used in filtering a continuous-time signal. The frequency components of interest in this signal are be-
tween 0 and 1 kHz, so we would like the filter to have a maximum passband attenuation of 3 dB within
that band. The undesirable components of the input signal occur beyond 2 kHz, and we like to attenuate
them by at least 10 dB. The maximum frequency present in the input signal is 5 kHz. Finally, we would
like the d.c. gain of the filter to be 10. Choose the Nyquist sampling frequency to process the input
signal. Use MATLAB to design the filter. Give the transfer function of the filter, plot its poles and zeros,
and its magnitude and unwrapped phase response using an analog frequency scale in kHz.
Answers: Low-pass, with fp = 1 kHz (fp = fhp ), αmax = 3 dB , fst = 2 kHz, αmin = 10 dB.
Solution: The desired filter is a low-pass, with fp = 1 kHz, αmax = 3 dB, so that fp = fhp ; fst = 2
kHz, αmin = 10 dB, and dc gain of 10 after the design with dc of 0 dB. Since fmax = 5 kHz, then we can
choose Fs = 10 kHz. The coefficients of the resulting filter are shown below.
% Pr. 12.11
Fs=10000;fp=1000;fst=2000;
wp=2*pi*fp/Fs/pi; ws=2*pi*fst/Fs/pi;
alphamax=3; alphamin=10;
[N,wn]=buttord(wp,ws,alphamax,alphamin)
[b,a]=butter(N,wp); % to guarantee that wp is half-power frequency (wn\ne wp)
G=10*sum(b)/sum(a);
b=G*b;
[H,w]=freqz(b,a);
figure(1)
subplot(221)
zplane(b,a)
subplot(222)
plot(w/pi*Fs/2, abs(H));grid; axis([0 Fs/2 -0.1 11])
subplot(223)
plot(w/pi*Fs/2,unwrap(angle(H)));grid; axis([0 Fs/2 -4 0])
1 10
0.5 8
Imaginary Part
2 6
0
4
−0.5
2
−1 0
−1 −0.5 0 0.5 1 0 1000 2000 3000 4000 5000
Real Part
−1
−2
−3
−4
0 1000 2000 3000 4000 5000
Figure 12.7: Problem 11: Poles/zeros and magnitude and phase responses.
• Using the fft function compute the DFT of the original signal, the noisy signal and the
noise, and plot their magnitudes. Is the cutoff frequency of the filters adequate to get rid
of the noise? Explain.
• Compute and plot the magnitude and the unwrapped phase, as well as the poles and ze-
ros for each of the three filters. Comment on the differences in the magnitude responses.
• Use the function filter to obtain the output of each of the filters, and plot the original
noiseless signal and the filtered signals. Compare them.
Solution: The following script designs the desired filters, computes and plots the output of
the filtered signals and different spectra:
% Pr. 12.12
clear all; clf
load train;y=y’; N=256;
noise=0.1*randn(1,N);
yc=y(1:N);
yn=yc+noise; % noisy signal
%[b,a]=butter(20,0.5); % butterworth
[b,a]=cheby2(20,60,0.5); % cheby2
%[b,a]=cheby1(20,0.01,0.5) % cheby1
figure(1)
[H,w]=freqz(b,a);
b=b/H(1);
[H,w]=freqz(b,a); % unit dc gain
subplot(311)
plot(w/pi, abs(H)); grid
hold on
plot(w/pi,0.707*ones(1,length(w)),’r’)
hold off
subplot(312)
plot(w/pi,unwrap(angle(H)))
subplot(313)
zplane(b,a)
k=0:N-1;w1=k*2*pi/N;w1=w1/pi-1;
Ns=fftshift(abs(fft(noise))); % noise spectrum
Yc=fftshift(abs(fft(yc))); % clean signal spectrum
Yn=fftshift(abs(fft(yn))); % noisy signal spectrum
figure(2)
subplot(311)
plot(w1,Yc); title(’clean signal spectrum’)
subplot(312)
plot(w1,Yn); title(’noisy signal spectrum’)
subplot(313)
plot(w1,Ns); title(’noise spectrum’)
yy=filter(b,a,yn);
n=0:N-1;
figure(3)
subplot(311)
plot(n,yn,’r’);title(’noisy/clean signal’)
hold on
plot(n,yc);legend(’noisy’,’clean’)
hold off
subplot(312)
plot(n,yy,’r’);legend(’clean signal’)
subplot(313)
plot(n,yc); legend(’filtered’,’clean’)
hold off
1 10
0.5 5
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
0 10
−10 5
−20 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
noise spectrum
1 4
Imaginary Part
0.5 3
0 2
−0.5 1
−1 0
−5 −4 −3 −2 −1 0 1 2 3 4 5 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
Real Part noisy/clean signal
0.4
0.2 noisy
0 clean
−0.2
−0.4
−0.6
0.4
0.2 clean signal
0
−0.2
−0.4
0.4
filtered
0.2
0
−0.2
−0.4
0 50 100 150 200 250
Figure 12.8: Problem 12: Chebyshev 2 case: Filter frequency response and poles zeros (top left), spectra of
different signals, signals (bottom)
12.13 Notch and all-pass filters — Notch filters is a family of filters that include the all-pass filter.
For the filter
(1 − α1 z −1 )(1 + α2 z −1 )
H(z) = K
(1 − 0.5z −1 )(1 + 0.5z −1 )
(a) Determine the values of α1 , α2 and K that would make H(z) an all-pass filter of unit
magnitude. Use MATLAB to compute and plot the magnitude response of H(z) using
the obtained values for α and K. Plot the poles and zeros of this filter.
(b) If we would like the filter H(z) to be a notch filter, of unit gain at at ω = π/2 (rad), and
notches at ω = 0 and π, determine the values of α and K to achieve this. Use MATLAB
functions to verify that the filter is a notch filter, and to plot the poles and zeros.
(c) Use MATLAB to show that when α1 = α2 = α and 1 ≤ α ≤ 2, the given H(z) correspond
to a family of notch filters with different attenuations. Determine K so that these filters
are unity gain at ω = π/2.
(d) Suppose we use the transformation z −1 = jZ −1 to obtain a filter H(Z) using H(z) ob-
tained in the previous item. Where are the notches of this new filter? What is the differ-
ence between the all-filters H(z) and H(Z)?
Answers: H(z) = K(1−α2 z −2 )/(1−0.25z −2 ), 1 ≤ α ≤ 2, notch filters with notches at ω = 0, π.
Solution: (a) For the filter to be all-pass the zeros must be the inverse of the poles 0.5 and
−0.5 or 2 and −2. Thus
(1 − 2z −1 )(1 + 2z −1 ) (1 − 4z −2 )
H(z) = K = K K>0
(1 − 0.5z −1 )(1 + 0.5z −1 ) (1 − 0.25z −2 )
As an all-pass H(z) has a constant magnitude (at any frequency) and in this case we want it
to be unity. For instance, for this filter to have unit gain at ω = π/2 we need that
1 − 4j 2
|H(ejπ/2 )| = K = K 5 = 1 ⇒ K = 1.25 = 0.25
1 − 0.25j 2 1.25 5
(b) The zeros of H(z) are α1 and −α2 which are real. For H(z) to be a notch filter |αi | = 1, and
so we let α1 = α2 = 1 for notches at ω = 0, π. Thus
(1 − z −1 )(1 + z −1 ) (1 − z −2 )
H(z) = K −1 −1
=K
(1 − 0.5z )(1 + 0.5z ) (1 − 0.25z −2 )
and for it to have unit gain at ω = π/2
1 − j2
|H(ejπ/2 )| = K = K 2 = 1 ⇒ K = 1.25 = 0.25
1 − 0.25j 2 1.25 2
(c) In general, the filter
(1 − α2 z −2 )
H(z) = K
(1 − 0.25z −2 )
for 1 ≤ α ≤ 2 is a family of notch filters with notches at ω = 0, π, including an all-pass filter
when α = 2. The gain K > 0 is
1 − α2 j 2 1 + α2
|H(e jπ/2
)| = K
=K = 1 ⇒ K = 1.25
1 − 0.25j 2 1.25 1 + α2
(1 − α2 j 2 Z −2 ) (1 + α2 Z −2 )
H(Z) = K = K
(1 − 0.25j 2 Z −2 ) (1 + 0.25Z −2 )
with poles at ±0.5ejπ/2 and zeros at ±αejπ/2 , i.e., shifted 90o with respect to the ones before.
For values of 1 ≤ α ≤ 2, H(Z) is a family of notch filters at ω = π/2, including an all-pass filter
when α = 2. When running the following script make the numerator and the denominator
correspond to H(z) or to H(Z) by deleting the comment symbol.
% Pr. 12.13
clear all;clf
%a=[1 0 -0.25]; %% denominator of H(z)
a=[1 0 0.25]; %% denominator of H(Z)
figure(1)
for alpha=2:-0.1:1;
alpha
K=1.25/(1+alphaˆ2)
%b=[K 0 -(alphaˆ2)*K]; %% numerator of H(z)
b=[K 0 (alphaˆ2)*K]; %% numerator of H(Z)
[H,w]=freqz(b,a);
subplot(211)
zplane(b,a); axis([-2 2 -2.5 2.5])
hold on
subplot(212)
plot(w/pi,abs(H));axis([0 1 0 1.5])
hold on
pause(1)
end
hold off
Imaginary Part
1
−1
−2
−2 0 2
Real Part
1.5
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
2
Imaginary Part
1
−1
−2
−2 0 2
Real Part
1.5
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 12.9: Problem 13: Loci of poles and zeros and magnitude responses of H(z) (top) and of H(Z)
(bottom).
12.14 Modulation property transformation for IIR filters — The modulation-based frequency trans-
formation of the DTFT is applicable to IIR filters. It is obvious in the case of FIR filters, but
requires a few more steps in the case of IIR filters. In fact, if we have that the transfer function
of the prototype IIR low-pass filter is H(z) = B(z)/A(z), with impulse response h[n], let the
transformed filter be Ĥ(z) = Z(2h[n] cos(ωo n) for some frequency ω0
Solution: (a) Modulating the impulse response h[n] of an IIR filter gives a new filter with
transfer function
∞
X
Ĥ(z) = 2h[n] cos(ω0 n)z −n
n=0
X∞ ∞
X
= h[n](ejω0 z −1 )n + h[n](e−jω0 z −1 )n
n=0 n=0
= H(e−jω0 z) + H(ejω0 z)
B(ejω0 z) B(e−jω0 z)
Ĥ(z) = + .
A(ejω0 z) A(e−jω0 z)
If a = 0.5 and ω0 = π/2, i.e., we want to convert the prototype low-pass filter into a band-pass
filter, we get
2 2z 2
Ĥ(z) = −2
= 2
1 + 0.25z z + 0.25
with two zeros at z = 0 and two poles at z = ±j0.5, so that the filter is a bandpass filter with
center frequency π/2.
(c) To get a high–pass filter we let ω0 = π. The following script computes and plots the
different frequency responses.
% Pr. 12.14
clear all; clf
b=1;a=[1 -0.5];
[H,w]=freqz(b,a);
alpha=0.5; w0=pi/2
b1=[2 -2*alpha*cos(w0)];a1=[1 -2*alpha*cos(w0) alphaˆ2];
w0=pi
b2=[2 -2*alpha*cos(w0)];a2=[1 -2*alpha*cos(w0) alphaˆ2];
[H1,w]=freqz(b1,a1);
[H2,w]=freqz(b2,a2);
figure(1)
lot(311)
plot(w/pi,abs(H));axis([0 1 0 2]);grid
subplot(312)
plot(w/pi,abs(H1));axis([0 1 0 3]);grid
subplot(313)
plot(w/pi,abs(H2));axis([0 1 0 4]);grid
2
1.5
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
4
3
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 12.10: Problem 14: Magnitude responses of low–pass prototype (top) and band–pass and high–pass
filters (middle and bottom).
12.15 Effect of size of filter in denoising —The size of an FIR averaging filter has significant effect
on denoising. Consider as the clean image the peppers.png, use the function rgb2gray and double
to get a double precision gray-level image x[m, n]. Suppose it is distorted by additive Gaussian
noise with zero mean and 0.1 variance (use the function imnoise). Consider 3 × 3 and 9 × 9
low-pass averaging filters (normalized so that the maximum value of the impulse response
is unity) to filter the noisy image. Call y3 [m, n] and y5 [m, n] of the 3 × 3 and 9 × 9 filters.
To determine which of the two filters is more efficient in denoising the image, determine the
mean-square errors X
i = (yi [m, n] − x[m, n])2 /Mi ,
m,n
2 2
for i = 3, 9, M3 = 3 and M9 = 9 . Indicate which of the two filters, under this criterion, is
more effective in denoising the image.
Solution
% Pr. 12.15
% Effect of size of filter in denoising
% 3x3 vs 9x9 averaging filtering of
% signal with added Gaussian noise
Figure 15 shows the results of filtering with a 3 × 3 and a 9 × 9 averaging FIR filters. The mean-
square errors are ε3 = 0.4890 and ε9 = 0.4337 a slightly better result for the 9 × 9 averaging
filter.
CLEAN NOISY
Generate an image of size 100 × 100 that is a block of ones of size 50 × 50 surrounded by zeros,
call it x[m, n].
(a) Filter x[m, n] and display it and the Laplacian filtered image y[m, n]. Use then Sobel’s
horizontal and vertical filters to filter x[m, n] and obtain a second image z[m, n] by adding
their outputs. Which of these two provide a better edge detection?
(b) Find the magnitude response of the Laplacian filter associated with HL .
(c) Apply the Laplacian filter mask HL to the peppers image and show the original and
filtered images.
Solution
% Pr. 12.16
% Laplacian Filtering
%
clc; clear all; close all;
x=zeros(100,100);
x(26:75,26:75)=ones(50,50);
figure(1)
colormap(’gray’)
subplot(221)
imagesc(x); title(’Original Image’)
% Laplacian
HL=[0 1 0; 1 -4 1; 0 1 0];
% Sobel
hx=[1 0 -1;2 0 -2;1 0 -1];
hy=hx’;
yx=imfilter(x,hx);
yy=imfilter(x,hy);
yS=sqrt(yx.ˆ2 + yy.ˆ2);
subplot(222)
imshow(yS);title(’Sobel Filtered Image’)
yHL=imfilter(x,HL);
subplot(223)
imshow(yHL);title(’Laplacian Filtered Image’)
HLf=fftshift(abs(fft2(HL,256,256)));
w=-1+1/128:1/128:1;
subplot(224)
mesh(w,w,HLf); title(’Frequency response of the Laplacian Filter’)
xlabel(’\omega_1/\pi’)
ylabel(’\omega_2/\pi’)
I=imread(’peppers.png’);
I2=im2double(I);
I3=rgb2gray(I2);
figure(2)
colormap(’gray’)
subplot(211)
imshow(I3);title(’Peppers Image’)
yL=imfilter(I3,HL);
subplot(212)
yL=yL>0.05;
imshow(yL);
title(’Laplacian Filtered Peppers Image’)
Peppers Image
Figure 12.12: Problem 16: Laplacian vs Sobel filtering (left), and edge detection using Laplacian
(right).
12.17 Circular filter Design a 2D FIR filter, using the function fwind1, that approximates the desired
magnitude response
jω1 jω2 1 r < 0.5
Hd (e , e ) =
0 otherwise
p
where r = f12 + f22 < 0.5 for normalized frequencies −1 ≤ f1 ≤ 1 and −1 ≤ f2 ≤ 1 and
fi = ωi /π, i = 1, 2. Use a Hamming window. Determine the magnitude response of the filter,
and plot the desired frequency magnitude response and that of the designed filter using the
function contour. Normalized the frequencies to be (−1, 1) × (−1, 1).
Solution
% Pr. 12.17
% Circular filter design using fwind1
clear all; clf
% filter specifications
N=125;
[f1,f2] = freqspace(N,’meshgrid’);
Hd = ones(N);
r = sqrt(f1.ˆ2 + f2.ˆ2);
% type of filter
Hd(r>0.5)=0; %low-pass
% Hamming window
h1 = fwind1(Hd,hamming(N));
H=freqz2(h1,124); w1=-1:2/124:1-2/124;w2=w1
figure(1)
colormap(’gray’)
subplot(221)
mesh(f1,f2,Hd); title(’Magnitude response of desired filter’)
subplot(222)
contour(f1,f2,Hd);grid
subplot(223)
mesh(f1,f2,Hd); title(’Magnitude response of designed filter’)
subplot(224)
contour(w1,w2,H); grid;
Figure 12.13: Problem 17:Magnitude responses of desired (top) and of designed filters (bottom).
Circular plots are shown of each of these.