Discrete Signal Processing (DSP):
Linear Time-Invariant (LTI) Systems
Table of Contents
1. Linear Time-Invariant (LTI) Systems¶
A discrete-time system $H$ is a transformation that maps an input signal $x[n]$ to an output signal $y[n]$:
$$ y[n] = H\{x[n]\} $$
System Examples
Below are some typical system behaviors:
$$ \begin{aligned} \text{Identity} &\qquad y[n] = x[n] \\ \text{Scaling} &\qquad y[n] = 2x[n] \\ \text{Offset} &\qquad y[n] = x[n] + 2 \\ \text{Nonlinearity} &\qquad y[n] = (x[n])^2 \\ \text{Shift} &\qquad y[n] = x[n + 2] \\ \text{Decimation} &\qquad y[n] = x[2n] \\ \text{Index Squaring} &\qquad y[n] = x[n^2] \end{aligned} $$
Each of these systems modifies the input signal according to a specific rule. The remaining sections will focus on those systems that are linear and time-invariant, as they are fundamental to signal processing theory.
1.1. Linear Systems¶
A system $H$ is linear if it satisfies:
(1) Scaling:
$$ H\{\alpha x[n]\} = \alpha H\{x[n]\}, \quad \forall \alpha \in \mathbb{C} $$
(2) Additivity:
If $ y_1 [n] = H\{x_1 [n]\} $ and $ y_2 [n] = H\{x_2 [n]\} $, then
$$ H\{x_1[n] + x_2[n]\} = y_1[n] + y_2[n] $$
Matrix Multiplication and Linear Systems
Matrix multiplication serves as a fundamental computational model for representing linear systems in signal processing. It provides a compact and expressive way to encode how input signals are linearly transformed into output signals.
Let $H$ denote a matrix representing the linear system, and let $x$ be the input vector (i.e., the signal). The output $y$ of the system is given by the matrix-vector product:
$$ \begin{aligned} y &= Hx \\\\ y[n] &= \sum_m H_{n,m} \, x[m] = \sum_m h_{n,m} \, x[m] \end{aligned} $$
Here, $h_{n,m} = H_{n,m}$ is the element in the $n$-th row and $m$-th column of the matrix $H$, indicating how input $x[m]$ contributes to output $y[n]$.
Surprisingly, any linear system can be fully described by a matrix multiplication. The matrix $H$ captures the influence of each input component on each output component.
Two Interpretations of Matrix Multiplication
Matrix multiplication can be understood through two complementary perspectives:
(1) Column-wise view
The output $y$ is a linear combination of the columns of $H$, weighted by the entries of the input vector $x$:
$$ y = x[0] \cdot H[:,0] + x[1] \cdot H[:,1] + \cdots + x[M-1] \cdot H[:,M-1] $$
This emphasizes how each column of the matrix contributes to the final output.
(2) Row-wise view
Each output component $y[n]$ is the inner product of the $n$-th row of $H$ with the input vector $x$:
$$ y[n] = \langle H[n,:],\, x \rangle $$
This highlights the computation of each output sample as a weighted sum of the input signal.
These two perspectives are essential for both understanding the structure of linear systems and implementing them efficiently in practice.
1.2. Time-Invariant Systems¶
A system $H$ is called time-invariant (or shift-invariant) if a shift in the input signal results in an identical shift in the output signal. This property ensures that the system's behavior does not depend on the absolute timing of the input. It reacts the same way regardless of when a signal is applied.
(1) Time-Invariance for Infinite-Length Signals
For systems operating on infinite-length (i.e., unbounded) discrete-time signals, time-invariance is formally defined as:
$$ H\{x[n - q]\} = y[n - q] \quad \text{for all integers } q $$
This means that if the input signal is delayed by $q$ samples, the output will also be delayed by $q$ samples, with no change in shape or amplitude.
(2) Time-Invariance for Finite-Length Signals
For systems that operate on signals of fixed length $N$, indexing is handled modulo-$N$ to preserve signal length after shifts. In this case, the time-invariance condition is expressed as:
$$ H\{x[(n - q)_N]\} = y[(n - q)_N] $$
Here, $(n - q)_N$ denotes modulo-$N$ circular indexing, which wraps the delayed indices around the signal length $N$.
This circular shift formulation is essential when defining circular convolution and performing system analysis in the frequency domain using the Discrete Fourier Transform (DFT).
Understanding time-invariance under both linear and circular shifting is crucial for designing and analyzing digital filters, especially when working with periodic or block-based signal processing systems.
1.3. Linear Time-Invariant (LTI) Systems¶
A system $H$ is said to be linear time-invariant (LTI) if it satisfies both linearity and time-invariance properties. These systems form the foundation of digital signal processing and enjoy powerful analytical tools. In this course, we will primarily focus on LTI systems.
Examples of System Properties
$$ \begin{aligned} \text{Identity} &\qquad y[n] = x[n] &\quad \text{Linear} &\quad \text{Time-Invariant} \\ \text{Scaling} &\qquad y[n] = 2\,x[n] &\quad \text{Linear} &\quad \text{Time-Invariant} \\ \text{Offset} &\qquad y[n] = x[n] + 2 &\quad \text{Nonlinear} &\quad \text{Time-Invariant} \\ \text{Square Signal} &\qquad y[n] = (x[n])^2 &\quad \text{Nonlinear} &\quad \text{Time-Invariant} \\ \text{Shift} &\qquad y[n] = x[n+2] &\quad \text{Linear} &\quad \text{Time-Invariant} \\ \text{Decimation} &\qquad y[n] = x[2n] &\quad \text{Linear} &\quad \text{Time-Variant} \\ \text{Nonlinear Index} &\qquad y[n] = x[n^2] &\quad \text{Linear} &\quad \text{Time-Variant} \end{aligned} $$
1.3.1. LTI Systems as Matrix Multiplications (Infinite-Length Signals)¶
For infinite-length signals, a general linear system is expressed as:
$$ y[n] = H\{x[n]\} = \sum_{m=-\infty}^{\infty} h_{n,m} \, x[m] $$
Time-invariance implies:
$$ H\{x[n - q]\} = y[n - q] = \sum_{m=-\infty}^{\infty} h_{n,m} \, x[m - q] $$
Using a change of variables: $n' = n - q$, $m' = m - q$
$$ H\{x[n']\} = \sum_{m=-\infty}^{\infty} h_{n'+q, m'+q} \, x[m'] $$
Matching expressions reveals the structure:
$$ h_{n,m} = h_{n+q, m+q} $$
$$H = \begin{bmatrix}&\vdots&\vdots&\vdots&\\ \cdots&h_{-1,-1}&h_{-1,0}&h_{-1,1}&\cdots\\ \cdots&h_{0,-1}&h_{0,0}&h_{0,1}&\cdots\\ \cdots&h_{1,-1}&h_{1,0}&h_{1,1}&\cdots\\ &\vdots&\vdots&\vdots&\\ \end{bmatrix}= \begin{bmatrix}&\vdots&\vdots&\vdots&\\ \cdots&h_{0,0}&h_{-1,0}&h_{-2,0}&\cdots\\ \cdots&h_{1,0}&h_{0,0}&h_{-1,0}&\cdots\\ \cdots&h_{2,0}&h_{1,0}&h_{0,0}&\cdots\\ &\vdots&\vdots&\vdots&\\ \end{bmatrix}$$
This implies $h_{n,m}$ depends only on the difference $n - m$, so we write:
$$ h_{n,m} = h[n - m] $$
The matrix representation becomes Toeplitz (constant along diagonals):
- The entire matrix is determined by the impulse response $h[n] = h_{n,0}$
- Diagonals are constant: $[H]_{n,m} = h[n - m]$
$$H = \begin{bmatrix} \ddots &\vdots&\vdots&\vdots&\cdots\\ \cdots&h_{0,0}&h_{-1,0}&h_{-2,0}&\cdots\\ \cdots&h_{1,0}&h_{0,0}&h_{-1,0}&\cdots\\ \cdots&h_{2,0}&h_{1,0}&h_{0,0}&\cdots\\ &\vdots&\vdots&\vdots&\ddots\\ \end{bmatrix}= \begin{bmatrix} \ddots & \vdots & \vdots & \vdots & \cdots \\ \cdots & h[0] & h[-1] & h[-2] & \cdots \\ \cdots & h[1] & h[0] & h[-1] & \cdots \\ \cdots & h[2] & h[1] & h[0] & \cdots \\ \cdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} $$
All of the entries in a Toeplitz matrix can be expressed in terms of the entries of its 0-th column:
$$h[n] = h_{n,0}$$
This is a remarkable observation: it aligns perfectly with the foundational fact that a linear time-invariant (LTI) system is completely characterized by its impulse response $h[n]$. Once the impulse response is known, we can construct the full system matrix $H$, which governs how the system acts on any input signal.
1.3.2. LTI Systems for Finite-Length Signals¶
For signals of fixed length $N$, the linear system is expressed as:
$$ y[n] = H\{x[n]\} = \sum_{m=0}^{N-1} h_{n,m} \, x[m] $$
Time-invariance for circular (modulo-$N$) signals implies:
$$ H\{x[(n - q)_N]\} = y[(n - q)_N] = \sum_{m=0}^{N-1} h_{n,m} \, x[(m - q)_N] $$
Change of variables: $n' = n - q$, $m' = m - q$
$$ H\{x[n']\} = \sum_{m=0}^{N-1} h_{(n'+q)_N, (m'+q)_N} \, x[m'] $$
Thus,
$$ h_{n,m} = h_{(n + q)_N, (m + q)_N} $$
$$H = \begin{bmatrix}h_{0,0}&h_{0,1}&h_{0,2}&\cdots&h_{0,N-1}\\ h_{1,0}&h_{1,1}&h_{1,2}&\cdots&h_{1,N-1}\\ h_{2,0}&h_{2,1}&h_{2,2}&\cdots&h_{2,N-1}\\ \vdots&\vdots&\vdots&&\vdots\\ h_{N-1,0}&h_{N-1,1}&h_{N-1,2}&\cdots&h_{N-1,N-1}\\ \end{bmatrix}= \begin{bmatrix}h_{0,0}&h_{N-1,0}&h_{N-2,0}&\cdots&h_{1,0}\\ h_{1,0}&h_{0,0}&h_{N-1,0}&\cdots&h_{2,0}\\ h_{2,0}&h_{1,0}&h_{0,0}&\cdots&h_{3,0}\\ \vdots&\vdots&\vdots&&\vdots\\ h_{N-1,0}&h_{N-2,0}&h_{N-3,0}&\cdots&h_{0,0}\\ \end{bmatrix}$$
This implies the structure:
$$ h_{n,m} = h[(n - m)_N] $$
The matrix $H$ is a circulant matrix, fully determined by the first column:
$$H = \begin{bmatrix}h_{0,0}&h_{N-1,0}&h_{N-2,0}&\cdots&h_{1,0}\\ h_{1,0}&h_{0,0}&h_{N-1,0}&\cdots&h_{2,0}\\ h_{2,0}&h_{1,0}&h_{0,0}&\cdots&h_{3,0}\\ \vdots&\vdots&\vdots&&\vdots\\ h_{N-1,0}&h_{N-2,0}&h_{N-3,0}&\cdots&h_{0,0}\\ \end{bmatrix}= \begin{bmatrix}h[0]&h[N-1]&h[N-2]&\cdots&h[1]\\ h[1]&h[0]&h[N-1]&\cdots&h[2]\\ h[2]&h[1]&h[0]&\cdots&h[3]\\ \vdots&\vdots&\vdots&&\vdots\\ h[N-1]&h[N-2]&h[N-3]&\cdots&h[0]\\ \end{bmatrix}$$
Row-$n$, column-$m$ entry is:
$$ [H]_{n,m} = h[(n - m)_N] $$
Notice the structure of the diagonals and the circular shifts.
This structure is fundamental to circular convolution, and its properties enable fast implementation via the Discrete Fourier Transform (DFT).
2. Impulse Response¶
In the study of linear time-invariant (LTI) systems, the impulse response plays a fundamental role. It provides a complete characterization of the system's behavior and allows us to construct the system matrix $H$ or compute the system's output via convolution.
Definition of the Impulse Signal
The discrete-time impulse signal $\delta[n]$, also called the delta function, is defined as:
$$ \delta[n] = \begin{cases} 1 & \text{if } n = 0 \\ 0 & \text{otherwise} \end{cases} $$
It serves as a probing signal that isolates the response of the system at a single time index.
2.1. Impulse Response of an LTI System¶
Let $H$ denote a discrete-time LTI system. The impulse response of the system, denoted $h[n]$, is the output when the input is $\delta[n]$:
$$ h[n] = H\{\delta[n]\} $$
This signal $h[n]$ contains all information about the system's behavior and can be used to determine the system output for any arbitrary input via convolution.
For linear time-invariant (LTI) systems, the impulse response uniquely characterizes the system. Specifically:
- The system output for any input $x[n]$ is given by the convolution:
$$ y[n] = \sum_{m=-\infty}^{\infty} h[n - m]\, x[m] $$
- When the input is $\delta[n]$, the output reduces to:
This confirms that $h[n]$ is the system's response to the impulse input.
Moreover, the impulse response determines the matrix representation of the system:
- The matrix $H$ can be constructed from $h[n]$, with each row (or column) being a shifted version of the impulse response
- For infinite-length signals, LTI systems correspond to Toeplitz matrices, where the entries satisfy:
$$ [H]_{n,m} = h[n - m] $$
That is, all elements along any given diagonal are equal.
Hence, the impulse response encapsulates all the information necessary to describe an LTI system, both analytically (via convolution) and structurally (via the system matrix).
2.2. Examples (Infinite-Length Signals)¶
In this example, we will compute the system matrix $H$ and, more importantly, visualize it in order to gain deeper insight into the impulse responses of various linear time-invariant (LTI) systems.
$$H= \begin{bmatrix}&\vdots&\color{red}{\vdots}&\vdots&\\ \cdots&h[0]&\color{red}{h[-1]}&h[-2]&\cdots\\ \cdots&h[1]&\color{red}{h[0]}&h[-1]&\cdots\\ \cdots&h[2]&\color{red}{h[1]}&h[0]&\cdots\\ &\vdots&\color{red}{\vdots}&\vdots&\\ \end{bmatrix} $$
(1) Impulse response of the scaling system
$$ \begin{align*} h[n] &= 2\delta[n]\\\\ H_{n,m} &= h[n-m] \end{align*} $$
%plot -s 560,200
N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];
x = zeros(size(nx))';
x(N+1) = 1;
h = zeros(size(nh))';
h(1) = 2;
subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,2])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,2])
%plot -s 560,420
hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);
imagesc(H); colormap('jet'); axis('square'); axis off
(2) Impulse response of the shift system
$$ \begin{align*} h[n] &= \delta[n-3]\\\\ H_{n,m} &= h[n-m] \end{align*} $$
%plot -s 560,200
N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];
x = zeros(size(nx))';
x(N+1) = 1;
h = zeros(size(nh))';
h(4) = 1;
subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])
%plot -s 560,420
hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);
imagesc(H); colormap('jet'); axis('square'); axis off
(3) Impulse response of the moving average system
$$ \begin{align*} h[n] &= \frac{1}{3}(\delta[n] + \delta[n-1] + \delta[n-2])\\\\ H_{n,m} = h[n-m] &= \frac{1}{3}(\delta[n-m] + \delta[n-m-1] + \delta[n-m-2]) \end{align*} $$
%plot -s 560,200
N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];
x = zeros(size(nx))';
x(N+1) = 1;
h = zeros(size(nh))';
h(1) = 1/3;
h(2) = 1/3;
h(3) = 1/3;
subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])
%plot -s 560,420
hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);
imagesc(H); colormap('jet'); axis('square'); axis off
(4) Impulse response of the recursive average system
$$y[n] = x[n] + \alpha y[n-1]$$
$$ \begin{align*} h[n] &= \alpha^n u[n]\\\\ H_{n,m} = h[n-m] &= \alpha^{n-m}u[n-m] \end{align*} $$
%plot -s 560,200
N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];
x = zeros(size(nx))';
x(N+1) = 1;
h = 0.8.^nh;
h = h';
subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])
%plot -s 560,420
hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);
imagesc(H); colormap('jet'); axis('square'); axis off
3.3. Examples (Finite-Length Signals)¶
We now repeat the visualization for finite-length signals. In this case, the resulting system matrix exhibits constant values along its diagonals, together with circular wraparound, which is characteristic of a circulant matrix.
$$h_{n,m} = h_{(n+q)_N (m+q)_N} \qquad \forall q \in \mathbb{Z}$$
$$ H = \begin{bmatrix}\color{red}{h[0]}&h[N-1]&h[N-2]&\cdots&h[1]\\ \color{red}{h[1]}&h[0]&h[N-1]&\cdots&h[2]\\ \color{red}{h[2]}&h[1]&h[0]&\cdots&h[3]\\ \vdots&\vdots&\vdots&&\vdots\\ \color{red}{h[N-1]}&h[N-2]&h[N-3]&\cdots&h[0]\\ \end{bmatrix} $$
(1) Impulse response of the shift system
$$ \begin{align*} h[n] &= \delta[n-2]\\\\ H_{n,m} &= h[(n-m)_N] \end{align*} $$
%plot -s 560,200
N = 8;
nh = [0:N-1];
nx = [0:N-1];
x = zeros(size(nx))';
x(1) = 1;
h = zeros(size(nh))';
h(3) = 1;
subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])
%plot -s 560,420
H = zeros(N,N);
for k = 0:N-1
H(:,k+1) = circshift(h,k);
end
imagesc(H); colormap('jet'); axis('square'); axis off
(2) Impulse response of the moving average system
$$ \begin{align*} h[n] &= \frac{1}{3}(\delta[n] + \delta[n-1] + \delta[n-2])\\\\ H_{n,m} &= h[(n-m)_N] \end{align*} $$
%plot -s 560,200
N = 8;
nh = [0:N-1];
nx = [0:N-1];
x = zeros(size(nx))';
x(1) = 1;
h = zeros(size(nh))';
h(1) = 1/3;
h(2) = 1/3;
h(3) = 1/3;
subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])
%plot -s 560,420
H = zeros(N,N);
for k = 0:N-1
H(:,k+1) = circshift(h,k);
end
imagesc(H); colormap('jet'); axis('square'); axis off
(3) Impulse response of the recursive average system
$$y[n] = x[n] + \alpha y[n-1]$$
$$ \begin{align*} h[n] &= \alpha^n u[n]\\\\ H_{n,m} &= h[(n-m)_N] \end{align*} $$
%plot -s 560,200
N = 8;
nh = [0:N-1];
nx = [0:N-1];
x = zeros(size(nx))';
x(1) = 1;
h = zeros(size(nh))';
h = 0.8.^nh;
subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])
%plot -s 560,420
H = zeros(N,N);
for k = 0:N-1
H(:,k+1) = circshift(h,k);
end
imagesc(H); colormap('jet'); axis('square'); axis off
3. Convolution¶
Convolution is a fundamental operation in linear systems theory. For discrete-time signals, convolution describes the output of an LTI system as the weighted sum of shifted and reversed versions of the input signal.
The convolution of an input signal $x[n]$ with a system’s impulse response $h[n]$ is defined as:
$$ y[n] = \sum_{m = -\infty}^{\infty} h[n - m]\, x[m] = x[n] * h[n] $$
This expression means that the output $y[n]$ is obtained by convolving the input $x[n]$ with the system's impulse response $h[n]$.
Commutativity: Convolution is commutative, which means the order of operands does not matter:
$$ x * h = h * x $$
This can be seen by a change of variables in the convolution sum:
$$ \begin{aligned} y[n] &= \sum_{m = -\infty}^{\infty} h[n - m]\, x[m] \\ &= \sum_{k = -\infty}^{\infty} h[k]\, x[n - k] \qquad \text{(where } k = n - m \text{)} \\ &= \sum_{k = -\infty}^{\infty} x[n - k]\, h[k] \end{aligned} $$
Thus, convolution is symmetric with respect to the input signal and the system impulse response: we may think of either as the "signal" and the other as the "system".
(MATLAB Code) Convolution
$$\begin{align*} x[n]=\{3,11,7,&\;0,-1,4,2\}\quad \quad -3\leq n \leq 3\\ &\uparrow \\\\ h[n]=\{2,&\;3,0,-5,2,1\}\quad \quad -1\leq n \leq 4\\ &\uparrow \end{align*} $$
%plot -s 560,200
%% matlab command 'conv'
x = [3,11,7,0,-1,4,2];
h = [2,3,0,-5,2,1];
y = conv(x,h);
stem(y,'filled'), ylim([-60 60])
% this is not correct
You have to include conv_m function file in the path.
function [y,ny] = conv_m(x,nx,h,nh)
% Modified convolution routine for signal processing
% [y,ny] = conv_m(x,nx,h,nh)
% y = convolution result
% ny = support of y
% x = first signal on support nx
% nx = support of x
% h = second signal on support nh
% nh = support of h
nyb = nx(1) + nh(1);
nye = nx(length(x)) + nh(length(h));
ny = [nyb:nye];
y = conv(x,h);
$$\begin{align*} x[n]=\{3,11,7,&\;0,-1,4,2\}\quad \quad -3\leq n \leq 3\\ &\uparrow \\\\ h[n]=\{2,&\;3,0,-5,2,1\}\quad \quad -1\leq n \leq 4\\ &\uparrow \end{align*} $$
%plot -s 560,200
x = [3,11,7,0,-1,4,2]; nx = [-3:3];
h = [2,3,0,-5,2,1]; nh = [-1:4];
[y,ny] = conv_m(x,nx,h,nh);
stem(ny,y,'filled'), axis tight, ylim([-60 60])
% this is what we want
%plot -s 560,600
subplot(3,1,1), stem(nx,x,'filled'); xlim([ny(1),ny(end)]); ylabel('x')
subplot(3,1,2), stem(nh,h,'filled'); xlim([ny(1),ny(end)]); ylabel('h')
subplot(3,1,3), stem(ny,y,'filled'); xlim([ny(1),ny(end)]); ylabel('y')
3.1. Convolution for Infinite-Length Signals¶
For systems operating on infinite-length discrete-time signals, convolution is mathematically equivalent to multiplying the input signal by a Toeplitz matrix.
The output signal $y[n]$ is given by the convolution sum:
$$ \begin{aligned} y[n] &= (x * h)[n] = \sum_{m = -\infty}^{\infty} h[n - m]\, x[m] \\\\ &= \cdots + h[n + 2]\, x[-2] + h[n + 1]\, x[-1] + h[n]\, x[0] + h[n - 1]\, x[1] + h[n - 2]\, x[2] + \cdots \end{aligned} $$
This expression represents an inner product between a shifted (and flipped) version of the impulse response $h$ and the input signal $x$. For each $n$, the vector $h[n - m]$ aligns with $x[m]$, and their pointwise product is summed to compute $y[n]$.
In matrix notation, the convolution can be written as:
$$ y[n] = H x $$
where $H$ is a Toeplitz matrix whose rows are shifted versions of the impulse response:
$$ \begin{align*} y[n] &= \begin{bmatrix} \cdots & \mid & \mid & \mid & \mid & \mid & \cdots\\ \cdots & h[n+2] & h[n+1] & h[n] & h[n-1] & h[n-2]\ & \cdots\\ \cdots & \mid & \mid & \mid & \mid & \mid & \cdots \end{bmatrix} \begin{bmatrix} \vdots\\ x[-2]\\ x[-1]\\ x[0]\\ x[1]\\ x[2]\\ \vdots \end{bmatrix}\\ \\ &= H x \end{align*} $$
Thus, convolution with an LTI system corresponds to multiplying the input vector $x$ by a structured matrix $H$ whose diagonal-constant structure reflects the system's time-invariance.
Convolution is product of matrix $H$ and $x$.
x = [3, 11, 7, 0, -1, 4, 2]';
h = [2, 3, 0, -5, 2, 1]';
hc = [h,zeros([length(x)-1,1])];
hr = [h(1),zeros([1,length(x)-1])];
H = toeplitz(hc,hr)
y = H*x
2.3. Convolution for Finite-Length Signals¶
In the case of finite-length signals of length $N$, convolution is performed using circular convolution, which assumes periodic (modulo-$N$) signal boundaries.
The circular convolution of two length-$N$ signals $x[n]$ and $h[n]$ is defined as:
$$ y[n] = x[n] \otimes h[n] = \sum_{m = 0}^{N-1} h\left[(n - m)_N\right] \, x[m] $$
Here, $(n - m)_N$ denotes modular arithmetic (i.e., reduction modulo $N$) to ensure the indices wrap around circularly.
Interpretation:
- The impulse response $h[n]$ is first circularly time-reversed.
- It is then circularly shifted by $n$ steps to the right.
- The output $y[n]$ is computed as the inner product of this shifted version of $h$ with the input vector $x$.
Circular convolution naturally arises in systems where periodic boundary conditions are assumed, and it is central to frequency-domain signal processing using the Discrete Fourier Transform (DFT).
%plot -s 560,200
N = 8;
n = 0:N-1;
h = [0 1 2 3 4 5 6 7]';
stem(n,h,'filled');
%% Create a circulent matrix H based on h
H = zeros(N,N);
for k = 0:N-1
H(:,k+1) = circshift(h,k);
end
H
%plot -s 800,300
subplot(1,2,1), imagesc(H); colormap('jet'); axis('square'); axis off
subplot(1,2,2), imagesc(h); colormap('jet'); axis equal; axis off
x = [1,1,0,0,0,0,0,0]';
y1 = H*x
x = [1,1,0,0,0,0,0,0];
h = [0,1,2,3,4,5,6,7];
y3 = cconv(x,h,8)
x1 = [1,2,2];
x2 = [1,2,3,4];
y = cconv(x1,x2,4)
% 5-point circular convolution
y = cconv(x1,x2,5)
% 6-point circular convolution
y = cconv(x1,x2,6)
%plot -s 560,500
N = 8;
xn = 0:N-1;
x = [1 0 -1 0 1 0 -1 0];
h = [0 1 2 3 4 0 0 0];
y = cconv(x,h,N);
subplot(3,1,1), stem(xn,x,'filled'), ylim([-2 4]), ylabel('x')
subplot(3,1,2), stem(xn,h,'filled'), ylim([-2 4]), ylabel('h')
subplot(3,1,3), stem(xn,y,'filled'), ylim([-2 4]), ylabel('y')
%plot -s 560,600
N = 8;
theta2 = (0:1/N:1)*2*pi;
theta = theta2(1:end-1);
re = cos(theta);
im = sin(theta);
subplot(3,1,1), stem3(re,im,x,'filled'); hold on;
plot(re,im,'k')
quiver(0,0,1,0,'k','maxheadsize',0.5)
quiver(0,0,0,1,'k','maxheadsize',0.5), hold off
xlabel('Real','fontsize',8), ylabel('Imaginary','fontsize',8), zlabel('x[n]')
zlim([-2,4])
subplot(3,1,2), stem3(re,im,h,'filled'); hold on;
plot(re,im,'k')
quiver(0,0,1,0,'k','maxheadsize',0.5)
quiver(0,0,0,1,'k','maxheadsize',0.5), hold off
xlabel('Real','fontsize',8), ylabel('Imaginary','fontsize',8), zlabel('h[n]')
zlim([-2,4])
subplot(3,1,3), stem3(re,im,y,'filled'); hold on;
plot(re,im,'k')
quiver(0,0,1,0,'k','maxheadsize',0.5)
quiver(0,0,0,1,'k','maxheadsize',0.5), hold off
xlabel('Real','fontsize',8), ylabel('Imaginary','fontsize',8), zlabel('y[n]')
zlim([-2,4])
4. Convolution Examples¶
Convolution is a versatile operation that can model a wide range of signal processing tasks. In this section, we explore concrete examples that demonstrate how convolution can be applied in practice. Each example illustrates a specific functionality of linear time-invariant systems, implemented through carefully designed impulse responses.
4.1. Denoising¶
Denoising is a fundamental signal processing task in which convolution is used to suppress noise while retaining the underlying structure of a signal. A common method for denoising is to apply a low-pass filter, which attenuates high-frequency components (often corresponding to noise) and preserves low-frequency trends.
One simple yet effective low-pass filter is the moving average filter. A 5-point moving average filter has the following impulse response:
$$ h[n] = \frac{1}{5} \left( \delta[n] + \delta[n - 1] + \delta[n - 2] + \delta[n - 3] + \delta[n - 4] \right) $$
Equivalently, the filter coefficients can be expressed explicitly as:
$$ h[0] = h[1] = h[2] = h[3] = h[4] = \frac{1}{5}, \quad h[n] = 0 \;\text{otherwise} $$
Convolving a noisy signal $x[n]$ with this impulse response smooths the signal by replacing each value with the average of itself and its four preceding neighbors. This operation can be written as:
$$ y[n] = (x * h)[n] = \frac{1}{5} \sum_{k=0}^{4} x[n - k] $$
This filter reduces the influence of abrupt changes or outliers in the signal, making it particularly effective for removing additive white noise. However, there is a trade-off: the stronger the smoothing effect, the more the fine details of the original signal may be blurred.
The 5-point moving average filter provides a good balance between noise reduction and detail preservation in many practical scenarios.
%plot -s 560,800
N = 50;
n = 0:N-1;
s = sin(pi/N*n) .* [ones(1,N/2), -ones(1,N/2)];
subplot(4,1,1)
stem(s,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal','fontsize',8)
% add noise to the signal
x = s + 0.1*randn(1,N);
subplot(4,1,2)
stem(x,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal + noise','fontsize',8)
% construct moving average filter impulse response of length M
M = 5;
h = ones(1,M)/M;
subplot(4,1,3)
stem(h,'b','Marker','none','LineWidth',2);
axis([1 N 0 1])
title('impulse response','fontsize',8)
% convolve noisy signal with impulse response
y = cconv(x,h,N);
subplot(4,1,4)
stem(y,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('output','fontsize',8)
4.2. Edge Detection¶
Edge detection is the process of identifying abrupt transitions or discontinuities in a signal. These transitions often correspond to important structural features, such as object boundaries in images or sudden changes in temporal signals.
A common approach is to use a high-pass filter that responds strongly to rapid changes while suppressing smooth or constant regions.
Consider the following impulse response:
$$ h[n] = \frac{1}{4} \, [-1, -1, 1, 1] $$
This filter is symmetric and spans four points. It computes the difference between the average of two points on the left and the average of two points on the right. When convolved with a signal $x[n]$, the output is:
$$ y[n] = \frac{1}{4} \left( -x[n - 1] - x[n - 2] + x[n + 0] + x[n + 1] \right) $$
This highlights regions where the signal changes from low to high or vice versa. Specifically:
- Flat or slowly varying regions yield near-zero output
- Sharp transitions produce strong positive or negative responses, depending on the direction of change
This filter effectively computes a local contrast across four adjacent samples:
- The negative weights on the left side suppress past values
- The positive weights on the right emphasize future values
This makes it sensitive to edges centered between the second and third elements of the window.
%plot -s 560,500
N = 50;
n = 0:N-1;
s = sin(pi/N*n) .* [ones(1,N/2), -ones(1,N/2)];
% add noise to the signal
x = s + 0.1*randn(1,N);
subplot(311)
stem(x,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal + noise','fontsize',8)
% haar wavelet edge detector
h = (1/4)*[-ones(1,2), ones(1,2)];
subplot(312)
stem(h,'b','Marker','none','LineWidth',2);
axis([1 N -0.5 0.5])
title('Haar wavelet edge detector','fontsize',8)
% % convolve noisy signal with impulse response
y = cconv(x,h,N);
subplot(313)
stem(y,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('output','fontsize',8)
4.3. Convolution on Audio¶
(Source: "ThinkDSP" by Prof. Allen Downey at Olin College)
We consider two real-world audio signals for this example:
Download wave file 1
Download wave file 2
Simulating Reverberation
To simulate how a violin might sound in a reverberant environment, we apply convolution between the violin signal and a gunshot waveform. The gunshot serves as a synthetic impulse response that characterizes the acoustic behavior of the simulated environment.
The convolution is expressed as:
$$ y[n] = x[n] * h[n] $$
where:
- $x[n]$ is the input signal (violin),
- $h[n]$ is the impulse response (gunshot),
- $y[n]$ is the output signal with simulated reverberation.
This example highlights a key application of convolution in digital audio: the ability to model and synthesize realistic acoustic effects by applying an impulse response to a clean signal. Such techniques are widely used in virtual acoustics, sound design, and music production.
[x, Fs] = audioread([pwd,'\data_files\violin_origional.wav']);
x = x/max(x); % normalize
sound(x, Fs); % play a wave file with sampling rate Fs
%plot -s 560,200
% plot wave file
xn = (1:length(x))/Fs;
plot(xn,x), xlabel('time in sec')
% impulse response in a closed room (by gunshot)
[h, Fs] = audioread([pwd,'\data_files\gunshot.wav']);
h = h(:,1); % stereo to mono
sound(h, Fs);
plot((1:length(h))/Fs, h), xlabel('sec')
% image how the music played in a closed room sounds like
y = cconv(x,h);
y = y/max(y);
sound(y, Fs);
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')