MATLAB Tutorial - ECOM 451 Digital Signal Processing

1. Getting Started with DSP in MATLAB

🎯 Course Learning Outcomes Addressed

CLO-1: Describe discrete-time signals and systems
CLO-2: Apply z-Transform for DSP system analysis
CLO-3: Apply DTFT in signal analysis
CLO-4: Apply DFT in signal analysis
CLO-5 & 6: Design FIR and IIR digital filters

Essential MATLAB Setup

Basic DSP Setup

% Clear workspace and figures
clear all; close all; clc;

% Define sampling parameters
Fs = 1000;           % Sampling frequency (Hz)
T = 1/Fs;            % Sampling period
N = 1024;            % Number of samples
n = 0:N-1;           % Sample indices
t = n*T;             % Time vector

% Create a test signal
f1 = 50;             % Signal frequency (Hz)
f2 = 120;            % Second frequency
x = sin(2*pi*f1*t) + 0.5*sin(2*pi*f2*t);

% Plot the signal
figure('Position', [100, 100, 800, 400]);
plot(t, x, 'b-', 'LineWidth', 1.5);
xlabel('Time (s)'); ylabel('Amplitude');
title('Test Signal: 50 Hz + 120 Hz Sinusoids');
grid on;

Key DSP Toolbox Functions

Function	Description	Example
`fft()`	Fast Fourier Transform	`X = fft(x)`
`ifft()`	Inverse FFT	`x = ifft(X)`
`freqz()`	Frequency response	`freqz(b, a)`
`filter()`	Apply filter	`y = filter(b, a, x)`
`conv()`	Convolution	`y = conv(x, h)`
`butter()`	Butterworth filter	`[b,a] = butter(n, Wn)`
`fir1()`	FIR filter design	`b = fir1(n, Wn)`

2. Discrete-Time Signals

Basic Signal Generation

Common Discrete Signals

%% Unit Impulse (Delta Function)
n = -10:10;
delta = (n == 0);   % 1 when n=0, 0 otherwise

figure;
stem(n, delta, 'b', 'LineWidth', 2, 'MarkerFaceColor', 'b');
xlabel('n'); ylabel('\delta[n]');
title('Unit Impulse \delta[n]');
grid on;

%% Unit Step Function
u = (n >= 0);       % 1 when n>=0, 0 otherwise

figure;
stem(n, u, 'r', 'LineWidth', 2, 'MarkerFaceColor', 'r');
xlabel('n'); ylabel('u[n]');
title('Unit Step u[n]');
grid on;

%% Exponential Sequence
a = 0.9;            % Decay factor (|a| < 1 for decay)
x_exp = a.^n .* (n >= 0);

figure;
stem(n, x_exp, 'g', 'LineWidth', 2, 'MarkerFaceColor', 'g');
xlabel('n'); ylabel('x[n]');
title(['Exponential Sequence: x[n] = (', num2str(a), ')^n u[n]']);
grid on;

Sinusoidal Sequences

Discrete Sinusoids

%% Discrete Sinusoid
n = 0:50;
omega0 = 0.2*pi;    % Digital frequency (rad/sample)
phi = 0;            % Phase

x = cos(omega0*n + phi);

figure;
subplot(2,1,1);
stem(n, x, 'b', 'LineWidth', 1.5, 'MarkerFaceColor', 'b');
xlabel('n'); ylabel('x[n]');
title(['Discrete Sinusoid: \omega_0 = ', num2str(omega0/pi), '\pi rad/sample']);
grid on;

% Show continuous envelope
subplot(2,1,2);
n_fine = 0:0.1:50;
x_envelope = cos(omega0*n_fine + phi);
plot(n_fine, x_envelope, 'r--', 'LineWidth', 1);
hold on;
stem(n, x, 'b', 'LineWidth', 1.5, 'MarkerFaceColor', 'b');
xlabel('n'); ylabel('x[n]');
title('Discrete Samples of Continuous Sinusoid');
legend('Continuous', 'Discrete');
grid on;

💡 Key Concept: Digital vs Analog Frequency

In DSP, we use digital frequency ω (rad/sample) related to analog frequency f (Hz) by:

ω = 2πf/Fs or ω = 2πfT

where Fs is the sampling frequency and T = 1/Fs is the sampling period.

3. Convolution & Correlation

Linear Convolution

Convolution Example

%% Linear Convolution
% Define two sequences
x = [1, 2, 3, 4];           % Input signal
h = [1, 1, 1];              % Impulse response

% Compute convolution
y = conv(x, h);

% Length of output: length(x) + length(h) - 1
fprintf('Length of x: %d\n', length(x));
fprintf('Length of h: %d\n', length(h));
fprintf('Length of y: %d\n', length(y));

% Plot results
figure('Position', [100, 100, 900, 600]);

subplot(3,1,1);
stem(0:length(x)-1, x, 'b', 'LineWidth', 2, 'MarkerFaceColor', 'b');
xlabel('n'); ylabel('x[n]');
title('Input Signal x[n]');
grid on;

subplot(3,1,2);
stem(0:length(h)-1, h, 'r', 'LineWidth', 2, 'MarkerFaceColor', 'r');
xlabel('n'); ylabel('h[n]');
title('Impulse Response h[n]');
grid on;

subplot(3,1,3);
stem(0:length(y)-1, y, 'g', 'LineWidth', 2, 'MarkerFaceColor', 'g');
xlabel('n'); ylabel('y[n]');
title('Output y[n] = x[n] * h[n] (Convolution)');
grid on;

Correlation

Cross-Correlation and Autocorrelation

%% Cross-Correlation
x = [1, 2, 3, 4, 3, 2, 1];
y = [0, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0];  % Shifted version

% Cross-correlation
[r_xy, lags] = xcorr(x, y);

% Find the lag of maximum correlation
[max_corr, idx] = max(r_xy);
lag_at_max = lags(idx);
fprintf('Maximum correlation at lag = %d\n', lag_at_max);

figure;
subplot(2,1,1);
plot(lags, r_xy, 'b-', 'LineWidth', 2);
hold on;
stem(lag_at_max, max_corr, 'r', 'LineWidth', 2, 'MarkerFaceColor', 'r');
xlabel('Lag'); ylabel('R_{xy}[k]');
title('Cross-Correlation');
grid on;

%% Autocorrelation
[r_xx, lags] = xcorr(x, 'biased');  % Biased estimate

subplot(2,1,2);
stem(lags, r_xx, 'g', 'LineWidth', 2, 'MarkerFaceColor', 'g');
xlabel('Lag'); ylabel('R_{xx}[k]');
title('Autocorrelation of x[n]');
grid on;

4. Discrete-Time Fourier Transform (DTFT)

Computing DTFT

DTFT Implementation

%% DTFT of a finite sequence
x = [1, 2, 3, 4, 3, 2, 1];  % Input sequence
n = 0:length(x)-1;          % Time indices

% Compute DTFT at many frequency points
omega = linspace(-pi, pi, 1000);
X_dtft = zeros(size(omega));

for k = 1:length(omega)
    X_dtft(k) = sum(x .* exp(-1j * omega(k) * n));
end

% Plot magnitude and phase
figure('Position', [100, 100, 800, 600]);

subplot(2,1,1);
plot(omega/pi, abs(X_dtft), 'b-', 'LineWidth', 2);
xlabel('\omega/\pi (rad/sample)');
ylabel('|X(e^{j\omega})|');
title('DTFT Magnitude Spectrum');
grid on;

subplot(2,1,2);
plot(omega/pi, angle(X_dtft), 'r-', 'LineWidth', 2);
xlabel('\omega/\pi (rad/sample)');
ylabel('\angle X(e^{j\omega}) (rad)');
title('DTFT Phase Spectrum');
grid on;

Using freqz() for System Analysis

Frequency Response with freqz()

%% Frequency Response of a Digital Filter
% Moving average filter (FIR)
M = 5;                      % Filter order
b = ones(1, M) / M;         % Numerator coefficients
a = 1;                      % Denominator (FIR)

% Compute frequency response
[H, w] = freqz(b, a, 1024);

figure('Position', [100, 100, 800, 600]);

subplot(2,1,1);
plot(w/pi, 20*log10(abs(H)), 'b-', 'LineWidth', 2);
xlabel('\omega/\pi (rad/sample)');
ylabel('Magnitude (dB)');
title(['Frequency Response: ', num2str(M), '-Point Moving Average Filter']);
grid on;

subplot(2,1,2);
plot(w/pi, unwrap(angle(H))*180/pi, 'r-', 'LineWidth', 2);
xlabel('\omega/\pi (rad/sample)');
ylabel('Phase (degrees)');
title('Phase Response');
grid on;

⚠️ Important: DTFT vs DFT

DTFT: Continuous function of frequency ω, for infinite or finite aperiodic sequences.

DFT: Sampled version of DTFT at N equally spaced frequencies, used for computation.

5. Discrete Fourier Transform (DFT) & FFT

Basic FFT Usage

FFT Spectrum Analysis

%% FFT Spectrum Analysis
Fs = 1000;                  % Sampling frequency
T = 1/Fs;                   % Sampling period
N = 1024;                   % Number of samples
t = (0:N-1)*T;              % Time vector

% Create signal with two frequencies
f1 = 50; f2 = 120;
x = sin(2*pi*f1*t) + 0.6*sin(2*pi*f2*t);

% Add some noise
x_noisy = x + 0.5*randn(size(x));

% Compute FFT
X = fft(x_noisy);

% Frequency axis
f = (0:N-1)*(Fs/N);         % Frequency vector
f_shifted = (-N/2:N/2-1)*(Fs/N);  % Centered frequency

% Two-sided spectrum
X_mag = abs(X)/N;

% One-sided spectrum (for real signals)
X_onesided = X_mag(1:N/2+1);
X_onesided(2:end-1) = 2*X_onesided(2:end-1);
f_onesided = f(1:N/2+1);

figure('Position', [100, 100, 900, 700]);

subplot(3,1,1);
plot(t, x_noisy, 'b-', 'LineWidth', 0.5);
xlabel('Time (s)'); ylabel('Amplitude');
title('Noisy Signal in Time Domain');
grid on;

subplot(3,1,2);
plot(f_onesided, X_onesided, 'r-', 'LineWidth', 1.5);
xlabel('Frequency (Hz)'); ylabel('|X(f)|');
title('One-Sided Amplitude Spectrum');
xlim([0 200]);
grid on;

subplot(3,1,3);
plot(f_shifted, fftshift(X_mag), 'g-', 'LineWidth', 1.5);
xlabel('Frequency (Hz)'); ylabel('|X(f)|');
title('Two-Sided Amplitude Spectrum (Centered)');
xlim([-200 200]);
grid on;

Zero-Padding for Better Resolution

Zero-Padding Example

%% Zero-Padding for Frequency Resolution
x = [1 1 1 1 1];            % Short sequence
N_original = length(x);
N_padded = 64;              % Zero-pad to this length

% FFT without padding
X1 = fft(x);
f1 = (0:N_original-1)/N_original;

% FFT with zero-padding
X2 = fft(x, N_padded);
f2 = (0:N_padded-1)/N_padded;

figure;
subplot(2,1,1);
stem(f1, abs(X1), 'b', 'LineWidth', 2, 'MarkerFaceColor', 'b');
xlabel('\omega/2\pi'); ylabel('|X[k]|');
title(['DFT without Zero-Padding (N = ', num2str(N_original), ')']);
grid on;

subplot(2,1,2);
stem(f2, abs(X2), 'r', 'LineWidth', 1.5, 'MarkerFaceColor', 'r');
xlabel('\omega/2\pi'); ylabel('|X[k]|');
title(['DFT with Zero-Padding (N = ', num2str(N_padded), ')']);
grid on;

💡 FFT Tips

Use N = 2^k (power of 2) for fastest FFT computation
Zero-padding interpolates the spectrum but doesn't add information
Use fftshift() to center the spectrum at DC
For real signals, only the first N/2+1 points are unique

6. z-Transform Analysis

Pole-Zero Analysis

Pole-Zero Plot

%% Pole-Zero Analysis
% Define system: H(z) = (z^2 - 0.5z) / (z^2 - 1.2z + 0.72)
b = [1, -0.5, 0];           % Numerator coefficients
a = [1, -1.2, 0.72];        % Denominator coefficients

% Find zeros and poles
zeros_sys = roots(b);
poles_sys = roots(a);

fprintf('Zeros:\n');
disp(zeros_sys);
fprintf('Poles:\n');
disp(poles_sys);

% Check stability (all poles inside unit circle)
if all(abs(poles_sys) < 1)
    fprintf('System is STABLE (all poles inside unit circle)\n');
else
    fprintf('System is UNSTABLE\n');
end

% Pole-Zero Plot
figure;
zplane(b, a);
title('Pole-Zero Plot');
grid on;

Impulse and Step Response

System Response Analysis

%% Impulse and Step Response
b = [1, 0.5];               % Numerator
a = [1, -0.8, 0.15];        % Denominator

N = 50;                     % Number of samples

% Impulse response
[h, n] = impz(b, a, N);

% Step response
[s, n_step] = stepz(b, a, N);

figure('Position', [100, 100, 800, 600]);

subplot(2,1,1);
stem(n, h, 'b', 'LineWidth', 2, 'MarkerFaceColor', 'b');
xlabel('n'); ylabel('h[n]');
title('Impulse Response');
grid on;

subplot(2,1,2);
stem(n_step, s, 'r', 'LineWidth', 2, 'MarkerFaceColor', 'r');
xlabel('n'); ylabel('s[n]');
title('Step Response');
grid on;

Partial Fraction Expansion

Inverse z-Transform via Partial Fractions

%% Partial Fraction Expansion
% H(z) = (2z^2 + 3z) / (z^2 - 0.5z - 0.5)
b = [2, 3, 0];
a = [1, -0.5, -0.5];

% Partial fraction expansion
[r, p, k] = residuez(b, a);

fprintf('Partial Fraction Expansion:\n');
fprintf('Residues (r):\n');
disp(r);
fprintf('Poles (p):\n');
disp(p);
fprintf('Direct terms (k):\n');
disp(k);

% Reconstruct the sequence
N = 20;
n = 0:N-1;
h_reconstructed = zeros(size(n));

for i = 1:length(r)
    h_reconstructed = h_reconstructed + r(i) * p(i).^n;
end

% Add direct terms
if ~isempty(k)
    h_reconstructed(1:length(k)) = h_reconstructed(1:length(k)) + k;
end

% Compare with impz
h_direct = impz(b, a, N);

figure;
stem(n, h_direct, 'bo-', 'LineWidth', 2, 'DisplayName', 'impz()');
hold on;
stem(n, real(h_reconstructed), 'rx--', 'LineWidth', 1.5, 'DisplayName', 'Partial Fractions');
xlabel('n'); ylabel('h[n]');
title('Impulse Response: Direct vs Partial Fraction');
legend('Location', 'best');
grid on;

7. Digital Filter Design

FIR Filter Design

FIR Filter Using Window Method

%% FIR Lowpass Filter Design
Fs = 1000;                  % Sampling frequency
Fc = 100;                   % Cutoff frequency
N = 50;                     % Filter order

% Normalized cutoff frequency (0 to 1, where 1 = Fs/2)
Wn = Fc / (Fs/2);

% Design FIR filter using different windows
b_rect = fir1(N, Wn, rectwin(N+1));
b_hamm = fir1(N, Wn, hamming(N+1));
b_black = fir1(N, Wn, blackman(N+1));
b_kaiser = fir1(N, Wn, kaiser(N+1, 5));

% Frequency response comparison
[H_rect, w] = freqz(b_rect, 1, 1024, Fs);
[H_hamm, ~] = freqz(b_hamm, 1, 1024, Fs);
[H_black, ~] = freqz(b_black, 1, 1024, Fs);
[H_kaiser, ~] = freqz(b_kaiser, 1, 1024, Fs);

figure('Position', [100, 100, 900, 600]);

subplot(2,1,1);
plot(w, 20*log10(abs(H_rect)), 'b-', 'LineWidth', 1.5); hold on;
plot(w, 20*log10(abs(H_hamm)), 'r-', 'LineWidth', 1.5);
plot(w, 20*log10(abs(H_black)), 'g-', 'LineWidth', 1.5);
plot(w, 20*log10(abs(H_kaiser)), 'm-', 'LineWidth', 1.5);
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
title('FIR Lowpass Filter - Window Comparison');
legend('Rectangular', 'Hamming', 'Blackman', 'Kaiser');
xlim([0 Fs/2]); ylim([-100 5]);
grid on;

subplot(2,1,2);
stem(0:N, b_hamm, 'b', 'LineWidth', 1.5, 'MarkerFaceColor', 'b');
xlabel('n'); ylabel('h[n]');
title('Filter Coefficients (Hamming Window)');
grid on;

IIR Filter Design

Butterworth and Chebyshev Filters

%% IIR Filter Design
Fs = 1000;
Fc = 100;                   % Cutoff frequency
Wn = Fc / (Fs/2);           % Normalized frequency

% Butterworth filter
order_butter = 4;
[b_butter, a_butter] = butter(order_butter, Wn);

% Chebyshev Type I (0.5 dB ripple)
order_cheby = 4;
[b_cheby1, a_cheby1] = cheby1(order_cheby, 0.5, Wn);

% Chebyshev Type II (40 dB stopband attenuation)
[b_cheby2, a_cheby2] = cheby2(order_cheby, 40, Wn);

% Elliptic filter (1 dB passband, 40 dB stopband)
[b_ellip, a_ellip] = ellip(order_cheby, 1, 40, Wn);

% Frequency responses
[H_butter, w] = freqz(b_butter, a_butter, 1024, Fs);
[H_cheby1, ~] = freqz(b_cheby1, a_cheby1, 1024, Fs);
[H_cheby2, ~] = freqz(b_cheby2, a_cheby2, 1024, Fs);
[H_ellip, ~] = freqz(b_ellip, a_ellip, 1024, Fs);

figure('Position', [100, 100, 900, 700]);

subplot(2,1,1);
plot(w, 20*log10(abs(H_butter)), 'b-', 'LineWidth', 2); hold on;
plot(w, 20*log10(abs(H_cheby1)), 'r-', 'LineWidth', 2);
plot(w, 20*log10(abs(H_cheby2)), 'g-', 'LineWidth', 2);
plot(w, 20*log10(abs(H_ellip)), 'm-', 'LineWidth', 2);
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
title('IIR Lowpass Filters Comparison (Order 4)');
legend('Butterworth', 'Chebyshev I', 'Chebyshev II', 'Elliptic');
xlim([0 Fs/2]); ylim([-80 5]);
grid on;

subplot(2,1,2);
zplane(b_butter, a_butter);
title('Butterworth Filter Pole-Zero Plot');
grid on;

Applying Filters to Signals

Signal Filtering Example

%% Filter Application Example
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
N = length(t);

% Create noisy signal: 50 Hz signal + high-frequency noise
signal = sin(2*pi*50*t);
noise = 0.5*sin(2*pi*300*t) + 0.3*randn(size(t));
x = signal + noise;

% Design lowpass filter (cutoff at 100 Hz)
[b, a] = butter(6, 100/(Fs/2));

% Apply filter
y_filter = filter(b, a, x);      % Causal filtering
y_filtfilt = filtfilt(b, a, x);  % Zero-phase filtering

figure('Position', [100, 100, 900, 800]);

subplot(4,1,1);
plot(t, signal, 'g-', 'LineWidth', 1.5);
xlabel('Time (s)'); ylabel('Amplitude');
title('Original Signal (50 Hz)');
grid on;

subplot(4,1,2);
plot(t, x, 'b-', 'LineWidth', 0.5);
xlabel('Time (s)'); ylabel('Amplitude');
title('Noisy Signal');
grid on;

subplot(4,1,3);
plot(t, y_filter, 'r-', 'LineWidth', 1.5);
xlabel('Time (s)'); ylabel('Amplitude');
title('After filter() - Note Phase Delay');
grid on;

subplot(4,1,4);
plot(t, y_filtfilt, 'm-', 'LineWidth', 1.5);
xlabel('Time (s)'); ylabel('Amplitude');
title('After filtfilt() - Zero-Phase Filtering');
grid on;

✅ Filter Design Summary

Filter Type	Advantages	Use When
FIR	Linear phase, always stable	Phase is critical
Butterworth	Maximally flat passband	Smooth response needed
Chebyshev I	Sharper cutoff	Passband ripple OK
Chebyshev II	Flat passband	Stopband ripple OK
Elliptic	Sharpest transition	Both ripples OK

8. Tips & Best Practices

Useful DSP Functions Quick Reference

Signal Generation

sin(), cos()
square(), sawtooth()
randn(), rand()
chirp()

Transform Analysis

fft(), ifft()
fftshift(), ifftshift()
freqz(), phasez()
grpdelay()

Filter Design

fir1(), fir2(), firpm()
butter(), cheby1(), ellip()
filter(), filtfilt()
freqs(), freqz()

System Analysis

tf(), zpk()
impz(), stepz()
zplane(), pzmap()
residuez()

Common Mistakes and Solutions

Error	Cause	Solution
Frequency axis wrong	Not accounting for Fs	Use `f = (0:N-1)*Fs/N`
Cutoff frequency out of range	Wn must be 0-1	Normalize: `Wn = Fc/(Fs/2)`
Unstable IIR filter	Poles outside unit circle	Reduce filter order
Spectral leakage	Non-integer periods in window	Apply window function
Phase distortion	Causal IIR filtering	Use `filtfilt()` for zero-phase

✅ Assignment Submission Checklist

✔ Header with name, ID, assignment number
✔ Code runs without errors from start to finish
✔ All variables properly defined
✔ Sampling frequency and time vectors correct
✔ Figures have titles, axis labels, and legends
✔ Frequency axes properly scaled (Hz or rad/sample)
✔ Results displayed with fprintf()
✔ Code properly commented

📚 Additional Resources

MATLAB Documentation: mathworks.com/help/matlab/
Signal Processing Toolbox: mathworks.com/help/signal/
DSP System Toolbox: mathworks.com/help/dsp/
Textbook: Sanjit K. Mitra, Digital Signal Processing, 4th Ed.
Course Materials: Check Blackboard for sample codes

Need Help?

Contact Dr. Al Bataineh during office hours (Tuesday/Thursday 11:00 AM - 12:00 PM)

Email: mffbataineh@uaeu.ac.ae | Office: F1-1175

← Back to ECOM 451 Course Page

ECOM 451 MATLAB Tutorial

📚 Table of Contents