Bangladesh University of Engineering and Technology
Department of Electrical and Electronic Engineering
Course No: EEE 312
Course Title: Digital Signal Processing I Laboratory
Experiment No: 03 and 04(Phase I)
Name of the Experiment: (a) Z- Z-transform and Its Application
(b) Frequency Domain Analysis of DT signals and
systems
Submitted to
Shahed Ahmed, Lecturer, EEE, BUET
Md. Obaidur Rahman
Submitted By
Md. Amir Hamza
Student ID:2006145
Department: EEE
Section: C1
Level:3 Term:1
Date of Submission :02/02/2024
Experiment 3
Exercise: C.4
After finding the z-transform of the system from LCCDE
For a stable system, ROC must include a unit circle, so, |z|>0.5
Matlab code for C.4
num = [1 -1/3];
den = [1, -1/2];
% Impulse response
figure;
subplot(2,1,1);
stem(impz(num, den), 'filled');
title('Impulse response');
subplot(2,1,2);
zplane(num, den);
title('Pole-zero plot of $H(z) = \frac{1 - \frac{1}{3}z^{-1}}{1
- 0.5z^{-1}}$, ROC $|z| > 0.5$', 'Interpreter', 'latex');
Plot:
Exercise: C.5
Matlab code for C.5
% Defining the system
numerator = [3, -4];
denominator = [1, -3.5, 1.5];
% Plotting pole-zero diagram
figure;
subplot(4, 1, 1);
zplane(numerator, denominator);
title('Pole-Zero Plot of $H(z) = \frac{3 - 4z^{-1}}{1 - 3z^{-1}
- 1.5z^{-2}}$', 'Interpreter', 'latex');
% Calculating residues and poles
[R, p, C] = residuez(numerator, denominator);
% Causal
n_causal = 0:10;
h_causal = R(1) * p(1).^n_causal + R(2) * p(2).^n_causal;
subplot(4, 1, 2);
stem(n_causal, h_causal, 'filled');
title('Causal with ROC $|z| > 3$', 'Interpreter', 'latex');
% Anticausal
n_anticausal = -10:-1;
h_anticausal = -R(1) * p(1).^n_anticausal - R(2) *
p(2).^n_anticausal;
subplot(4, 1, 3);
stem(n_anticausal, h_anticausal, 'filled');
title('Anticausal with ROC $|z| < 0.5$', 'Interpreter',
'latex');
% Noncausal
n_noncausal = -10:10;
h1 = R(2) * p(2).^n_noncausal(n_noncausal >= 0);
h2 = R(1) * p(1).^n_noncausal(n_noncausal < 0);
h = [h2, h1];
subplot(4, 1, 4);
stem(n_noncausal, h, 'filled');
title('Noncausal with ROC $0.5 < |z| < 3$', 'Interpreter',
'latex');
Plot:
Exercise: C.6
After using residuez() function, the partial fraction of the system is
ROC:
Matlab Code for C.6
% Defining the system
numerator = [2, 16, 44, 56, 31];
denominator = [3, 3, -15, 18, -12];
% Calculating residues and poles
[R, p, C] = residuez(numerator, denominator);
% Plotting pole-zero diagram
figure;
zplane(numerator, denominator);
title('Pole-Zero Plot of $G(z) = \frac{2 + 16z^{-1} + 44z^{-2} +
56z^{-3} + 31z^{-4}}{3 + 3z^{-1} - 15z^{-2} + 18z^{-3} - 12z^{-
4}}$', 'Interpreter', 'latex');
Plot:
Exercise: C.7
Transfer function for the system is derived from the given
Difference equation is
Matlab code for C.7
% Defining the system
numerator = [1];
denominator = [1, -5/6, 1/6];
% Defining the input signal
n = 0:20;
x = zeros(1, length(n));
x(1) = 1;
x(2) = -1/3;
% Calculating the system response
y = filter(numerator, denominator, x);
% Plotting pole-zero diagram
figure;
subplot(3, 1, 1);
zplane(numerator, denominator);
title('Pole-Zero Plot of $H(z) = \frac{1}{1 - \frac{5}{6}z^{-1} + \frac{1}
{6}z^{-2}}$ ROC $|z| > 1$', 'Interpreter', 'latex');
% Plotting input and output response
subplot(3, 1, 2);
stem(n, y, 'filled', 'DisplayName', 'y');
hold on;
stem(n, x, 'filled', 'DisplayName', 'x');
legend;
title('Input and Output Response');
% Reconstruct the system response from the output
h = zeros(1, length(n));
h(1) = y(1);
for i = 2:length(n)
h(i) = y(i) + (1/3) * h(i-1);
end
% Plotting reconstructed system response
subplot(3, 1, 3);
stem(n, h, 'filled');
title('Reconstructed System Response from Output');
Plot:
Exercise: C.8
The system function has poles at 0.5 and 2. To stabilize the system, we can
eliminate the poles that correspond to a specific causality condition. For the
system to be both causal and stable, all the poles must be strictly within the
unit circle, with a region of convergence (ROC) defined by |z| > r, where r is
less than 1. This means we can combine H(z) with another system,
H1(z) = 1 - 2z^(-1)
to achieve this. Similarly, for the system to be anti-causal and stable, all the
poles must be outside the unit circle, with an ROC defined by |z| < r, where r
is greater than 1. In this case, we can combine H(z) with another system,
H2(z) = 1 - 0.5z^(-1).
Matlab Code for C.8
% Defining the system
numerator = [1];
denominator = [1, -2.5, 1];
% Plotting pole-zero diagram
figure;
subplot(3, 2, 1);
zplane(numerator, denominator);
title('Pole-Zero Plot of $H(z) = \frac{1}{1 + 2.5z^{-1} - z^{-2}}$ ROC $|z| >
2$', 'Interpreter', 'latex');
% Plotting corresponding impulse response
subplot(3, 2, 2);
stem(impz(numerator, denominator));
title('Corresponding Impulse Response');
% Assuming the system to be causal and canceling the pole at 2
numerator1 = conv(numerator, [1, -2]);
subplot(3, 2, 3);
zplane(numerator1, denominator);
title('Pole-Zero Plot of $H(z) = \frac{1(1 - 2z^{-1})}{1 + 2.5z^{-1} - z^{-
2}}$ ROC $|z| > 0.5$', 'Interpreter', 'latex');
% Plotting corresponding causal impulse response
subplot(3, 2, 4);
stem(impz(numerator1, denominator));
title('Corresponding Causal Impulse Response');
% Assuming the system to be anti-causal and canceling the pole at 1/2
numerator2 = conv(numerator, [1, -1/2]);
subplot(3, 2, 5);
zplane(numerator2, denominator);
title('Pole-Zero Plot of $H(z) = \frac{1(1 - 0.5z^{-1})}{1 + 2.5z^{-1} - z^{-
2}}$ ROC $|z| < 2$', 'Interpreter', 'latex');
% Plotting corresponding anti-causal impulse response
subplot(3, 2, 6);
n = -20:-1;
h = -2.^n;
stem(n, h);
title('Corresponding Anticausal Impulse Response');
Plot:
Schur_Cohn_stability test
This MATLAB program is designed to test the stability of a system by using
the Schur-Cohn stability test. It takes the coefficients of the denominator of
a system's transfer function as input. When an unstable system function is
tested using the Schur-Cohn stability test, it is found that not all K values are
less than zero. This result confirms that the input system function is indeed
unstable.
Matlab Code : Schur Cohn Stability Test
% Defining the coefficients of the polynomial
coefficients = [2, 3, -1, 0.3, 1, 0.4];
% Plotting the pole-zero diagram
figure;
zplane(1, coefficients);
% Checking stability using Schur-Cohn stability test
stabilityResult = SchurCohnStability(coefficients);
function stability = SchurCohnStability(coefficients)
% Initialize variables
n = length(coefficients);
K = zeros(1, n);
% Performing Schur-Cohn stability test
for m = n:-1:1
coefficients = coefficients / coefficients(1);
K(m) = coefficients(end);
coefficients = (coefficients - K(m) * fliplr(coefficients)) / (1 -
K(m)^2);
coefficients = coefficients(1:end-1);
end
% Checking if the system is stable
if all(K < 0)
stability = 'The system is stable.';
else
stability = 'The system is unstable.';
end
end
Plot:
Experiment: 4(Phase I)
For problem 1-4 the index for the discrete time signal will be
n= -145:145
using this index performing the code:
Problem 1 & 2:
We're computing the discrete Fourier transform (DFT) of a
finite-length (-145:145) sinc signal. Because the DFT is a
periodic signal with a period of 2π, adjusting the interval of the
frequency grid highlights this periodic behavior.
Matlab Code
% Defining the frequency variable for the sinc pulse
f_c = 1/8;
% Defining the index for the sinc pulse
n = -145:145;
% Generating the sinc pulse
x = sinc(2 * f_c * n);
% Creating a figure with three subplots
figure;
% Plotting the sinc pulse x(n) against n
subplot(3, 1, 1);
stem(n, x, '.');
xlabel('n');
ylabel('x(n)');
title('sinc(2f_c*n) pulse');
% Defining the number of points in the digital frequency grid
M = 201;
% Defining the frequency ranges
w_list = {[0, 2], [-2, 2]};
% Loop through each frequency range
for i = 1:length(w_list)
% Creating a digital frequency grid
w = linspace(pi * w_list{i}(1), pi * w_list{i}(2), M);
% Calculating the resolution of the digital frequency
dw = w(2) - w(1);
% Initialize the DTFT of x(n)
X = zeros(1, M);
% Calculating the DTFT of x(n) for each frequency point in w
for i1 = 1:M
for i2 = 1:length(x)
X(i1) = X(i1) + x(i2) * exp(-1i * w(i1) * n(i2));
end
end
% Plotting the magnitude of X(w) against the normalized frequency w/pi
subplot(3, 1, i + 1);
plot(w / pi, abs(X));
xlabel('Frequency (w/pi) (rad/sec)');
ylabel('X(w)');
title(sprintf('Frequency Spectra for w = [%d pi, %d pi]', w_list{i}(1),
w_list{i}(2)));
end
Plot:
Problem 3 & 4:
We set M = 201 as the count of points in the DTFT frequency
grid, ranging from 0 to π. In our MATLAB code, we're effectively
performing a DFT operation to reconstruct the original signal
with a finely spaced frequency grid, aiming to approximate the
continuous frequency nature of DTFT. Consequently, when we
choose a wider range for x[n], surpassing the number of
frequency grid points, the reconstructed signal seems to exhibit
periodic behavior with a period of M.
Matlab Code
% Defining the number of points in the digital frequency grid
M = 201; % Initializing the number of points in the digital frequency grid
w = linspace(0, 2*pi, M); % Creating a linearly spaced vector of digital
frequencies
dw = w(2) - w(1); % Calculating the frequency resolution
% Defining a list of ranges for n_re
n_re_list = {[-145, 100], [-145, 80], [-240, 240]}; % Creating a cell array
of ranges for n_re
% Iterating through each range in n_re_list
for i = 1:length(n_re_list)
n_re = n_re_list{i}(1):n_re_list{i}(2); % Extracting the current range
for n_re
x_re = zeros(1, length(n_re)); % Initializing the reconstructed signal
vector
% Iterating through each point in the digital frequency grid
for i1 = 1:M
% Iterating through each point in the current range for n_re
for i2 = 1:length(x_re)
% Calculating the reconstructed signal at the current point in
n_re
x_re(i2) = x_re(i2) + 1/(2*pi) * X(i1) * exp(-1i * w(i1) *
n_re(i2)) * dw;
end
end
% Plotting the reconstructed signal for the current range of n_re
figure(2); % Creating a new figure
subplot(3, 1, i); % Creating a subplot
stem(n_re, x_re, '.'); % Plotting the reconstructed signal
xlabel('n'); ylabel('$x_{re}(n)$'); % Setting the axis labels
title(sprintf('$x_{re}(n)$ with index n = %d:%d', n_re_list{i}(1),
n_re_list{i}(2))); % Setting the title
end
Plot:
From the above figure, we can easily understand the periodicity
of the reconstructed signal.
Problem 5: DTFT of a CT signal
Matlab Code
% Setting the sampling frequency
Fs = 2e3; % Sampling Frequency
dt = 1 / Fs; % Time step
% Defining signal parameters
f = 100; % Frequency of the sinusoidal pulse
t = -0.02:dt:0.02; % Time interval
N = length(t); % Number of samples
n = t / dt; % Index array
x = zeros(1, N); % Initializing the signal vector
% Generating the sinusoidal pulse
pulse_indices = abs(n) <= 0.005 / dt; % Indices where the pulse is active
x(pulse_indices) = sin(2 * pi * f * t(pulse_indices)); % Creating the
sinusoidal pulse
% Plotting the continuous and sampled signals
figure(3);
subplot(221);
plot(t, x);
xlabel('t (s)'), ylabel('x(t)');
title('Continuous-time signal');
subplot(222);
stem(n, x, '.');
xlabel('n'), ylabel('x(n)');
title(sprintf('Sampled signal (Fs = %d Hz)', Fs));
% Calculating the Discrete Fourier Transform (DFT)
M = 500; % Total number of points in the frequency range
w = (-M/2:M/2) * 2 * pi / M; % Frequency grid
W = exp(-1i * w' * n); % Matrix formation
X = W * x'; % Calculating the DFT
% Plotting the magnitude spectrum
subplot(223);
plot(w * Fs, abs(X));
xlabel('Analog frequency, $\omega$ (rad/sec)');
ylabel('$|X(\omega)|$');
title('Magnitude Spectrum');
% Plotting the phase spectrum
subplot(224);
plot(w * Fs, angle(X) * 180 / pi);
xlabel('Analog frequency, $\omega$ (rad/sec)');
ylabel('angle ($X(\omega)$)');
title('Phase Spectrum');
Plot:
Problem 6: Duality of DTFT
Matlab Code
% Setting the sampling frequency and time step
Fs = 50e3; % Sampling frequency
dt = 1 / Fs; % Time step
% Defining the pulse parameters
T = 1e-3; % Period of the pulse
t = -T/2 + dt:dt:T/2; % Time interval for a period
n = t / dt; % Index for data points in a period
% Initializing the signals
x_rect = zeros(1, length(t));
x_sinc = zeros(1, length(t));
% Generating a single rectangular pulse
D = 0.5; % Duty cycle
PW = D * T; % Pulse width
x_rect(abs(n) <= (PW / dt) / 2) = 1.1; % Generating a single rectangular
pulse
% Generating a single sinc pulse
f_sinc = 0.5e4;
x_sinc = sinc(2 * f_sinc * t); % Generating a single sinc pulse
% Plotting the discrete-time signals
figure(4);
subplot(221);
stem(n, x_rect, '.');
xlabel('n');
ylabel('x(n)');
title('Discrete-time rectangular signal');
subplot(223);
stem(n, x_sinc, '.');
xlabel('n');
ylabel('x(n)');
title('Discrete-time sinc signal');
% Calculating the Discrete Time Fourier Transform (DTFT)
M = 500; % Total number of points in the frequency range
w = (-M/2:M/2) * 2 * pi / M; % Frequency grid
W = exp(-1i * w' * n); % Matrix formation
X_rect = W * x_rect';
X_sinc = W * x_sinc';
% Plotting the DTFT
subplot(222);
plot(w / pi, X_rect);
xlabel('Digital frequency, $\omega / \pi$');
ylabel('$X(\omega)$');
title('DTFT of rectangular pulse');
subplot(224);
plot(w / pi, X_sinc);
xlabel('Digital frequency, $\omega / \pi$');
ylabel('$X(\omega)$');
title('DTFT of sinc pulse');
Plot:
