EE-384: Digital Signal Processing Spring 2023

Laboratory 4: Discrete-time Fourier Analysis

Instructor: Mr Ammar Naseer EE UET New Campus

Aims
We can represent any arbitrary discrete signal as a linear combination of basis signals. Each basis signal
set provides a new signal representation. Each representation has some advantages and some
disadvantages depending upon the type of system under consideration. However, when the system is
linear and time-invariant, only one representation stands out as the most useful. It is based on the
complex exponential signal set and is called the discrete-time Fourier transform. This lab session
discusses discrete-time signal and system representation in the frequency domain.

Pre-Lab:
The Discrete-time Fourier Transform (DTFT)
If x(n) is absolutely summable, that is∑∞
−∞ 𝑥(𝑛) < ∞ then its discrete-time Fourier transform is given by

∞
𝑗𝜔
𝑋(𝑒 ) = ℱ[𝑥(𝑛)] = ∑ 𝑥(𝑛)𝑒 𝑗𝜔𝑛 (5.1)
𝑛=−∞

The inverse discrete-time Fourier transform (IDTFT) of 𝑋(𝑒 𝑗𝜔 ) is given by

1 𝜋
𝑥(𝑛) = ℱ −1 [𝑋(𝑒 𝑗𝜔 )] = ∫ 𝑋(𝑒 𝑗𝜔 )𝑒 𝑗𝜔𝑛 𝑑𝜔 (5.2)
2𝜋 −𝜋

The operator ℱ[ ] transforms a discrete signal x(n) into a complex-valued continuous function 𝑋(𝑒 𝑗𝜔 ) of
real variable 𝜔, called a digital frequency, which is measured in radians/sample.

Since 𝑋(𝑒 𝑗𝜔 ) is a complex-valued function, we will have to plot its magnitude and its angle (or the real
and the imaginary part) with respect to separately to visually describe 𝑋(𝑒 𝑗𝜔 ). Now 𝜔 is a real variable
between -∞ and ∞, which would mean that we can plot only a part of the 𝑋(𝑒 𝑗𝜔 ) function using
MATLAB. Using two important properties of the discrete-time Fourier transform, we can reduce this
domain to the [0, π] interval for real-valued sequences.

4A. Evaluate 𝑋(𝑒 𝑗𝜔 ) of x(n) = (0.5)nu(n) at 501 equispaced points between [0, π] and plot its magnitude,
angle, real, and imaginary parts.

Page 1
Laboratory 4: The Discrete-time Fourier Analysis 4.2

>> w = [0:1:500]*pi/500; %[0, pi] axis divided into 501 points.

>> X = exp(j*w) ./ (exp(j*w) - 0.5*ones(1,501));
>> magX = abs(X);
>> angX = angle(X);
>> realX = real(X);
>> imagX = imag(X);
>> subplot(2,2,1); plot(w/pi,magX); grid; % Plot magnitude of X
>> xlabel('frequency in pi units'); title('Magnitude Part');
ylabel('Magnitude');
>> subplot(2,2,3); plot(w/pi,angX); grid; % Plot angle of X
>> xlabel('frequency in pi units'); title('Angle Part');
ylabel('Radians')
>> subplot(2,2,4); plot(w/pi,imagX); grid; % Plot imaginary of X
>> xlabel('frequency in pi units'); title('Imaginary Part');
ylabel('Imaginary')
>> subplot(2,2,2); plot(w/pi,realX); grid; % Plot real of X
>> xlabel('frequency in pi units'); title('Real Part'); ylabel('Real')

Note that we divided the w array by pi before plotting so that the frequency axes are in the units of π
and therefore easier to read. This practice is strongly recommended.

If 𝑥(𝑛) is of finite duration, then MATLAB can be used to compute 𝑋(𝑒 𝑗𝜔 ) numerically at any frequency
𝜔. If, in addition, we evaluate 𝑋(𝑒 𝑗𝜔 ) at equispaced frequencies between [0, π], then (5.1) can be
implemented as a matrix-vector multiplication operation. To understand this, let us assume that the
sequence 𝑥(𝑛) has N samples between n1 < n < nN (i.e., not necessarily between [0, N1]) and that we
want to evaluate 𝑋(𝑒 𝑗𝜔 ) at
𝜋
𝜔𝑘 = 𝑘 𝑘 = 0, 1 … . . 𝑀
𝑀

which are (M + 1) equispaced frequencies between [0,𝜋]. Then (6.1) can be written as
𝑁
𝜋
𝑗𝜔𝑘 −𝑗( )𝑘𝑛𝑖
𝑋(𝑒 ) = ∑𝑒 𝑀 𝑥(𝑛𝑖 ) 𝑘 = 0, 1 … . . 𝑀
𝑙=1

When {𝑥(𝑛𝑖 )} and {𝑋(𝑒 𝑗𝜔𝑘 )} are arranged as column vectors x and X, respectively, we have

𝑿 = 𝑾𝑥 (5.3)

where W is an (M + 1) × N matrix given by

𝜋⁄ )𝑘𝑛
𝑾 = {𝑒 −𝑗( 𝑀 1; 𝑛
1 ≤ 𝑛 ≤ 𝑛𝑁 , 𝑘 = 0, 1, … , 𝑀}

In addition, if we arrange {𝑘} and {𝑛𝑖 } as row vectors k and n respectively, then
𝜋 𝑇
𝑾 = [exp (−𝑗 𝒌 𝒏)]
𝑀
Laboratory 4: The Discrete-time Fourier Analysis 4.3

In MATLAB we represent sequences and indices as row vectors; therefore taking the transpose of (5.3),
we obtain
𝜋 𝑇
𝑿𝑇 = 𝒙𝑇 [exp (−𝑗 𝒌 𝒏] (5.4)
𝑀
Note that nTk is an N × (M + 1) matrix. Now (6.4) can be implemented in MATLAB as follows.

>> k = [0:M]; n = [n1:n2];

>> X = x * (exp(-j*pi/M)) .^ (n'*k);

4B. Numerically compute the discrete-time Fourier transform of the sequence x(n) ={1, 2, 3, 4, 5} at 501
equispaced frequencies between [0, π].

>> n = -1:3;
>> x = 1:5;
>> k = 0:500;
>> w = (pi/500)*k;
>> X = x * (exp(-j*pi/500)) .^ (n'*k);
>> magX = abs(X);
>> angX = angle(X);
>> realX = real(X);
>> imagX = imag(X);
>> % use subplot to plot the results

𝑗𝜋
4C. Let (𝑛) = (0.9 exp ( 3 )) 𝑛 , 0 ≤ n ≤ 10. Determine 𝑋(𝑒 𝑗𝜔 ) and investigate its periodicity.

Since x(n) is complex-valued, 𝑋(𝑒 𝑗𝜔 ) satisfies only the periodicity property. Therefore it is uniquely
defined over one period of 2π. However, we will evaluate and plot it at 401 frequencies over two
periods between [2π; 2π] to observe its periodicity.

>> n = 0:10;
>> x = (0.9*exp(j*pi/3)) .^n;
>> k = -200:200;
>> w = (pi/100)*k;
>> X = x * (exp(-j*pi/100)) .^ (n'*k);
>> magX = abs(X);
>> angX =angle(X);
>> % use subplot to plot the results

It can be observed that 𝑋(𝑒 𝑗𝜔 ) is periodic in 𝜔 but is not conjugate-symmetric.

4D. Let 𝑥(𝑛) = 0.9𝑛 , -10 ≤ n ≤ 10. Investigate the conjugate-symmetry property of its discrete-time
Fourier transform.

Once again compute and plot 𝑋(𝑒 𝑗𝜔 ) over two periods to study its symmetry property. We observe that
𝑋(𝑒 𝑗𝜔 ) is not only periodic in 𝜔 but is also conjugate-symmetric. Therefore for real sequences we will
plot their Fourier transform magnitude and angle graphs from 0 to π.
Laboratory 4: The Discrete-time Fourier Analysis 4.4

In this example we will verify the linearity property using real-valued finite-duration sequences. Let x1(n)
and x2(n) be two random sequences uniformly distributed between [0; 1] over 0 ≤ n ≤ 10. Then we can
use our numerical discrete-time Fourier transform procedure as follows.

>> x1 = rand(1,11);
>> x2 = rand(1,11);
>> n = 0:10;
>> alpha = 2;
>> beta = 3;
>> k = 0:500;
>> w = (pi/500)*k;
>> X1 = x1 * (exp(-j*pi/500)).^(n'*k);
>> X2 = x2 * (exp(-j*pi/500)).^(n'*k);
>> x = alpha*x1 + beta*x2;
>> X = x * (exp(-j*pi/500)).^(n'*k);
>> X_check = alpha*X1 + beta*X2;
>> error = max(abs(X-X_check))

Since the maximum absolute error between the two Fourier transform arrays is less than 10 -14, the two
arrays are identical within the limited numerical precision of MATLAB.

Frequency Domain Representation

We earlier stated that the Fourier transform representation is the most useful signal representation for
LTI systems. It is due to the following result.

Let 𝑥(𝑛) = 𝑒 𝑗𝜔0 𝑛 be the input to an LTI system represented by the impulse response ℎ(𝑛). Then
∞
𝑗𝜔0 𝑛
𝑦(𝑛) = 𝑒 ∗ ℎ(𝑛) = ∑ ℎ(𝑘)𝑒 𝑗𝜔0 (𝑛−𝑘)
−∞

= [∑ ℎ(𝑘)𝑒 −𝑗𝜔0 𝑘 ] 𝑒 𝑗𝜔0 𝑛

−∞

= [ℱ[ℎ(𝑛)]|𝜔=𝜔0 ] 𝑒 𝑗𝜔0 𝑛 (5.5)

The discrete-time Fourier transform of an impulse response is called the frequency response (or transfer
function) of an LTI system and is denoted by
∞
𝑗𝜔𝑛
𝐻(𝑒 ) = ∑ ℎ(𝑛)𝑒 𝑗𝜔𝑛 (5.6)
−∞

Hence the output sequence is the input exponential sequence modified by the response of the system
at frequency 𝜔0 . This justifies the definition of 𝐻𝑒 𝑗𝜔0 𝑛 as a frequency response because it is what the
complex exponential is multiplied by to obtain the output 𝑦(𝑛).
Laboratory 4: The Discrete-time Fourier Analysis 4.5

In general, the frequency response 𝐻𝑒 𝑗𝜔0 𝑛 is a complex function of 𝜔. The magnitude |𝐻𝑒 𝑗𝜔0 𝑛 | of
𝐻𝑒 𝑗𝜔0 𝑛 is called the magnitude (or gain) response function, and the angle < 𝐻𝑒 𝑗𝜔0 𝑛 is called the phase
response function as we shall see below.

From 5.5, we can represent the system by

𝑥(𝑛) = 𝑒 𝑗𝜔0 𝑛 𝐻(𝑒 𝑗𝜔 ) 𝑦(𝑛) = 𝐻(𝑒 𝑗𝜔0 ) × 𝑒 𝑗𝜔0 𝑛 (5.7)

Equation 5.7 can be generalized to arbitrary absolutely summable sequences. Let 𝑋(𝑒 𝑗𝜔 ) =
ℱ[𝑥(𝑛)] and 𝑌(𝑒 𝑗𝜔 ) = ℱ[𝑦(𝑛)]; then using the convolution property, we have

𝑌(𝑒 𝑗𝜔 ) = 𝐻(𝑒 𝑗𝜔 )𝑋(𝑒 𝑗𝜔 ) (5.8)

4E. Determine the frequency response 𝐻(𝑒 𝑗𝜔 ) of a system characterized by ℎ(𝑛) = (0.9)𝑛 𝑢(𝑛). Plot
the magnitude and the phase responses.

Using (5.6)
∞ ∞
𝐻(𝑒 𝑗𝜔 ) = ∑ ℎ(𝑛)𝑒 −𝑗𝜔𝑛 = ∑ (0.9)𝑛 𝑒 −𝑗𝜔𝑛
−∞ 0

∞ 1
= ∑ (0.9𝑒 −𝑗𝜔 )𝑛 =
0 1 − 0.9𝑒 −𝑗𝜔

hence

1
|𝐻(𝑒 𝑗𝜔 )| = √
(1 − 0.9 𝑐𝑜𝑠 𝜔)2 + (0.9 𝑠𝑖𝑛 𝜔)2

and

0.9 𝑠𝑖𝑛 𝜔
< 𝐻(𝑒 𝑗𝜔 ) = −arctan [ ]
1 − 0.9 𝑐𝑜𝑠 𝜔

To plot these responses, we can either implement the |𝐻(𝑒 𝑗𝜔 )| and < 𝐻(𝑒 𝑗𝜔 ) functions or the
frequency response 𝐻(𝑒 𝑗𝜔 ) and then compute its magnitude and phase. The latter approach is more
useful from a practical viewpoint.

>> w = [0:1:500]*pi/500; % [0, pi] axis divided into 501 points.

>> H = exp(j*w) ./ (exp(j*w) - 0.9*ones(1,501));
>> magH = abs(H); angH = angle(H);
>> subplot(2,1,1); plot(w/pi,magH); grid;
>> xlabel('frequency in pi units'); ylabel(' |H| ');
>> title('Magnitude Response');
>> subplot(2,1,2); plot(w/pi,angH/pi); grid
>> xlabel('frequency in pi units'); ylabel('Phase in pi Radians');
>> title('Phase Response');

When an LTI system is represented by the difference equation

Laboratory 4: The Discrete-time Fourier Analysis 4.6

𝑁 𝑀

𝑦(𝑛) + ∑ 𝑎1 𝑦(𝑛 − 𝑙) = ∑ 𝑏𝑚 𝑥(𝑛 − 𝑚) (5.9)

𝑖=1 𝑚=0
Then
𝑁 𝑀
𝑗𝜔 𝑗𝜔𝑛 𝑗𝜔 𝑗𝜔(𝑛−𝑙)
𝐻(𝑒 )𝑒 + ∑ 𝑎𝑙 𝐻(𝑒 )𝑒 = ∑ 𝑏𝑚 𝑒 𝑗𝜔(𝑛−𝑚)
𝑖=1 𝑚=0

∑𝑀
𝑚=0 𝑏𝑚 𝑒
−𝑗𝜔𝑚
𝐻(𝑒 𝑗𝜔 ) = (5.10)
1 + ∑𝑁
𝑖=0 𝑎𝑖 𝑒
−𝑗𝜔𝑙

This equation can easily be implemented in MATLAB, given the difference equation parameters.

4F. An LTI system is specified by the difference equation

𝑦(𝑛) = 0.8𝑦(𝑛 − 1) + 𝑥(𝑛)

We calculate and plot the steady-state response 𝑦𝑠𝑠 (𝑛) to

𝑥(𝑛) = cos (0.5𝜋𝑛)𝑢(𝑛)

>> b = 1; a = [1,-0.8];
>> n=[0:100];x = cos(0.05*pi*n);
>> y = filter(b,a,x);
>> subplot(2,1,1); stem(n,x);
>> xlabel('n'); ylabel('x(n)'); title('Input sequence')
>> subplot(2,1,2); stem(n,y);
>> xlabel('n'); ylabel('y(n)'); title('Output sequence')

In practice the difference equations are of large order and hence we need a compact procedure to
implement the general expression (5.10). This can be done using a simple matrix-vector multiplication. If
we evaluate 𝐻(𝑒 𝑗𝜔 ) at k= 0, 1,..., K equispaced frequencies over [0,π], then

∑𝑀𝑚=0 𝑏𝑚 𝑒
−𝑗𝜔𝑘 𝑚
𝐻(𝑒 𝑗𝜔𝑘 ) = , 𝑘 = 0, 1, … . , 𝐾 (5.11)
1 + ∑𝑁
𝑖=0 𝑎𝑖 𝑒
−𝑗𝜔𝑘 𝑙

If we let {𝑏𝑚 }, {𝑎𝑙 } (with 𝑎0 = 1), {m = 0, …,M}, {l = 0, …, N}, and {𝜔𝑘 }be arrays (or row vectors), then
the numerator and the denominator of (6.11) become

𝑏𝑒𝑥𝑝(−𝑗𝑚𝑇 𝜔); 𝑎𝑒𝑥𝑝(𝑗𝑙 𝑇 𝜔)

respectively. Now the array 𝐻(𝑒 𝑗𝜔𝑘 ) in (5.11) can be computed using a ./ operation. This procedure can
be implemented in a MATLAB function to determine the frequency response function, given {𝑏𝑚 } and
{𝑎𝑙 } arrays.

4G. A 3rd-order low-pass filter is described by the difference equation

Laboratory 4: The Discrete-time Fourier Analysis 4.7

𝑦(𝑛) = 0.0181𝑥(𝑛) + 0.0543𝑥(𝑛 − 1) + 0.0543𝑥(𝑛 − 2) + 0.0181𝑥(𝑛 − 3) + 1.76𝑦(𝑛 − 1)

− 1.1829𝑦(𝑛 − 2) + 0.2781𝑦(𝑛 − 3)

We will plot the magnitude and the phase response of this filter

>> b = [0.0181, 0.0543, 0.0543, 0.0181];

>> a = [1.0000, -1.7600, 1.1829, -0.2781];
>> m = 0:length(b)-1; l = 0:length(a)-1;
>> K = 500; k = 0:1:K;
>> w = pi*k/K;
>> num = b * exp(-j*m'*w);
>> den = a * exp(-j*l'*w);
>> H = num ./ den;
>> magH = abs(H); angH = angle(H);
>> subplot(2,1,1); plot(w/pi,magH); grid; axis([0,1,0,1])
>> xlabel('frequency in pi units'); ylabel(' | H | ');
>> title('Magnitude Response');
>> subplot(2,1,2); plot(w/pi,angH/pi); grid
>> xlabel('frequency in pi units'); ylabel('Phase in pi Radians');
>> title('Phase Response');

Main Lab
4H: Let 𝑥(𝑛) be a random sequence uniformly distributed between [0, 1] over 0 ≤ n ≤ 10 and let 𝑦(𝑛) =
𝑥(𝑛 − 2). Verify the sample shift property.

4I: Verify the folding property, let 𝑥(𝑛) be a random sequence over -5 ≤ n ≤ 10 uniformly distributed
between [0, 1].r -π ≤ 𝜔 ≤ π.

4J: For the linear time invariant systems described by the impulse response,
ℎ(𝑛) = 𝑠𝑖𝑛𝑐(0.2𝑛)[𝑢(𝑛) − 𝑢(𝑛 − 40)]
determine the frequency response function H (ejw) and plot the magnitude response and the phase
response