Understanding Discrete-Time Convolution in MATLAB: A Step-by-Step Guide
Convolution is a fundamental operation in signal processing, used to combine two discrete-time sequences to produce a third sequence that represents how one sequence modifies another. In this blog post, we'll walk through a MATLAB script that computes and visualizes the convolution of two user-defined discrete-time sequences. This guide is designed for students and engineers learning signal processing, offering a clear, hands-on approach to implementing convolution with professional-quality plots. We'll explore the code's functionality, its key components, and the visualization techniques used to make the results intuitive.
What is Discrete-Time Convolution?
Convolution of two discrete-time sequences ( $x[n]$ ) and ( $h[n]$ ) produces an output sequence ( $y[n]$ ), defined as:
$$y[n] = x[n] * h[n] = \sum_{k=-\infty}^{\infty} x[k] h[n-k]$$For finite sequences, the summation is computed over the valid range of indices where both sequences have non-zero values. In MATLAB, the conv function efficiently performs this operation, but we need to handle user inputs, compute correct time indices, and visualize the results effectively.
The MATLAB Script: Overview
The MATLAB script we’ll discuss prompts the user to input two sequences, ( $x[n]$ ) and ( $h[n]$ ), along with their time indices. It computes their convolution ( $y[n]$ ), calculates the appropriate output indices, and generates three stem plots to visualize the input sequences and the convolution result. The plots use LaTeX formatting for labels and titles, and the axes are customized to place the y-axis at the left border with the x-axis zero intersecting it, creating a clean, professional appearance.
Let’s break down the script’s key components, using a Socratic approach to encourage deeper understanding.
Step 1: User Input and Validation
The script begins by clearing the workspace and prompting the user for two sequences:
- First sequence: Values of ( $x[n]$ ) and corresponding time indices ( $nx$ ).
- Second sequence: Values of ( $h[n]$ ) and corresponding time indices ( $nh$ ).
Code Snippet:
disp('Enter the first sequence x[n]:');
x = input('Sequence values (e.g., [1 2 3]): ');
nx = input('Corresponding time indices (e.g., [0 1 2]): ');
if length(x) ~= length(nx)
error('Error: The number of values in x[n] must equal the number of time indices.');
end
The script validates that the number of values matches the number of indices for each sequence, preventing errors in subsequent computations. Why is this validation important? Consider what might happen if the indices and values don’t align—how would it affect the convolution or plotting?
Step 2: Convolution Computation
The convolution is computed using MATLAB’s conv function:
y = conv(x, h);
For sequences of lengths ( $L$ ) (for $x[n]$ ) and ( $M$ ) (for $h[n]$ ), the output ( $y[n]$ ) has length ( $L + M - 1$ ). This reflects the number of valid overlaps as one sequence slides over the other. For example, if ( $x[n] = [1, 2, 3]$ ) (length 3) and ( $h[n] = [1, 1]$ ) (length 2), the convolution output has length ( $3 + 2 - 1 = 4$ ). Can you visualize why the output length is longer than the input sequences?
Step 3: Output Time Indices
To plot $y[n]$ , we need its time indices ( $ny$ ). The script calculates these based on the input indices:
output_length = length(x) + length(h) - 1;
start_index = nx(1) + nh(1);
end_index = start_index + (output_length - 1);
ny = start_index : end_index;
The start index is the sum of the first indices of $x[n]$ and $h[n]$, as this is the earliest point where the convolution produces a non-zero result. The output spans $L + M - 1$ points. For instance, if $nx = [1, 2, 3]$ and $nh = [1, 2]$ , the output starts at $1 + 1 = 2$ and ends at $2 + (3 + 2 - 1 - 1) = 5$ , giving $ny = [2, 3, 4, 5]$ . Why does the start index depend on the first indices of the inputs?
Step 4: Plotting the Sequences
The script creates a figure with three subplots to visualize:
- $x[n]$ - blue stem plot.
- $h[n]$ - red stem plot
- $y[n]$ - black stem plot
Each subplot uses $\text{LaTeX}$ formatting for labels and titles to render mathematical notation professionally:
subplot(3,1,1);
stem(nx, x, 'b', 'filled', 'LineWidth', 1.5);
title('Input Sequence $x[n]$', 'Interpreter', 'latex');
xlabel('$n$', 'Interpreter', 'latex');
ylabel('$x[n]$', 'Interpreter', 'latex');
The stem function is ideal for discrete-time sequences, as it displays values as markers on vertical lines, emphasizing their discrete nature. The 'Interpreter', 'latex' property ensures labels like $n$ and $x[n]$ render with proper subscripts and font styles. Why is LaTeX formatting preferred for mathematical plots?
Step 5: Customizing Axes
To make the plots intuitive, the script customizes the axes:
- Y-axis at the left border:
ax.YAxisLocation = 'left'places the y-axis on the left edge of each subplot. - X-axis at zero:
ax.XAxisLocation = 'origin'aligns the x-axis with ( y = 0 ). - Consistent x-axis limits:
ax.XLim = [min([nx, nh, ny]), max([nx, nh, ny])]ensures all subplots share the same x-axis range, making it easier to compare indices across sequences. - Zero included in y-axis:
ax.YLim = [min([x, h, y, 0]), max([x, h, y, 0])]ensures the y-axis includes zero, aligning the x-axis zero with the y-axis zero at the left border. - Minor grid lines:
set(gca, 'XMinorGrid', 'on')andset(gca, 'YMinorGrid', 'on')add finer grid lines for readability. - Box off:
box offremoves the top and right axes lines for a cleaner look.
These settings ensure the zero points of both axes intersect at the left border, creating a visually appealing and mathematically accurate presentation. What effect does aligning the axes at zero have on interpreting the plots?
Step 6: Adding a Super Title
The script uses sgtitle to add a figure-wide title:
sgtitle('Discrete-Time Convolution', 'Interpreter', 'latex', 'FontSize', 14);
In MATLAB R2024a/b, sgtitle is supported and integrates seamlessly with the subplot layout. The LaTeX interpreter ensures the title renders professionally, though plain text is used here for simplicity. How does a super title enhance the figure’s clarity?
Step 7: Testing the Script
To test the script, try inputs like:
- $x[n] = [1, 2, 3, 4, 5]$ , $nx = [1, 2, 3, 4, 5]$
- $h[n] = [9, 8, 7, 6, 5]$ , $nh = [1, 3, 5, 6, 8]$.
The convolution ( $y[n]$ ) will have length $5 + 5 - 1 = 9$ , with indices ( $ny = [2, 3, ..., 10]$ ). The plots will show:
- $x[n]$ at indices $[1, 2, 3, 4, 5]$.
- $h[n]$ at sparse indices $[1, 3, 5, 6, 8]$.
- $y[n]$ at indices $[2, 3, ..., 10]$.
The axes will intersect at $(0, 0)$ on the left border, with consistent x-axis ranges across subplots. Try running the script with different inputs—how do the output indices and plot ranges adjust?
Learning Takeaways
This script teaches several key concepts:
- Convolution: Understand how two sequences combine to produce a new sequence, with the output length and indices determined by the input lengths and starting points.
- MATLAB Plotting: Learn to create professional stem plots with LaTeX-formatted labels and customized axes.
- Input Handling: Validate user inputs to ensure robust code execution.
- Visualization: Use consistent axis limits and grid lines to make plots intuitive and comparable.
To deepen your understanding, experiment with:
- Consecutive vs. sparse indices for ( nh ).
- Sequences with negative indices or values.
- Adjusting plot properties (e.g., font sizes, colors).
This script is a powerful tool for learning signal processing and MATLAB visualization. Run it, explore the results, and consider how convolution applies to real-world signals like audio or image processing!
Experimental MATLAB code
% Clear workspace and command window
clear all;
clc;
% Input first sequence
disp('Enter the first sequence x[n]:');
x = input('Sequence values (e.g., [1 2 3]): ');
nx = input('Corresponding time indices (e.g., [0 1 2]): ');
% Validate input for first sequence
if length(x) ~= length(nx)
error('Error: The number of values in x[n] must equal the number of time indices.');
end
% Input second sequence
disp('Enter the second sequence h[n]:');
h = input('Sequence values (e.g., [1 1]): ');
nh = input('Corresponding time indices (e.g., [0 1]): ');
% Validate input for second sequence
if length(h) ~= length(nh)
error('Error: The number of values in h[n] must equal the number of time indices.');
end
% Compute convolution
y = conv(x, h);
% Determine time indices for the output sequence
output_length = length(x) + length(h) - 1;
start_index = nx(1) + nh(1);
end_index = start_index + (output_length - 1);
ny = start_index : end_index;
% Create figure
figure('Name', 'Convolution of Discrete-Time Sequences', 'NumberTitle', 'off');
% Plot first sequence
subplot(3,1,1);
stem(nx, x, 'b', 'filled', 'LineWidth', 1.5);
title('Input Sequence $x[n]$', 'Interpreter', 'latex');
xlabel('$n$', 'Interpreter', 'latex');
ylabel('$x[n]$', 'Interpreter', 'latex');
grid on;
ax = gca;
ax.XAxisLocation = 'origin'; % Move x-axis to origin
ax.YAxisLocation = 'left'; % Move y-axis to origin
set(gca,'XMinorGrid','on');
set(gca,'YMinorGrid','on');
ax.XLim = [min([nx-1, nh-1, ny-1,0]), max([nx+1, nh+1, ny+1])]; % Set x-axis limits to include all indices
ax.YLim = [min([x-1,0]), max(x*1.25)]; % Ensure y-axis includes zero
box off; % Remove top and right axes lines
% Plot second sequence
subplot(3,1,2);
stem(nh, h, 'r', 'filled', 'LineWidth', 1.5);
title('Input Sequence $h[n]$', 'Interpreter', 'latex');
xlabel('$n$', 'Interpreter', 'latex');
ylabel('$h[n]$', 'Interpreter', 'latex');
grid on;
ax = gca;
ax.XAxisLocation = 'origin'; % Move x-axis to origin
ax.YAxisLocation = 'left'; % Move y-axis to origin
set(gca,'XMinorGrid','on');
set(gca,'YMinorGrid','on');
ax.XLim = [min([nx-1, nh-1, ny-1,0]), max([nx+1, nh+1, ny+1])]; % Set x-axis limits to include all indices
ax.YLim = [min([h-1,0]), max(h*1.25)]; % Ensure y-axis includes zero
box off; % Remove top and right axes lines
% Plot convolution result
subplot(3,1,3);
stem(ny, y, 'k', 'filled', 'LineWidth', 1.5);
title('Convolution Output $y[n] = x[n] * h[n]$', 'Interpreter', 'latex');
xlabel('$n$', 'Interpreter', 'latex');
ylabel('$y[n]$', 'Interpreter', 'latex');
grid on;
ax = gca;
ax.XAxisLocation = 'origin'; % Move x-axis to origin
ax.YAxisLocation = 'left'; % Move y-axis to origin
set(gca,'XMinorGrid','on');
set(gca,'YMinorGrid','on');
ax.XLim = [min([nx-1, nh-1, ny-1,0]), max([nx+1, nh+1, ny+1])]; % Set x-axis limits to include all indices
ax.YLim = [min([y-1,0]), max(y*1.25)]; % Ensure y-axis includes zero
box off; % Remove top and right axes lines
% Adjust figure size
set(gcf, 'Position', [100, 100, 800, 600]);
% Add LaTeX-formatted super title
sgtitle('Discrete-Time Convolution', 'Interpreter', 'latex', 'FontSize', 14);
Plot Results
