Filter Design with Linear Algebra

The FIR Filter as a Matrix-Vector Product

A finite impulse response (FIR) filter with coefficients h = [h₀, h₁, ..., h_{M−1}] computes each output sample as a weighted sum of the M most recent input samples. For a block of N input samples collected into a vector x, the entire output vector y is produced by the matrix-vector product y = Hx, where H is the Toeplitz convolution matrix.

This matrix perspective reframes filter design as a question: which matrix H — equivalently, which coefficient vector h — produces the best output? "Best" depends on our criterion. The most tractable and widely used criterion is least squares: minimize the sum of squared differences between the actual output and some desired output.

Why Linear Algebra?

Once the filter is cast as a matrix, the entire machinery of linear algebra — projections, orthogonality, eigendecompositions — can be brought to bear on filter design. Optimal filters emerge as solutions to linear systems, not as ad-hoc design recipes.

Specifying What We Want: The Desired Response

In supervised filter design, we have access to a desired signal d[n] — what the output should ideally be — and an input signal x[n]. For each time index n, we stack M past input samples into the regressor vector:

Regressor Vector

\mathbf{x}[n] = [x[n],\, x[n-1],\, \ldots,\, x[n-M+1]]^T \in \mathbb{R}^M

The regressor x[n] is the sliding window of M past input samples at time n. The filter output is the inner product y[n] = hᵀx[n], where h is the coefficient vector. Assembling all N regressors as rows of a matrix X gives the batch problem: y = Xh.

Collecting N pairs of regressors and desired samples into the matrix X (N×M) and vector d (N×1), the goal is to find the coefficient vector h that makes Xh as close to d as possible in the least-squares sense.

Least Squares Filter Design

The least squares problem minimizes the total squared error between the filter output and the desired signal over N samples:

Least Squares Criterion

J(\mathbf{h}) = \|\mathbf{d} - X\mathbf{h}\|^2 = \sum_{n=0}^{N-1}\bigl(d[n] - \mathbf{h}^T\mathbf{x}[n]\bigr)^2

The cost J(h) is a quadratic function of h (a paraboloid in M-dimensional space), so it has a unique global minimum when XᵀX is invertible — which requires N ≥ M and sufficiently rich input signals.

Setting the gradient ∂J/∂h = 0 yields the normal equations — the hallmark of least-squares problems. These are M linear equations in M unknowns, with a positive semidefinite coefficient matrix.

Normal Equations

(X^T X)\,\mathbf{h}^* = X^T \mathbf{d}

XᵀX is the M×M empirical autocorrelation (Gram) matrix, and Xᵀd is the M×1 cross-correlation vector. The solution h* = (XᵀX)⁻¹Xᵀd is the least-squares filter — the orthogonal projection of d onto the column space of X.

Geometric Interpretation

The least-squares solution projects the desired vector d onto the subspace spanned by the columns of X (i.e., the set of all achievable filter outputs). The error vector d − Xh* is orthogonal to every column of X — the hallmark of orthogonal projection.

The Wiener Filter: Optimal in MSE Sense

When signals are modeled as stationary random processes, the optimal filter minimizes the mean squared error (MSE) E[|d[n] − y[n]|²]. Taking expectations replaces the empirical matrices with their statistical counterparts:

R = E[x[n]xᵀ[n]] — the M×M autocorrelation matrix of the input
p = E[d[n]x[n]] — the M×1 cross-correlation vector between the desired signal and the input

The MSE-optimal filter satisfies the Wiener-Hopf equation:

Wiener-Hopf Equation

R\,\mathbf{h}_{\mathrm{opt}} = \mathbf{p}\quad\Longrightarrow\quad \mathbf{h}_{\mathrm{opt}} = R^{-1}\mathbf{p}

R is a symmetric positive definite Toeplitz matrix (since the input is stationary). Its Toeplitz structure can be exploited by the Levinson-Durbin algorithm to solve for h_opt in O(M²) operations rather than the O(M³) of general Gaussian elimination. The minimum achievable MSE is ξ_min = σ_d² − pᵀh_opt, where σ_d² = E[d²[n]].

The Wiener filter is the statistical analog of the least-squares filter — the two coincide in the limit of large N. In practice, the true statistics R and p are unknown and must be estimated from data, leading to the sample (data-driven) Wiener filter.

Eigenstructure of the Autocorrelation Matrix

Because R is symmetric positive definite, it admits an eigendecomposition R = QΛQᵀ. The eigenvectors Q are the principal directions of the input signal's power distribution, and the eigenvalues λ₁ ≥ λ₂ ≥ ... ≥ λ_M > 0 are the powers in those directions. The Wiener filter solution can be expressed in the eigenbasis as:

Eigen-Domain Wiener Filter

\mathbf{h}_{\mathrm{opt}} = \sum_{k=1}^{M} \frac{\mathbf{q}_k^T \mathbf{p}}{\lambda_k}\,\mathbf{q}_k

In the eigenbasis, each component of the Wiener filter scales the projection of p onto the k-th eigenvector by 1/λₖ. Eigenvectors with small eigenvalues (low input power) contribute large components to h_opt, making the filter sensitive to estimation noise — the "ill-conditioning" problem motivating regularization.

Adaptive Filters: LMS and RLS

In practice, signal statistics change over time (non-stationarity), or we may have insufficient data to estimate R and p reliably. Adaptive filters update their coefficients online, tracking changing statistics without needing to store or invert large matrices.

Least Mean Squares (LMS)

The LMS algorithm approximates the gradient of the MSE cost with an instantaneous estimate, using a single input-error sample pair to update h:

LMS Update Rule

\mathbf{h}[n+1] = \mathbf{h}[n] + \mu\, e[n]\,\mathbf{x}[n]

e[n] = d[n] − hᵀ[n]x[n] is the instantaneous error. μ is the step size (learning rate). Each update costs only O(M) — no matrix inversions. LMS converges in the mean if 0 < μ < 2/λ_max, where λ_max is the largest eigenvalue of R. The condition number κ(R) = λ_max/λ_min determines convergence speed: ill-conditioned R leads to slow convergence.

Recursive Least Squares (RLS)

RLS minimizes the weighted sum of all past squared errors, updating the inverse autocorrelation matrix P = (XᵀX)⁻¹ recursively using the matrix inversion lemma (Sherman-Morrison-Woodbury formula). RLS converges in exactly M steps (for exact arithmetic) and tracks non-stationarity much faster than LMS, at the cost of O(M²) operations per update instead of O(M).

Beamforming as a Linear Algebra Problem

An array of K antennas receives the same signal from a direction θ, each with a different phase shift. The received vector at time n is z[n] = a(θ)s[n] + n[n], where a(θ) is the steering vector (the array response to a signal from direction θ), s[n] is the desired signal, and n[n] is noise plus interference.

A beamformer applies a weight vector w to the antenna array output: y[n] = wᴴz[n]. The goal is to choose w to pass the signal from direction θ₀ while suppressing noise and interference from other directions.

MVDR Beamformer

\min_{\mathbf{w}}\; \mathbf{w}^H R_z \mathbf{w} \quad \text{subject to} \quad \mathbf{w}^H \mathbf{a}(\theta_0) = 1

The MVDR (Minimum Variance Distortionless Response) beamformer minimizes the output power (noise + interference) subject to maintaining unit gain in the look direction. The solution is w_opt = R_z⁻¹a / (aᴴR_z⁻¹a), where R_z = E[zz ᴴ] is the array covariance matrix. This requires inverting the K×K covariance matrix — a core linear algebra operation.

Beamforming is thus a constrained least-squares problem: minimize quadratic form wᴴR_zw subject to the linear constraint wᴴa = 1. The solution follows directly from Lagrange multipliers and the matrix inverse — a beautiful application of the linear algebra we've built throughout this course.

Key Takeaways

FIR filters are Toeplitz matrix-vector products; designing them optimally means solving a least-squares problem via the normal equations (XᵀX)h = Xᵀd. The statistical analog is the Wiener filter, given by the Wiener-Hopf equation Rh = p, where R is the input autocorrelation matrix and p is the cross-correlation with the desired signal. The eigenstructure of R governs both filter performance and conditioning. Adaptive algorithms (LMS, RLS) track changing statistics online: LMS approximates the gradient at O(M) cost; RLS uses the matrix inversion lemma at O(M²) cost for exact least-squares updates. Beamforming is a constrained quadratic optimization over antenna weights, solvable with a single matrix inverse — the MVDR beamformer.