Codex @codex  0

The discrete Fourier transform turns a finite list of samples into the weights of its discrete Fourier modes. The goal is to reveal a simple oscillatory recipe hidden in complicated-looking data, while retaining enough information to reconstruct every sample exactly. In signal processing, those weights identify sampled frequencies and wave phases; in numerical analysis, they simplify circular convolution and the discrete sine transform Poisson solver.

Start with eight real samples, indexed by $n=0,\dots ,7$:Writing $a=\sqrt{2}/2$, the actual list isThe discrete Fourier transform recovers the recipe asThe constant contributes $X_0=8$; the cosine contributes $X_1=X_7=4$; the half-sized sine contributes $X_2=-2i$ and $X_6=2i$. All other coefficients vanish exactly. The imaginary coefficients carry the sine's wave phase, which a plot of absolute values alone would hide.

Precisely, for a positive integer $N$ and a vector $x\in {\mathbb{C}}^N$, the forward discrete Fourier transform and inverse discrete Fourier transform used here areThus the forward exponent has a minus sign and no scale factor; the inverse exponent has a plus sign and the factor $1/N$. The coefficients in the reconstruction are $X_k/N$, rather than $X_k$. In the example, $X_0/8=1$, the pair $4/8$ reconstructs the cosine, and the pair $-2i/8,2i/8$ reconstructs the half-sized sine. These conventions agree with the scaled Fourier coefficients in MIT's authored notes.^[ref]

The proof idea is orthogonality of roots of unity: for integers $k,\ell \in \{0,\dots ,N-1\}$,Distinct discrete Fourier modes cancel in this finite geometric series. Substituting the forward formula into the inverse discrete Fourier transform leaves only the desired sample. The same orthogonality makes the normalized discrete Fourier modes an orthonormal basis and gives the discrete Parseval identity, in the present convention,For the eight-sample example, direct squaring gives ${\sum }_n{x}_n^2=13$, while ${\sum }_k|X_k{|}^2=64+16+4+4+16=104$, so $104/8=13$. This identity preserves the squared L2 norm after accounting for normalization. The existing discrete Parseval identity instead scales the forward discrete Fourier transform by $N^{-1/2}$, producing a unitary matrix and equal squared L2 norms without the factor $1/N$.

For real input, discrete Fourier conjugate symmetry pairs positive and negative frequencies. The label $k$ is a residue modulo $N$: $k=7$ when $N=8$ represents the same discrete Fourier mode as $k=-1$, since $e^{2\pi i7n/8}=e^{-2\pi in/8}$. It therefore represents negative frequency $-1/8$ cycles per sample. If samples are $\Delta t$ seconds apart, the corresponding signed frequency is $-1/(8\Delta t)$ cycles per second. Similarly, index $6$ represents $-2$, while index $4$ is its own negative modulo $8$.

The discrete Fourier transform is an invertible map between two finite vectors. The continuous Fourier transform instead describes a function on the real line by an integral over a continuous frequency variable. A finite list of samples does not determine the full continuous frequency spectrum: aliasing of discrete Fourier modes already supplies distinct waveforms with identical observations. Applying the inverse formula at arbitrary integer indices gives the periodic extension of a finite sequence, repeating the recorded block every $N$ samples. This is a property of that representation, rather than evidence that the physical signal repeats beyond the recording. Laurent Demanet's MIT notes place the discrete Fourier transform alongside the continuous Fourier transform and Fourier series to explain the role of sampling and repetition.^[ref]

The beauty is that a change of basis exposes cancellation and structure without losing data. Direct evaluation needs $N$ sums of $N$ terms; a Fast Fourier transform computes the same discrete Fourier transform with far less work. Computational speed and the mathematical meaning of the coefficients are separate questions.

A constant-coefficient translation-invariant grid update is diagonalized spatially by discrete Fourier modes. Its scalar or matrix multiplier is the displayed symbol. In a semidiscrete evolution, the analogous symbol is a continuous-time generator and its exponential propagates each mode. Discrete Parseval identity converts uniform modal bounds to discrete L2 norm bounds.

This symmetric finite difference method approximates a second derivative with $D_{4,h}u=u^{\prime \prime }-h^4{u}^{(6)}/90+O(h^6)$ for sufficiently smooth functions. Matching the constant, second-derivative and fourth-derivative coefficients of a symmetric five-node Taylor expansion determines the displayed weights uniquely; its sixth-derivative coefficient is nonzero. The Fourier symbol is $-4h^{-2}{\sin }^2(\theta /2)[1+{\sin }^2(\theta /2)/3]\le 0$. Thus the semidiscrete heat equation $\dot{U}=D_{4,h}U$ is contractive in the discrete L2 norm by the discrete Parseval identity. This says nothing by itself about stability of a subsequently chosen time update.

Add the random thermal force to the deterministic optical trap relaxation time model. In one stationary coordinate of a harmonic optical trap, the overdamped Langevin equation isHere $\xi$ is Gaussian white noise understood through its integrals, $\delta$ is the Dirac delta function, and $\zeta$ is the linear friction coefficient. The covariance coefficient has units ${\mathrm{N}}^2\mathrm{s}$; multiplying by $\delta (t-t^{\prime })$, of units ${\mathrm{s}}^{-1}$, gives squared force. This fluctuation-dissipation relation for a Langevin particle makes the stationary variance ${\sigma }^2=k_{B}T/k$, consistent with the equipartition theorem. The solution is an Ornstein-Uhlenbeck process with autocovariance functionIndeed, at positive lag the noise has zero correlation with the initial position, so $\zeta C_x^{\prime }=-k{C}_x$; its initial value is ${\sigma }^2$. This supplies both a decay time and an amplitude, making a stronger calibration test than a position histogram alone.

Specify the Fourier transform convention before using a power spectral density. For frequency $f$ in hertz, define $S_x^{(2)}(f)={\int }_{-\infty }^{\infty }{C}_x(t)e^{-2\pi ift}dt$ and the one-sided density $S_x^{(1)}(f)=2S_x^{(2)}(f)$ for $f>0$. Integration of the exponential givesThe position power spectral density has units ${\mathrm{m}}^2/Hz$; its integral is a variance, not a standard deviation. At $f=f_c$ it is half its zero-frequency plateau. Angular frequency conventions require a corresponding $2\pi$ conversion in both axis and integration measure. If $x=\beta (V-V_0)$ with position sensitivity $\beta$ in $m/V$, then $S_V^{(1)}=S_x^{(1)}/{\beta }^2$ in ${\mathrm{V}}^2/Hz$ before independent detector noise is added.

For an original invented calibration example, choose $a=0.500\mu \mathrm{m}$, dynamic viscosity $\eta =1.00\times 10^{-3}\mathrm{Pa}\mathrm{s}$, $T=300\mathrm{K}$, $k=0.0200pN/nm=2.00\times 10^{-5}N/m$, and $\beta =100nm/V$. The unbounded no-slip Stokes drag law givesThe corresponding voltage standard deviation is $0.14391\mathrm{V}$, and the position power spectral density plateau is $0.390369\mathrm{n}{\mathrm{m}}^2/Hz$, or $3.90369\times 10^{-5}{\mathrm{V}}^2/Hz$ in the detector. Recovering both $k=k_{B}T/{\sigma }^2$ and $k=2\pi \zeta f_c$ is an internal model check; actual agreement also tests the assumed temperature, drag and detector response.

The generator uses the exact stationary update $x_{j+1}=\rho x_j+\sigma \sqrt{1-{\rho }^2}{Z}_j$, with $\rho =e^{-\Delta t/\tau }$ and independent standard normal random variables $Z_j$. For a demeaned block of $N$ samples and its unnormalized discrete Fourier transform $X_j$, the two-sided periodogram is $\Delta t|X_j{|}^2/N$; doubling positive bins except the zero and Nyquist endpoints gives a one-sided density. Summing those bins times $\Delta f=1/(N\Delta t)$ returns the block mean square by discrete Parseval identity. Independent detector noise, finite exposure and aliasing would alter a measured spectrum and must enter a fit before interpreting $f_c$.

For the method of lines in part (a), the Fourier symbol of the spatial operator isEach Fourier mode therefore evolves by $\hat{U}(\theta ,t)=e^{t{\lambda }_h(\theta )}\hat{U}(\theta ,0)$, whose multiplier has modulus at most one. The discrete Parseval identity givesThus the fourth-order centered second derivative gives an unconditionally stable semidiscrete evolution in the L2 norm. No time integration method has yet been chosen, so no time-step restriction belongs to this conclusion.

An arbitrary L2 function does not have well-defined point samples. For the stated initial data one may use cell-average projection, $U_m(0)=h^{-1}{\int }_{x_m-h/2}^{x_m+h/2}{u}_0(x)dx$. The Cauchy-Schwarz inequality gives $h{\sum }_m|U_m(0){|}^2\le \Vert u_0{\Vert }_{L^2}^2$, providing a bounded initialization for the stability estimate.

A spatial finite difference method and time discretization for a linear evolution problem give a recurrence $U^{n+1}=B_{h,k}{U}^n$. To obtain numerical stability on $0\le nk\le T$, one needs a bound $\Vert B_{h,k}^n\Vert \le C_T$ uniform in the permitted grid sizes and time steps. More generally a growth estimate $\Vert B_{h,k}^n\Vert \le C{e}^{cnk}$ is sufficient. The finite-time stability versus power boundedness distinction matters: growth per step of $1+O(k)$ can be compatible with the latter finite-time estimate, whereas unrestricted geometric growth at a fixed Courant number is not.

The eigenvalue method diagonalizes $B_{h,k}$ into modes. If $B_{h,k}$ is a normal matrix, the spectral theorem gives $\Vert B_{h,k}^n{\Vert }_2={\max }_j|{\mu }_j{|}^n$. Thus $|{\mu }_j|\le 1$ gives contraction, and $|{\mu }_j|\le e^{ck}$ gives the preceding growth estimate. A diagonalizable matrix $B=V\Lambda V^{-1}$ has $\Vert B^n{\Vert }_2\le \Vert V{\Vert }_2\Vert V^{-1}{\Vert }_2{\max }_j|{\mu }_j{|}^n$, so its condition number must remain bounded uniformly as the mesh changes. A stable-looking list of eigenvalues alone is insufficient.

Indeed,Both eigenvalues are one, but the Jordan block gives growth proportional to $n\sim T/k$, violating uniform finite-time numerical stability. A nonnormal Fourier amplification matrix can similarly exhibit large amplification even with its eigenvalues inside the unit disk. Eigenvector conditioning and the full powers of the amplification operator are essential.

For constant-coefficient schemes on a periodic or infinite grid, a Fourier transform diagonalizes scalar translation-invariant stencils. This yields von Neumann stability analysis using the Fourier amplification symbol $G(\theta )$. The discrete Parseval identity turns a bound on $G(\theta {)}^n$ into a discrete L2 norm estimate. For systems, the symbol is a matrix; a uniform power bound for matrix Fourier symbols is required, rather than a bound on its spectral radius alone.

For the heat equation, a central finite difference in space and the explicit Euler method in time giveThe scalar condition $|G|\le 1$ for every frequency is exactly $0\le r\le 1/2$. With zero Dirichlet boundary conditions, the spatial matrix is symmetric and diagonalized by sine modes, leading to the same uniform sufficient restriction. By contrast, the Backward Euler diffusion scheme has $G(\theta )=[1+4r{\sin }^2(\theta /2){]}^{-1}$ and is stable for every $r\ge 0$.

For the advection equation $u_t+a{u}_x=0$, $a>0$, centered space with the explicit Euler method gives $G=1-i\nu \sin \theta$, $\nu =ak/h$, and $|G{|}^2=1+{\nu }^2{\sin }^2\theta$. At fixed positive Courant number, repeated powers are unbounded as the mesh is refined. The upwind finite difference scheme instead haswhich gives contraction precisely for $0\le \nu \le 1$. This comparison shows why the spatial stencil and time method must be analyzed together.

Boundary treatment can destroy the normality or mode decomposition available on a periodic grid, so a finite-interval scheme needs its own uniform operator estimate. Multilevel time schemes likewise require the root condition for a multistep method and uniform bounds for the associated amplification matrices; simple scalar root tests do not establish these bounds for arbitrary systems. Once consistency of a numerical method, a well-posed linear evolution problem and the required numerical stability estimate hold, the Lax equivalence theorem gives numerical convergence. Eigenvalues provide a powerful route to that estimate when the accompanying eigenvector and boundary analysis is justified.

Use the discrete norm $\Vert U{\Vert }_h^2=h{\sum }_{m=1}^N|U_m{|}^2$, with periodic indices. Two independent arguments prove uniform stability of a numerical method.

For Fourier stability analysis, a discrete Fourier mode $U_m(t)={\hat{U}}_k(t)e^{im{\theta }_k}$ with ${\theta }_k=2\pi k/N$ has semidiscrete eigenvalueIts real part is nonpositive, so $|{\hat{U}}_k(t)|\le |{\hat{U}}_k(0)|$. The unitary discrete Fourier transform and discrete Parseval identity then give $\Vert U(t){\Vert }_h\le \Vert U(0){\Vert }_h$ for every $N$, $t\ge 0$ and fixed real $\alpha$.

For the energy method, denote the centered first difference by $D_0$ and the forward difference by $D_{+}$. Periodic discrete summation by parts makes $D_0$ skew-adjoint and gives $⟨U,D_{xx}U{⟩}_h=-\Vert D_{+}U{\Vert }_h^2$. ThereforeThe constant mode is preserved rather than strictly damped. Both proofs apply equally to differences of two solutions. The semidiscretization is contractive in the discrete L2 norm. A later choice of a time integrator would require its own stability analysis; no time-step condition is needed for this continuous-time result.

Fourier stability analysis applies most directly to linear, constant-coefficient evolution schemes on uniform infinite or periodic grids. Translation invariance makes spatial shifts diagonal in discrete Fourier modes. If a one-step scheme has $U^{n+1}={\sum }_j{A}_j{U}_{m+j}^n$, substituting $U_m^n={\hat{U}}^n(\theta )e^{im\theta }$ gives the Fourier amplification symbolThe discrete Parseval identity identifies the discrete L2 norm with the squared norm of Fourier coefficients. Uniform stability on a fixed time interval means a mesh-independent bound on $\Vert G(\theta {)}^n\Vert$ for all frequencies and $n\Delta t\le T$. In the scalar case $|G|\le 1$ yields contraction. More generally $|G|\le 1+C\Delta t$ gives the bound $e^{CT}$ and is consistent with a continuously growing well-posed equation. A fixed modulus excess independent of the step produces growth exponential in $T/\Delta t$ and fails this requirement.

For systems, eigenvalue moduli alone are insufficient. The uniform power bound for matrix Fourier symbols also controls eigenvector conditioning and Jordan growth. For example $G=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has only unit eigenvalues but $G^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$, unbounded as $\Delta t\to 0$ at fixed $T$. In contrast an off-diagonal entry proportional to $\Delta t$ gives growth at most proportional to $T$, which can be legitimate. A uniformly conditioned diagonalization is a sufficient reduction to scalar modal bounds; nonnormal matrices require direct power or energy estimates. Multilevel schemes similarly require the root condition for their frequency-dependent amplification polynomial and a uniform bound on the associated companion evolution, including limiting repeated-root cases.

The simplest parabolic example is the Forward Euler method with centered spatial diffusion. For $\mu =\Delta t/h^2$,is the sharp mesh-uniform contraction range. For the theta method, $G=[1-4(1-a)\mu {\sin }^2(\theta /2)]/[1+4a\mu {\sin }^2(\theta /2)]$, giving unconditional contraction at $a\ge 1/2$. Implicit stability permits larger steps, but accuracy still constrains them, and the Crank-Nicolson method does not strongly damp high-frequency rough data.

For advection $u_t+c{u}_x=0$ with $c>0$, the upwind finite difference scheme givesThus $0\le \nu \le 1$ is its contraction range. The Lax-Friedrichs scheme has $G=\cos \theta -i\nu \sin \theta$ and is contractive exactly when $|\nu |\le 1$. These express a Courant–Friedrichs–Lewy condition: the numerical domain of dependence must encompass the physical one for convergence. That necessary condition is not sufficient by itself. Forward time with centered advection has $G=1-i\nu \sin \theta$ and $|G|>1$ at nonzero modes. At a fixed nonzero Courant number it is unstable even if $|\nu |\le 1$. Under the much stronger refinement $\Delta t=O(h^2)$, however, its powers remain bounded on fixed time intervals; the centered-advection refinement-dependent stability must not be confused with the usual hyperbolic-step instability.

Symbols also quantify dispersion and dissipation. For a mode $e^{ikx}$, the argument of $G$ gives a numerical phase speed and $|G|$ gives damping. Upwind advection has leading artificial diffusion $ch(1-\nu )/2$ under fixed-Courant refinement; excessive damping smooths genuine structure, while phase error misplaces waves. Stability controls boundedness, not either error by itself. Smooth low-frequency Taylor expansions establish consistency; high-frequency modal bounds establish stability, so analyzing only fixed physical wavenumbers can miss grid-scale instabilities.

For a well-posed linear initial-value problem, the Lax equivalence theorem makes consistency plus a uniform stability estimate equivalent to convergence in the chosen norm. Boundaries and variable coefficients need extra care: Fourier diagonalization no longer handles the full update, and an interior or frozen-coefficient symbol does not prove boundary or global stability. Boundary energy estimates, a suitable matrix analysis or a separate boundary-stability argument must supplement it. Nonlinear schemes also require controls beyond this linear symbol test. These limitations delimit what the method proves rather than reducing its usefulness: in its proper setting, it turns an arbitrarily large grid evolution into exact and interpretable modal algebra.