Formula & theorem reference

Sampling f at spacing \Delta t is modeled by multiplying f(t) by a train of unit impulses spaced \Delta t apart, \sum_n \delta(t - n\Delta t). Multiplication in time is (by the Fourier convolution theorem's dual statement) convolution in frequency, and the Fourier transform of an impulse train is itself an impulse train, spaced \omega_s = 2\pi f_s = 2\pi/\Delta t apart in \omega. So sampling replaces \hat f(\omega) with infinitely many shifted copies of itself, \sum_k \hat f(\omega - k\omega_s), spaced \omega_s apart — the spectrum gets periodized.

Since \hat f is supported only on (-2\pi B, 2\pi B), the shifted copies are non-overlapping exactly when the spacing \omega_s = 2\pi f_s exceeds the width 4\pi B of that support, i.e. f_s > 2B. Under that condition, the original \hat f(\omega) can be recovered from the sampled spectrum by an ideal low-pass filter: multiply by a rectangular window that is 1 on (-2\pi B, 2\pi B) (or, more simply, on (-\omega_s/2, \omega_s/2), any gap between copies) and 0 elsewhere, which passes the central copy unchanged and zeroes out every shifted replica. Multiplying by that rectangle in frequency corresponds to convolving with its inverse transform in time, and the inverse transform of an ideal rectangular low-pass filter is exactly a \operatorname{sinc} function scaled to that cutoff — convolving the sample train with this \operatorname{sinc} reproduces f pointwise, which is exactly the stated reconstruction formula.

If instead f_s \le 2B, the shifted copies of \hat f overlap, and the overlapping parts of \hat f(\omega - k\omega_s) for different k add together in a way that cannot be undone by any filter — the periodized spectrum no longer determines the original \hat f uniquely, so no reconstruction formula, sinc-based or otherwise, can recover f from those samples in general.

The proof is by induction on the order M, using one key fact about how a single lattice stage relates an order-m polynomial to its order-(m-1) reduction.

Base case (M=1). A_1(z) = 1 + k_1z^{-1} has its single root at z=-k_1, which lies inside the unit circle exactly when |k_1|<1. The claim holds trivially.

Inductive step. Assume the claim holds for every order less than m. Define the reciprocal polynomial A_m^*(z) = z^{-m}A_m(z^{-1}) (the coefficients of A_m read in reverse). A direct computation from the step-down formula shows A_m(z) = A_{m-1}(z) + k_mz^{-1}A_{m-1}^*(z), so A_m's roots are controlled jointly by A_{m-1}'s roots and k_m. The Schur-Cohn criterion for this recursive polynomial family (a standard result relating a polynomial's root locations to a scaled combination of a lower-order polynomial and its own reciprocal) states that A_m(z) has every root strictly inside the unit circle if and only if |k_m|<1 and A_{m-1}(z) (the very polynomial the step-down recursion produces from A_m) itself has every root strictly inside the unit circle. By the inductive hypothesis, A_{m-1}(z)'s roots are all inside the unit circle exactly when |k_1|,\dots,|k_{m-1}|<1. Combining the two conditions gives exactly |k_1|,\dots,|k_m|<1, completing the induction.

Applying this chain of equivalences down from M to 1 gives the full theorem: A_M(z)'s roots are all strictly inside the unit circle if and only if every reflection coefficient k_1,\dots,k_M has magnitude strictly less than 1.

A stable, causal system has an impulse response h[n] satisfying h[n]=0 for n<0 and \sum_n |h[n]| < \infty. Consider the causal z-transform H(z) = \sum_{n=0}^{\infty} h[n]z^{-n}; stability places the unit circle inside its region of convergence, and causality means this series has no negative powers of z at all — H(z) is analytic (has a convergent power series) throughout |z|>r for some r<1, in particular on and outside the unit circle.

For H(e^{j\omega}) to vanish identically on a set of nonzero measure on the unit circle while H(z) remains analytic in the region containing that circle would force H(z)\equiv 0 everywhere, by the identity theorem for analytic functions (an analytic function that is zero on a set with an accumulation point is zero on its entire connected domain). A system function that is identically zero is not a meaningful filter, so any genuine causal, stable H(z) that isn't trivially zero must have |H(e^{j\omega})|>0 almost everywhere.

The integral condition sharpens this qualitative fact into a precise quantitative one: writing H(z)'s log-magnitude in terms of its (causal) cepstrum coefficients via the discrete-time Fourier series of \ln H(e^{j\omega}), causality of h[n] forces the corresponding one-sided causal structure on the cepstrum, and this structure is only consistent with a function whose log-magnitude is integrable in this sense — a magnitude response can decay steeply near a band edge, but an EXACT, band-limited zero (an infinitely deep null over a whole sub-interval) makes \ln|H(e^{j\omega})|=-\infty there, making the integral diverge and violating the condition. No causal, stable system can therefore realize a magnitude response with a truly flat zero over any nonzero-width band — the mathematical reason an ideal brick-wall filter is impossible, not merely impractical.

Sufficiency and optimality. Suppose E(\omega) alternates L+2 times as described, but some other degree-L trigonometric polynomial A'(\omega) achieves a strictly smaller maximum weighted error \delta' < \delta. Consider the difference Q(\omega) = A(\omega) - A'(\omega), itself a trigonometric polynomial of degree at most L. At each alternation point \omega_i, E(\omega_i)=\pm\delta while the corresponding error for A' has magnitude at most \delta'<\delta, so Q(\omega_i)=W(\omega_i)^{-1}\big(E(\omega_i) - E'(\omega_i)\big) must share the same sign as E(\omega_i) at every one of the L+2 alternation points (the smaller \delta' cannot cancel the full swing of \pm\delta). Since the sign of E(\omega_i) itself alternates across the L+2 points, Q(\omega) changes sign at least L+1 times between them, forcing at least L+1 zeros of Q — but a nonzero trigonometric polynomial of degree at most L can have at most 2L zeros in one period, and more specifically an argument bounding zeros consistent with L+1 sign changes across L+2 ordered points is only possible if Q\equiv 0, contradicting that A'\ne A achieves a genuinely smaller error. So no strictly better approximation exists, and A(\omega) is optimal; the same zero-counting argument applied to any two distinct alternating-optimal solutions shows they must coincide, giving uniqueness.

Necessity. Conversely, if A(\omega) is the true minimax solution, standard results from Chebyshev approximation theory (the same machinery underlying ordinary polynomial minimax approximation, adapted here to trigonometric polynomials) guarantee that fewer than L+2 alternations would leave room to perturb A(\omega) slightly — in the direction opposite the error's sign at every existing extremum — and strictly reduce the maximum error everywhere, contradicting optimality. So the optimal solution must alternate at least L+2 times.

Property 1 (decorrelation). By the spectral theorem for real symmetric matrices, R_x=\Phi\Lambda\Phi^T for the diagonal eigenvalue matrix \Lambda=\text{diag}(\lambda_1,\dots,\lambda_N) and orthogonal \Phi (its columns the orthonormal eigenvectors). The transform coefficients' own covariance is R_y = E[YY^T] = E[\Phi^TXX^T\Phi] = \Phi^TR_x\Phi = \Phi^T(\Phi\Lambda\Phi^T)\Phi = \Lambda, using \Phi^T\Phi=I (orthogonality) twice. R_y=\Lambda is diagonal, so \text{Cov}(Y_i,Y_j)=\lambda_i\delta_{ij} exactly — the KLT coefficients are uncorrelated, with Y_i's own variance equal to \lambda_i.

Property 2 (optimal truncation), the case M=N-1. Reconstructing X from an orthogonal transform Y'=\Psi^TX while discarding one coefficient Y'_k (replacing it with its mean, here 0) introduces mean-square error exactly E[(Y'_k)^2] = \psi_k^TR_x\psi_k, the Rayleigh quotient of R_x at the discarded direction \psi_k. Since R_x's eigenvalues bound its Rayleigh quotient, \lambda_N \le \psi_k^TR_x\psi_k \le \lambda_1 for any unit vector \psi_k, the discarded direction minimizing this error is exactly \phi_N (the eigenvector of the smallest eigenvalue) — equivalently, keeping the N-1 coefficients of largest eigenvalue and discarding the smallest is optimal, with minimum achievable error \lambda_N.

General M. Applying this one-coefficient argument repeatedly — at each step, among the directions not yet fixed, discarding the one with the smallest remaining Rayleigh quotient is optimal, and by the same bounding argument this is always the next-smallest eigenvector — shows by induction that discarding the N-M smallest-eigenvalue directions (equivalently, keeping the M largest) minimizes the total mean-square truncation error among every orthogonal transform, with minimum achievable error \sum_{i=M+1}^{N}\lambda_i.

The exponential kernel factors as a product of a purely-m term and a purely-n term: e^{-i2\pi(km/M + ln/N)} = e^{-i2\pi km/M}\cdot e^{-i2\pi ln/N}. Substituting this into the defining double sum and pulling the m-independent factor e^{-i2\pi km/M} outside the inner sum over n (a finite sum, so reordering is unconditionally valid): X_{k,l} = \sum_{m=0}^{M-1} e^{-i2\pi km/M} \sum_{n=0}^{N-1} x_{m,n}\,e^{-i2\pi ln/N} = \sum_{m=0}^{M-1} e^{-i2\pi km/M}\,Y_{m,l}, which is exactly the stated two-stage computation: the inner sum over n is a length-N DFT applied to row m (producing Y_{m,l}), and the outer sum over m is a length-M DFT applied down column l of Y. Symmetric factoring — pulling out the n-independent factor first instead — gives the columns-then-rows order with the identical final sum X_{k,l}, confirming the order doesn't matter.

Consequence for computation. Rather than evaluating the double sum directly for every (k,l) pair, at O(M^2N^2) total cost, the 2-D DFT can be computed as M independent length-N 1-D DFTs (one per row) followed by N independent length-M 1-D DFTs (one per column). Using the FFT for each 1-D transform, this costs O(MN\log N) for the row pass plus O(MN\log M) for the column pass, i.e. O(MN\log(MN)) overall — the same divide-and-conquer speedup the 1-D FFT gets, applied once per dimension.