Convolution as Matrix Multiplication

Signals and Linear Operations

Convolution is the fundamental operation of linear time-invariant (LTI) systems. Given an input signal x and an impulse response (filter) h, their convolution produces an output signal y. Mathematically, (x * h)[n] = Σₖ x[k] h[n−k].

Convolution is a linear operation: it is linear in its inputs — double the input, double the output; sum of inputs gives sum of outputs. This is a crucial observation because it means convolution can be written as a matrix-vector multiplication. Wherever you see a linear operation on a finite vector, there is a matrix hiding behind it.

Core Insight

Every linear operation from ℝⁿ to ℝᵐ is equivalent to multiplication by some m×n matrix. Convolution is linear, so it must have a matrix. The particular structure of convolution — the filter coefficients appearing repeatedly — determines the specific shape of that matrix.

The Toeplitz Matrix

When you write out the convolution of a length-N signal x with a length-M filter h, you see that each output sample y[n] is a dot product of x with a shifted copy of h. Collecting all output samples into a matrix equation reveals the Toeplitz matrix H — a matrix where every diagonal (top-left to bottom-right) has a constant value.

For a filter h = [h₀, h₁, h₂, ...], the Toeplitz matrix H looks like:

Main diagonal: all h₀
First sub-diagonal: all h₁
Second sub-diagonal: all h₂
... and so on

The convolution output is simply the matrix-vector product y = Hx. Each row of H is a shifted copy of the filter h, with zeros padded as needed.

Toeplitz Convolution

\mathbf{y} = H\mathbf{x}, \quad H_{ij} = h[i-j]

The convolution y = x * h is equivalent to the matrix-vector product y = Hx, where H is a Toeplitz matrix. Each diagonal of H holds the same filter coefficient — the defining property of a Toeplitz matrix. Direct computation requires O(N²) multiplications for length-N signals.

Circular Convolution and Circulant Matrices

Linear (aperiodic) convolution with zero-padding gives Toeplitz matrices. But if we impose periodic boundary conditions — treating the signal as if it wraps around — we get circular convolution, and its matrix is a circulant matrix.

In a circulant matrix, each column is a cyclic shift of the column to its left. The first column completely determines the entire matrix. For example, if the first column is [h₀, h₁, h₂, h₃]ᵀ, the second column becomes [h₃, h₀, h₁, h₂]ᵀ, the third [h₂, h₃, h₀, h₁]ᵀ, and so on. Entrywise, C_ij = h[(i−j) mod N] — the same rule as the Toeplitz matrix above, now wrapped around.

Circulant matrices have a remarkable algebraic property: they share the same set of eigenvectors regardless of the specific values in the first column. Those eigenvectors are the DFT basis vectors — the complex exponentials e^(j2πkn/N) for k = 0, 1, ..., N−1.

Why Circular Convolution?

Circular convolution arises naturally when working with periodic signals and the DFT. To use FFT-based methods for linear convolution, we zero-pad signals to sufficient length and then perform circular convolution, which gives the same result as linear convolution in the non-wrapped region.

DFT Diagonalizes Circulant Matrices

Because the DFT basis vectors are the eigenvectors of every circulant matrix, the DFT matrix F simultaneously diagonalizes all circulant matrices. This is the deepest reason why convolution becomes pointwise multiplication in the frequency domain.

Let C be an N×N circulant matrix and F be the N-point DFT matrix (whose (k,n) entry is ω^(kn)/√N where ω = e^(−j2π/N)). Then C admits the eigendecomposition:

DFT Diagonalization

C = F^{-1} \Lambda F

Every circulant matrix C is diagonalized by the DFT matrix F. The diagonal matrix Λ contains the eigenvalues of C on its diagonal, which are exactly the DFT of the first column of C — equivalently, the DFT of the filter h. F⁻¹ = F* (the conjugate transpose), since F is unitary. Note the order: it is F⁻¹ΛF, not FΛF⁻¹. Building the circulant from the first column is what puts the inverse on the left; the row-shifted matrix is Cᵀ, which performs circular correlation rather than convolution.

The eigenvalues λₖ = H[k] are just the frequency-domain representation of the filter. So circular convolution Cx = y becomes, in the frequency domain, Λ(Fx) = Fy — which is pointwise multiplication of the DFT of x by the DFT of h, giving the DFT of y.

FFT-Based Fast Convolution

Direct computation of the N-point convolution via matrix multiplication requires O(N²) operations. But the DFT diagonalization shows us a faster route: compute the DFTs of x and h, multiply them pointwise, then take the inverse DFT. Using the FFT algorithm, each DFT requires only O(N log N) operations.

FFT-Based Convolution

\mathbf{y} = \text{IFFT}\!\left(\text{FFT}(\mathbf{x}) \odot \text{FFT}(\mathbf{h})\right)

The symbol ⊙ denotes pointwise (element-wise) multiplication. Three FFT-scale operations (two forward FFTs and one inverse FFT) replace the O(N²) matrix multiplication, reducing total complexity to O(N log N). For large N this is an enormous speedup: at N = 10⁶, the ratio N/log₂N ≈ 50,000×.

This is one of the most impactful algorithm design insights in all of engineering. The FFT-based convolution underpins audio and video codecs, wireless communication signal processing, image filtering, radar pulse compression, gravitational wave detection, and countless other applications.

Overlap-Add and Overlap-Save

When the input signal x is very long (or of unknown/infinite length), it is impractical to apply a single giant FFT. The overlap-add and overlap-save methods partition the input into manageable blocks and combine the results to achieve the same output as a single linear convolution.

Overlap-Add

Segment x into non-overlapping blocks x₁, x₂, x₃, .... Zero-pad each block to length N + M − 1 (where N is the block length and M−1 is the filter order). Compute the circular convolution of each block with h via FFT, producing output blocks that overlap by M−1 samples. Add the overlapping tails of adjacent output blocks to reconstruct the full linear convolution.

Overlap-Save

Collect input blocks of length N that overlap the previous block by M−1 samples. Compute the circular convolution of each block with h via FFT. Discard the first M−1 samples of each output block (which are corrupted by circular wrap-around), and save only the valid N − M + 1 samples. Concatenate the saved portions to recover the linear convolution.

Both methods achieve the same O(N log N) per-output-sample complexity as the direct FFT approach, but require only O(N) memory at any given time — critical for real-time processing of streaming signals.

Key Takeaways

Convolution is a linear operation, so it corresponds to multiplication by a matrix — specifically a Toeplitz matrix (shifted filter coefficients on each diagonal). Circular convolution corresponds to a circulant matrix, which has the special property that the DFT matrix diagonalizes it: C = F⁻¹ΛF. The eigenvalues Λ are the DFT of the filter h. This diagonalization is the algebraic foundation of why convolution equals pointwise multiplication in the frequency domain. FFT-based convolution exploits this to reduce complexity from O(N²) to O(N log N). Overlap-add and overlap-save extend this approach to long or streaming signals, maintaining efficiency while keeping memory usage bounded.