Discrete Cosine Transform

The discrete cosine transform (DCT) rewrites a finite sequence of samples, such as a row of pixels or an $8 \times 8$ block, as a weighted sum of sampled cosines of increasing frequency. It is real, orthogonal when suitably scaled, and computed quickly with the FFT. Its importance in imaging comes from one property: for smooth, strongly correlated blocks, most of the energy lands in a few low-frequency coefficients. The JPEG standard transforms every $8 \times 8$ block with it before quantizing (CCITT, 1992). Compared with the Fourier transform, which treats a block as one period of a periodic signal, the DCT implicitly extends the block by mirror reflection, so no artificial jump appears at its borders.

There are several DCTs, which differ in how the sequence is assumed to continue beyond its ends. “The” DCT usually means the type II transform, and “the” inverse DCT means type III (SciPy 1.16.3 documentation). Strang (1999) dates the DCT to a 1974 paper by Ahmed, Natarajan and Rao, and describes that paper as deriving the DCT-II basis as an approximation to the eigenvectors of the covariance matrix $\rho^{|j-k|}$ of strongly correlated signals, for which the exact eigenvectors form the optimal Karhunen–Loève basis.

Definition

The DCT of a sequence $x_0, \dots, x_{N-1}$ is the sequence of its inner products with $N$ vectors whose components are cosines sampled at $N$ evenly spaced points. The different types choose whether the samples sit on the grid points or half-way between them, and whether the frequencies are whole or half-integer multiples of a base frequency; each choice is equivalent to a boundary condition at the two ends of the sequence. With the right scale factors the cosine vectors are orthonormal, so the transform is a rotation of $\mathbb{R}^N$: it loses nothing and preserves energy.

Intuition

A DFT basis vector repeats every $N$ samples, so the DFT sees a block as one tile of a periodic pattern. A row of pixels rarely ends where it began, and the periodic continuation then has a jump at every tile boundary. Strang (1999, section 7) points out the consequence: a jump makes the Fourier coefficients decay slowly and produces Gibbs oscillation, seen in images as ringing. Reflecting the block across its boundary removes the jump, since the reflected copy starts where the original ends, and a periodic signal that is continuous needs far fewer high frequencies. The cosines of the DCT-II are exactly the basis of such mirror-extended sequences: the $k$-th vector is a cosine that completes $k/2$ cycles across the block, so it is symmetric about both block edges.

The first vector is flat, a constant grey level, and its coefficient is the block’s mean up to a scale factor (the DC coefficient, in JPEG’s terms). The next vectors are a half-cycle ramp, a full cycle, and so on. A smooth block is mostly its mean plus a gentle slope, so most of its energy falls in the first few coefficients and the rest are small. This is energy compaction, and it is what makes the DCT useful for compression and denoising: small coefficients can be coarsely quantized or discarded with little damage.

Formal Definition

The DCT-II and its inverse

Strang (1999, section 2) writes the $k$-th DCT-II basis vector as having $j$-th component $\cos\big((j + \tfrac12) k \tfrac{\pi}{N}\big)$, with the $k = 0$ vector divided by $\sqrt2$; all these vectors have length $\sqrt{N/2}$. Scaling them to unit length gives the orthonormal DCT-II,

$$X_k = \alpha_k \sum_{n=0}^{N-1} x_n \cos\!\left( \frac{\pi (2n+1) k}{2N} \right), \qquad \alpha_0 = \sqrt{\tfrac{1}{N}}, \quad \alpha_k = \sqrt{\tfrac{2}{N}} \;\; (k \ge 1),$$

for $k = 0, \dots, N-1$. Written as a matrix, $X = Cx$ with $C_{kn} = \alpha_k \cos(\pi(2n+1)k/2N)$, and $C$ is orthogonal, so the inverse is $x = C^\top X$:

$$x_n = \sum_{k=0}^{N-1} \alpha_k X_k \cos\!\left( \frac{\pi (2n+1) k}{2N} \right).$$

Strang notes that the DCT-III matrix is the transpose of the DCT-II matrix, which is why the inverse DCT-II is a DCT-III.

Two dimensions and JPEG’s $8 \times 8$ form

The two-dimensional DCT applies the one-dimensional transform to every row and then to every column: it is the tensor product of two 1D transforms (Strang, 1999, section 7). For an $N \times N$ block $x_{mn}$,

$$X_{uv} = \alpha_u \alpha_v \sum_{m=0}^{N-1} \sum_{n=0}^{N-1} x_{mn} \cos\!\left( \frac{\pi (2m+1) u}{2N} \right) \cos\!\left( \frac{\pi (2n+1) v}{2N} \right),$$

and its basis images are products of a vertical and a horizontal cosine. If the block is flattened row by row into a vector of $N^2$ values, the 2D transform matrix is the Kronecker product $C \otimes C$.

The JPEG standard groups each image component into $8 \times 8$ blocks, subtracts $2^{P-1}$ from every sample of precision $P$ bits (128 for 8-bit images: the level shift), and transforms each block into 64 coefficients, one DC and 63 AC coefficients (CCITT, 1992, sections 4.3 and A.3.1). Its section A.3.3 gives the forward DCT as

$$S_{vu} = \tfrac14 C_u C_v \sum_{x=0}^{7} \sum_{y=0}^{7} s_{yx} \cos\frac{(2x+1)u\pi}{16} \cos\frac{(2y+1)v\pi}{16}, \qquad C_u, C_v = \tfrac{1}{\sqrt2} \text{ for } u, v = 0, \text{ and } 1 \text{ otherwise},$$

where $x$ and $u$ run horizontally and $y$ and $v$ vertically. With $N = 8$, $\alpha_u \alpha_v = \tfrac14 C_u C_v$, so this is the orthonormal 2D DCT-II above; the first example checks it numerically. The standard calls these equations an informative “ideal functional definition” and leaves the accuracy that real implementations must reach to Part 2 of the specification. Everything JPEG does after the transform, quantizing each coefficient with a table, scanning in zig-zag order and entropy coding, belongs to image compression, not to the transform; the standard specifies no default quantization tables (CCITT, 1992, section 4.3).

Normalizations in software

Libraries differ in scale factors, and the differences matter when coefficients are compared across tools. SciPy 1.16.3 implements types I to IV (its documentation notes that eight types exist in theory). With the default norm="backward", its DCT-II is unnormalized,

$$y_k = 2 \sum_{n=0}^{N-1} x_n \cos\!\left(\frac{\pi k (2n+1)}{2N}\right),$$

and the inverse carries the scaling. With norm="ortho", both directions are scaled by the same factor and, by default (orthogonalize=True), $y_0$ is divided by $\sqrt2$, which makes the transform matrix orthonormal: the $C$ of the formula above. The documentation adds that dct(x, norm="ortho") equals MATLAB’s dct(x), and warns that for types I to III the orthonormal form no longer corresponds directly to the DFT unless orthogonalize=False is passed (scipy.fft.dct docstring, SciPy 1.16.3). dctn and idctn apply the transform along several axes.

The family of types

Strang (1999) derives every DCT from one observation: the cosine vectors are the eigenvectors of symmetric second-difference matrices, which have rows $-1, 2, -1$ inside and differ only in their first and last rows. Those boundary rows encode the boundary condition. At each end the discrete analogue of $u'(0) = 0$ can be centred on a grid point (whole-sample symmetry, $u_{-1} = u_1$) or on a midpoint (half-sample symmetry, $u_{-1} = u_0$), and at the right end a zero-value (Dirichlet) condition is also possible (sections 3 and 4). With the same centring at both ends this gives the four standard types:

Type	Centres (Strang’s term)	Component $j$ of vector $k$	Boundary rows of $A$
DCT-I	$j = 0$ and $N-1$	$\cos\big(jk\frac{\pi}{N-1}\big)$	$(2, -2)$ and $(-2, 2)$
DCT-II	$j = -\frac12$ and $N - \frac12$	$\cos\big((j+\frac12)k\frac{\pi}{N}\big)$	$(1, -1)$ and $(-1, 1)$
DCT-III	$j = 0$ and $N$	$\cos\big(j(k+\frac12)\frac{\pi}{N}\big)$	$(2, -2)$ and $(-1, 2)$
DCT-IV	$j = -\frac12$ and $N - \frac12$	$\cos\big((j+\frac12)(k+\frac12)\frac{\pi}{N}\big)$	$(1, -1)$ and $(-1, 3)$

The eigenvalues are $2 - 2\cos\theta$ with $\theta = k\pi/(N-1)$, $k\pi/N$, $(k + \frac12)\pi/N$ and $(k + \frac12)\pi/N$ respectively. Types I and III need a $\sqrt2$ rescaling of end components to make their bases orthogonal. The final 3 in the DCT-IV matrix comes from a zero value imposed at a midpoint, which makes the extension antisymmetric at the right end: the DCT-IV vectors are even at the left end and odd at the right (Strang, 1999, section 4). The other four combinations, with one end centred on a grid point and the other on a midpoint, give DCT-5 to DCT-8, whose cosines have the half-integer periods $N - \frac12$ or $N + \frac12$; Strang (section 5) writes that these “odd” transforms are not used in signal processing; his introduction credits the complete set of eight to Wang in 1985, citing a paper by Wang and Hunt. Replacing the zero-slope condition at the left end by a zero value gives eight discrete sine transforms in the same way.

Properties

Orthogonality and energy. The orthonormal DCT of every type is an orthogonal matrix, so $\lVert X \rVert = \lVert x \rVert$ and errors keep their size between the pixel and coefficient domains. Mean Squared Error lists this among the MSE’s properties: the MSE between two images equals the MSE between their orthonormal DCTs.

Relation to the DFT. Mirror-extend $x$ to length $2N$, $\tilde x = (x_0, \dots, x_{N-1}, x_{N-1}, \dots, x_0)$, the half-sample symmetric extension. Its DFT is

$$\tilde X_k = \sum_{n=0}^{N-1} x_n \left( e^{-i\pi k n/N} + e^{i\pi k (n+1)/N} \right) = 2\, e^{i\pi k/2N} \sum_{n=0}^{N-1} x_n \cos\!\left(\frac{\pi k (2n+1)}{2N}\right),$$

so SciPy’s unnormalized DCT-II is $y_k = e^{-i\pi k/2N} \tilde X_k$ for $k < N$ (our derivation; the first example checks it). The DCT-II is therefore a DFT of twice the length on data whose periodic extension has no jump, which makes precise the Fourier Transform article’s remark that the DCT extends a block by mirror reflection rather than tiling. The DCT is not the real part of an $N$-point DFT: the DFT’s cosines $\cos(2\pi jk/N)$ complete whole cycles, while DCT-II index $k$ is $k$ half-cycles.

Fast computation. It follows from this link that an FFT computes the DCT in $O(N \log N)$ operations instead of the $N^2$ of a matrix product. Strang (1999) notes that the DCT-II and DCT-IV have FFT implementations, and reports a factorization due to Wickerhauser that obtains the DCT-IV and DST-IV together from one DFT of order $2N$. SciPy computes its DCTs with the pocketfft library.

Separability. The 2D transform is a 1D transform of the rows followed by one of the columns, so an $N \times N$ block costs $2N$ 1D transforms.

Energy compaction and the Karhunen–Loève transform. For a random vector with covariance $\Sigma$, the orthonormal transform that leaves the coefficients uncorrelated is the matrix of eigenvectors of $\Sigma$, the Karhunen–Loève transform (KLT), which is principal component analysis applied to a known covariance. Among all orthonormal transforms it puts the largest possible share of the expected energy into the first $k$ coefficients, for every $k$ (the sum of the $k$ largest diagonal entries of $T\Sigma T^\top$ cannot exceed the sum of the $k$ largest eigenvalues of $\Sigma$). In our reading, its drawback is that it depends on $\Sigma$, which must be known or estimated, whereas the DCT is fixed.

The link between the two is the first-order autoregressive, or AR(1), model, whose correlation between samples $|j - k|$ apart is $\rho^{|j-k|}$. Strang (1999, section 2) describes this matrix as the covariance of a useful class of signals, with $\rho$ near 1 measuring the correlation between neighbours, and says that the DCT vectors are close to the optimal Karhunen–Loève basis for it while being independent of $\rho$. The reason can be read off Strang’s construction (our derivation). The inverse of the AR(1) correlation matrix $R$ is tridiagonal:

$$(1 - \rho^2)\, R^{-1} = \rho A_2 + (1 - \rho) E, \qquad E = \mathrm{diag}(1, 1-\rho, \dots, 1-\rho, 1),$$

where $A_2$ is the second-difference matrix with boundary rows $(1, -1)$ and $(-1, 1)$ whose eigenvectors are the DCT-II vectors. $R$ and $R^{-1}$ have the same eigenvectors, and as $\rho \to 1$ the correction $(1 - \rho)E/\rho$ vanishes while the eigenvalues $2 - 2\cos(k\pi/N)$ of $A_2$ stay distinct, so the KLT vectors tend to the DCT-II vectors, with angles of first order in $(1 - \rho)/\rho$. The third example measures this convergence.

To summarize compaction by one number, the examples use the ratio of the arithmetic to the geometric mean of the coefficient variances $\sigma_i^2$,

$$G = \frac{\frac1n \sum_i \sigma_i^2}{\left( \prod_i \sigma_i^2 \right)^{1/n}},$$

reported in decibels, $10 \log_{10} G$. Strang (1999, section 5) uses the term “coding gain” for comparing transforms, without defining it in the preprint; $G$ is the ratio we use under that name. It is 1 (0 dB) when all coefficients have the same variance, as for the identity transform on a stationary signal. Among orthonormal transforms the arithmetic mean is fixed, since it is $\mathrm{tr}(\Sigma)/n$, and Hadamard’s inequality bounds the product of the diagonal entries of $T\Sigma T^\top$ from below by $\det\Sigma$, with equality only when $T\Sigma T^\top$ is diagonal; so the KLT has the largest $G$ (our reasoning). For a separable 2D covariance $R \otimes R$ and a separable transform $T \otimes T$, the coefficient variances are all products of the 1D variances, so the 2D ratio is the square of the 1D ratio and its value in decibels doubles.

Symmetric convolution. Circular convolution is diagonalized by the DFT. Strang (1999, section 6) reports Martucci’s answer to the corresponding question for cosine transforms: extend the filter and the signal symmetrically, convolve them cyclically, and keep the original $N$ outputs; this symmetric convolution corresponds to multiplying the cosine series. Since each end can be extended about a sample or a midpoint, the filter’s extension can differ from the signal’s.

Blocks and their borders. Mirror extension removes the seam inside each block’s transform, but a block transform does nothing for the borders between blocks. When the coefficients of neighbouring blocks are quantized independently, the reconstructed blocks no longer join smoothly: Strang (1999, section 7) writes that JPEG’s blocking becomes severe as the compression ratio grows (the preprint does not include the reconstructions his text refers to). The remedy he describes is overlap between blocks. In Malvar’s modulated lapped transform, as Strang describes it, each of the $N$ basis vectors has length $2N$ and overlaps the next block: its $j$-th component is $p(j)\cos\big[(k + \frac12)(j + \frac{N+1}{2})\frac{\pi}{N}\big]$, a window $p$ modulated by the DCT-IV, and the transform stays orthogonal provided $p^2(j) + p^2(j + N) = 1$.

Examples

The script below checks the definitions, then measures energy compaction on a synthetic image with a known covariance, then follows the KLT as the correlation approaches 1. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 (for one $8 \times 8$ block of the bundled camera image), and runs in under a second on one thread. Nothing is quantized: quantization belongs to compression codecs. Arrays are indexed row first; blocks are flattened row by row, so the 2D transforms are Kronecker products of 1D matrices whose rows are the basis vectors.

The first part checks the conventions: SciPy’s orthonormal DCT-II against Strang’s cosines, the orthogonality of all four types, the four second-difference matrices of the table above, the DFT of the mirrored row, JPEG’s formula written out as a double sum, and the tridiagonal form of the AR(1) inverse.

import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"  # single-threaded BLAS, set before NumPy is imported
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy.fft import dct, dctn, fft
from scipy.signal import lfilter
from skimage import data

N = 8
j = np.arange(N)

# ---------- Part A: definitions and identities ----------
C = dct(np.eye(N), norm="ortho", axis=0)          # row k = k-th orthonormal DCT-II basis vector
strang = np.array([np.cos((j + 0.5) * k * np.pi / N) for k in range(N)])
strang[0] /= np.sqrt(2)
print(f"A1 SciPy ortho DCT-II = cos((j+1/2)k pi/N) basis:   {np.abs(C - strang * np.sqrt(2 / N)).max():.1e}")
for t in (1, 2, 3, 4):
    M = dct(np.eye(N), type=t, norm="ortho", axis=0)
    print(f"A2 type {t}: max |M M^T - I| = {np.abs(M @ M.T - np.eye(N)).max():.1e}")

def second_difference(top, bottom):               # rows -1, 2, -1 with the given boundary rows
    A = 2 * np.eye(N) - np.eye(N, k=1) - np.eye(N, k=-1)
    A[0, :2], A[-1, -2:] = top, bottom
    return A

bc = {1: ((2, -2), (-2, 2), lambda k: k * np.pi / (N - 1), lambda th: np.cos(j * th)),
      2: ((1, -1), (-1, 1), lambda k: k * np.pi / N, lambda th: np.cos((j + 0.5) * th)),
      3: ((2, -2), (-1, 2), lambda k: (k + 0.5) * np.pi / N, lambda th: np.cos(j * th)),
      4: ((1, -1), (-1, 3), lambda k: (k + 0.5) * np.pi / N, lambda th: np.cos((j + 0.5) * th))}
for t, (top, bottom, theta, vec) in bc.items():
    A = second_difference(top, bottom)
    res = max(np.abs(A @ vec(theta(k)) - (2 - 2 * np.cos(theta(k))) * vec(theta(k))).max()
              for k in range(N))
    print(f"A3 DCT-{t} cosines are eigenvectors of A{t}: residual {res:.1e}")

block = data.camera().astype(float)[200:208, 240:248] - 128   # level shift, T.81 A.3.1
x = block[0]
mirror = np.real(np.exp(-1j * np.pi * j / (2 * N)) * fft(np.concatenate([x, x[::-1]]))[:N])
print(f"A4 DCT-II = phase-shifted DFT of the mirrored row:   {np.abs(dct(x) - mirror).max():.1e}")
S = np.zeros((N, N))
for v in range(N):
    for u in range(N):
        cu, cv = (1 / np.sqrt(2) if u == 0 else 1), (1 / np.sqrt(2) if v == 0 else 1)
        S[v, u] = 0.25 * cu * cv * np.sum(block * np.cos((2 * j[None, :] + 1) * u * np.pi / 16)
                                          * np.cos((2 * j[:, None] + 1) * v * np.pi / 16))
print(f"A5 T.81 FDCT formula = dctn(norm='ortho'):          {np.abs(S - dctn(block, norm='ortho')).max():.1e}")

def ar1_corr(rho):
    return rho ** np.abs(j[:, None] - j[None, :])

for rho in (0.9, 0.99, 0.999):
    E = np.diag([1] + [1 - rho] * (N - 2) + [1])
    lhs = (1 - rho ** 2) * np.linalg.inv(ar1_corr(rho))
    print(f"A6 rho={rho}: (1-rho^2) R^-1 = rho A2 + (1-rho) E: "
          f"{np.abs(lhs - (rho * second_difference((1, -1), (-1, 1)) + (1 - rho) * E)).max():.1e}")

A1 SciPy ortho DCT-II = cos((j+1/2)k pi/N) basis:   4.4e-16
A2 type 1: max |M M^T - I| = 2.2e-16
A2 type 2: max |M M^T - I| = 2.2e-16
A2 type 3: max |M M^T - I| = 2.2e-16
A2 type 4: max |M M^T - I| = 2.2e-16
A3 DCT-1 cosines are eigenvectors of A1: residual 8.9e-16
A3 DCT-2 cosines are eigenvectors of A2: residual 1.8e-15
A3 DCT-3 cosines are eigenvectors of A3: residual 3.6e-15
A3 DCT-4 cosines are eigenvectors of A4: residual 2.9e-15
A4 DCT-II = phase-shifted DFT of the mirrored row:   8.9e-16
A5 T.81 FDCT formula = dctn(norm='ortho'):          4.3e-14
A6 rho=0.9: (1-rho^2) R^-1 = rho A2 + (1-rho) E: 2.2e-15
A6 rho=0.99: (1-rho^2) R^-1 = rho A2 + (1-rho) E: 4.6e-14
A6 rho=0.999: (1-rho^2) R^-1 = rho A2 + (1-rho) E: 1.5e-13

Every identity holds to rounding error. The SciPy matrix is Strang’s basis rescaled to unit length, and all four orthonormal types are orthogonal matrices; the four second-difference matrices have the pure cosines of their types as eigenvectors, so the types really are distinguished by their boundary rows alone; and the unnormalized DCT-II of a row equals the DFT of its mirrored, double-length copy after the half-sample phase shift. JPEG’s FDCT on a level-shifted camera block agrees with dctn(..., norm="ortho") to rounding error, so a JPEG coefficient and a SciPy orthonormal coefficient are the same number. The last three lines confirm the identity behind the AR(1) argument at three correlations. All six planned checks held, against a tolerance of $10^{-12}$ for A1 to A4 and $10^{-10}$ for A5 and A6.

The second part builds a $1024 \times 1024$ Gaussian field whose correlation is exactly $\rho^{|\Delta r|}\rho^{|\Delta c|}$ with $\rho = 0.95$: white noise is filtered by the stationary AR(1) recursion down the columns and then along the rows. Any $8 \times 8$ block of it therefore has the covariance $\Sigma = R \otimes R$, and the exact KLT is the Kronecker product of the eigenvectors of $R$. The script compares the KLT, the DCT, the DFT and the identity (raw pixels), both exactly from $\Sigma$ and empirically over the 16,384 non-overlapping blocks. The DFT needs a choice the sources leave open: an $N$-point real signal has complex coefficients that come in conjugate pairs. “DFT real” uses Strang’s real basis, the normalized cosine and sine vectors $\cos(2\pi jk/N)$ and $\sin(2\pi jk/N)$ (section 1); “DFT complex” uses the unitary complex DFT and the variances $E|X|^2$ of its 64 coefficients, which counts each conjugate pair as two coefficients. The script prints both.

# ---------- Part B: energy compaction on a separable AR(1) field ----------
def klt(rho):                                     # exact KLT: eigenvectors of R, largest variance first
    w, V = np.linalg.eigh(ar1_corr(rho))
    return V[:, np.argsort(w)[::-1]].T

def real_dft():                                   # Strang's cos and sin vectors, normalized
    rows = []
    for k in range(N // 2 + 1):
        rows.append(np.cos(2 * np.pi * j * k / N))
        if 0 < k < N // 2:
            rows.append(np.sin(2 * np.pi * j * k / N))
    return np.array([r / np.linalg.norm(r) for r in rows])

def coef_var(T, cov):                             # variances of the coefficients T x
    return np.real(np.einsum("ij,jk,ik->i", T, cov, T.conj()))

def energy(var, ks=(1, 2, 4, 8, 16, 32)):         # share of energy in the k largest-variance coefficients
    c = np.cumsum(np.sort(var)[::-1]) / var.sum()
    return np.array([c[k - 1] for k in ks])

def gain_db(var):                                 # arithmetic over geometric mean of the variances
    return 10 * np.log10(var.mean() / np.exp(np.mean(np.log(var))))

rho = 0.95
rng = np.random.default_rng(0)
w = rng.standard_normal((1024, 1024))
s = np.sqrt(1 - rho ** 2)
for axis in (0, 1):                               # y0 = w0, y_n = rho y_(n-1) + s w_n, along each axis
    w = np.moveaxis(w, axis, 0).copy()
    w[0] /= s
    w = np.moveaxis(lfilter([s], [1, -rho], w, axis=0), 0, axis)
field = w
print(f"\nB0 field: variance {field.var():.4f}, lag-1 corr rows {np.mean(field[:, 1:] * field[:, :-1]) / np.mean(field ** 2):.4f}"
      f" cols {np.mean(field[1:] * field[:-1]) / np.mean(field ** 2):.4f}, diagonal"
      f" {np.mean(field[1:, 1:] * field[:-1, :-1]) / np.mean(field ** 2):.4f} (rho^2 = {rho ** 2:.4f})")
blocks = field.reshape(128, N, 128, N).transpose(0, 2, 1, 3).reshape(-1, N * N)   # 16,384 blocks

R = ar1_corr(rho)
Sigma = np.kron(R, R)                             # exact covariance of a row-major 8 x 8 block
F1 = fft(np.eye(N), axis=0, norm="ortho")         # unitary DFT matrix
one_d = {"KLT": klt(rho), "DCT": C, "DFT real": real_dft(), "DFT complex": F1, "identity": np.eye(N)}
T = {name: np.kron(M, M) for name, M in one_d.items()}

P = T['KLT'] @ Sigma @ T['KLT'].T
print(f"B1 KLT coefficients uncorrelated: max |off-diagonal| {np.abs(P - np.diag(np.diag(P))).max():.1e}")
print("\nexact          k=1      k=2      k=4      k=8     k=16     k=32   gain dB  empirical gain dB")
ex_E, ex_G, emp_G, dE = {}, {}, {}, 0.0
for name, M in T.items():
    var = coef_var(M, Sigma)
    emp = np.mean(np.abs(blocks @ M.T) ** 2, axis=0)
    ex_E[name], ex_G[name], emp_G[name] = energy(var), gain_db(var), gain_db(emp)
    dE = max(dE, np.abs(energy(emp) - ex_E[name]).max())
    print(f"{name:12s}" + "".join(f" {e:8.5f}" for e in ex_E[name]) + f" {ex_G[name]:8.4f}  {emp_G[name]:8.4f}")
print(f"B2 max |empirical - exact| energy share {dE:.5f}, gain "
      f"{max(abs(emp_G[n] - ex_G[n]) for n in T):.4f} dB")
print(f"B3 KLT - DCT gain gap: exact {ex_G['KLT'] - ex_G['DCT']:.5f} dB, empirical (same blocks) "
      f"{emp_G['KLT'] - emp_G['DCT']:.5f} dB")
print(f"B4 G2D - 2 G1D, max over transforms: "
      f"{max(abs(ex_G[n] - 2 * gain_db(coef_var(M, R))) for n, M in one_d.items()):.1e} dB")
ge = lambda a, b: bool(np.all(ex_E[a] >= ex_E[b] - 1e-12))
print("B5 energy KLT >= DCT >= both DFTs >= identity at every k:",
      ge("KLT", "DCT") and ge("DCT", "DFT real") and ge("DCT", "DFT complex")
      and ge("DFT real", "identity") and ge("DFT complex", "identity"))

B0 field: variance 0.9626, lag-1 corr rows 0.9484 cols 0.9482, diagonal 0.8992 (rho^2 = 0.9025)
B1 KLT coefficients uncorrelated: max |off-diagonal| 1.2e-14

exact          k=1      k=2      k=4      k=8     k=16     k=32   gain dB  empirical gain dB
KLT           0.77227  0.83544  0.91710  0.95914  0.99048  0.99857  17.6924   17.5269
DCT           0.77109  0.83419  0.91632  0.95896  0.99031  0.99855  17.6518   17.4883
DFT real      0.77109  0.81691  0.88636  0.94067  0.98787  0.99688  15.4079   15.2396
DFT complex   0.77109  0.80581  0.87526  0.94067  0.98671  0.99644  15.1747   15.0050
identity      0.01562  0.03125  0.06250  0.12500  0.25000  0.50000   0.0000    0.0002
B2 max |empirical - exact| energy share 0.00964, gain 0.1697 dB
B3 KLT - DCT gain gap: exact 0.04060 dB, empirical (same blocks) 0.03857 dB
B4 G2D - 2 G1D, max over transforms: 1.6e-13 dB
B5 energy KLT >= DCT >= both DFTs >= identity at every k: True

The field has the planned statistics: unit variance up to sampling error, neighbour correlations of about 0.95 along both axes and about $0.95^2$ diagonally, as a separable model requires; the three planned data checks held (variance within 0.1 of 1, each correlation within 0.01 of its target). The KLT coefficients are uncorrelated to rounding error, as they must be, and the B1 check held against its tolerance of $10^{-12}$. Before the first run, with the editor’s approval, B1 was changed to read the off-diagonal entries of a single matrix; the piloted version subtracted a second computation of the diagonal, so its printed value included rounding differences between the two. The exact table shows the ordering the theory predicts, and the B5 line confirms it at every $k$: the KLT captures the most energy, the DCT trails it by at most 0.0013 at every $k$, and both DFT accountings fall behind from $k = 2$ on. At $k = 1$ the DCT and both DFTs tie, because all three share the flat DC vector. In coding gain the KLT exceeds the DCT by 0.041 dB, while the DCT beats the DFT by about 2.2 dB with the real basis and 2.5 dB with the complex accounting. The DFT’s loss is the cost of the periodic seam: its basis is built for a block whose ends meet, and an AR(1) block’s ends do not. The three planned hypotheses held: the ordering of energies at every $k$ and of gains, a KLT–DCT gain gap below 0.1 dB, and a DCT–DFT gap above 1 dB under both accountings. B4 confirms that the 2D gains are exactly twice the 1D gains in decibels for every transform; that check held against its tolerance of $10^{-9}$ dB.

The empirical columns show why the comparison uses the exact covariance. Sampling error on one field moves each empirical energy share by less than 0.01 and each empirical gain by up to 0.17 dB from the exact value, four times the KLT–DCT difference; in our reading this is because the long-range correlation leaves neighbouring blocks far from independent. The error is shared, though: on the same blocks, the empirical KLT–DCT gap is 0.039 dB, against 0.041 dB exactly. The three planned sampling checks held: energy shares within 0.015 and gains within 0.25 dB of their exact values, and the paired gap within 0.005 dB of the exact gap. These tolerances were fixed at planning from a run on the same field; the planning value behind the gain tolerance left out the complex DFT, whose row has the largest gain error in the printed table.

The third part works in one dimension with $N = 8$ and sweeps $\rho$. For each value it checks that the KLT vector with the $k$-th largest eigenvalue is the one closest to DCT vector $k$ (“matched”), then prints the largest angle between corresponding vectors and the 1D gains.

# ---------- Part C: the DCT basis approaches the KLT as rho -> 1 ----------
print("\nrho     max angle deg  KLT gain  DCT gain  KLT-DCT gap  DFT real  DFT complex  matched")
angles = {}
for r in (0.5, 0.8, 0.9, 0.95, 0.99, 0.999):
    V, Rr = klt(r), ar1_corr(r)
    matched = bool(np.all(np.argmax(np.abs(C @ V.T), axis=0) == j))   # KLT vector k pairs with DCT vector k
    angles[r] = np.degrees(np.arccos(np.clip(np.abs(np.sum(C * V, axis=1)), 0, 1))).max()
    g = [gain_db(coef_var(M, Rr)) for M in (V, C, real_dft(), F1)]
    print(f"{r:<7} {angles[r]:12.4g} {g[0]:9.4f} {g[1]:9.4f} {g[0] - g[1]:12.5f} {g[2]:9.4f} {g[3]:12.4f}  {matched}")
print(f"C1 angle(0.99) / angle(0.999) = {angles[0.99] / angles[0.999]:.4f}; "
      f"first-order prediction {((1 - 0.99) / 0.99) / ((1 - 0.999) / 0.999):.4f}")

rho     max angle deg  KLT gain  DCT gain  KLT-DCT gap  DFT real  DFT complex  matched
0.5            14.35    1.0932    1.0499      0.04330    0.9387       0.9261  True
0.8            6.448    3.8824    3.8314      0.05097    3.2455       3.1872  True
0.9            3.352    6.3109    6.2761      0.03479    5.3616       5.2684  True
0.95           1.708    8.8462    8.8259      0.02030    7.7040       7.5873  True
0.99          0.3465   14.8850   14.8804      0.00460   13.5704      13.4315  True
0.999        0.03475   23.6179   23.6174      0.00047   22.2623      22.1180  True
C1 angle(0.99) / angle(0.999) = 9.9687; first-order prediction 10.0909

As $\rho$ rises, the largest angle between a KLT vector and the matching DCT vector falls at every step, from 14.35° at $\rho = 0.5$ to 0.035° at $\rho = 0.999$, with every KLT vector matched to its DCT counterpart at every $\rho$, and from $\rho = 0.99$ to $0.999$ it shrinks by a factor of 9.97, within about 1% of the factor 10.09 that the first-order rate $(1 - \rho)/\rho$ predicts. The KLT–DCT gap in coding gain falls at every step from $\rho = 0.8$ to $0.999$, where it is below 0.001 dB. It is smaller at $\rho = 0.5$ than at $0.8$ because it must also vanish as $\rho \to 0$, where $R$ becomes the identity and every orthonormal transform is a KLT. At every correlation the DCT’s gain exceeds both DFT accountings, by more than 1 dB from $\rho = 0.95$ on. The gains themselves grow quickly with $\rho$: the more correlated the pixels, the more a transform can concentrate. The matching check held at every $\rho$, and all four planned hypotheses held: the angle falls at every step, the ratio is within the planned 10% of the first-order prediction, the gap falls from $\rho = 0.8$ on and ends below 0.001 dB, and the DCT’s gain exceeds both DFTs’ at every $\rho$. The correlations were chosen at planning to approach 1, so this part illustrates the derivation above rather than testing it independently.

These results hold for the AR(1) model, which describes a smooth texture with exponentially decaying correlation and no edges. For blocks with other statistics the DCT is not guaranteed to be near the KLT; its appeal is that it is fixed, fast and close to optimal for the smooth, highly correlated content that the model captures.

Common Misconceptions

• “The DCT is the real part of the DFT.” The real part of an $N$-point DFT uses cosines that complete whole cycles over the block and that, without the sines, do not form a basis. The DCT-II corresponds to the DFT of the mirrored sequence of length $2N$, with a half-sample phase shift, and its index $k$ counts half-cycles.
• “The DCT compresses an image.” The orthonormal DCT is invertible and preserves energy; by itself it saves nothing. Compression comes from quantizing the coefficients and entropy coding the result, the stages that follow the transform in JPEG (CCITT, 1992, section 4.3). The transform only makes those stages effective, by concentrating energy.
• “scipy.fft.dct is orthonormal.” Its default, norm="backward", is unnormalized, with a factor 2 in every coefficient; norm="ortho" gives the orthonormal transform that JPEG’s formula defines. Mixing the two changes coefficients by factors of $\sqrt{2N}$ and, for the DC term, $2\sqrt N$.
• “The DCT is the optimal transform for images.” The KLT is optimal for a known covariance, and the DCT approaches it only for highly correlated AR(1)-like statistics, as the examples show. The DCT’s advantages are that it is fixed, real and fast.
• “All DCTs are the same.” The types differ in their boundary conditions and their bases, and SciPy’s type argument selects among four of them. Applications use the DCT-II for blocks and the DCT-IV in lapped transforms (Strang, 1999).

Where It Is Used

• Image and video compression. JPEG’s baseline and extended lossy processes transform $8 \times 8$ blocks with the DCT (CCITT, 1992), as Digital Camera and Image Formation describe in the camera pipeline; quantization, zig-zag scanning and rate belong to image compression. Strang’s 1999 preprint also reports that Malvar’s modulated lapped transform was heavily used by the Sony MiniDisc and Dolby AC-3 and included in the MPEG-4 standard for video.
• Image quality experiments. Image Quality Assessment, Mean Squared Error, SSIM and LPIPS each build a distortion by quantizing orthonormal $8 \times 8$ block DCTs and score it with their metrics: the first three level-shift by 128 and use scaled copies of T.81’s example luminance table (Table K.1, which the standard gives as an example only), and LPIPS uses one step for every coefficient.
• Denoising. Thresholding transform coefficients removes noise when the signal is compact in the transform and the noise is spread over all coefficients; Image Denoising describes this for wavelet shrinkage and BM3D’s 3D transform, and Non-Local Means reports its authors’ comparison with a DCT empirical Wiener filter.
• Relatives. Principal component analysis computes the data-dependent transform that the DCT approximates, and the Fourier transform is the periodic counterpart from which the DCT is computed.