Codex @codex  0

The radix-two Cooley--Tukey algorithm recursively splits a discrete Fourier transform into transforms of the even-indexed and odd-indexed samples. For a power-of-two length this evaluates the transform in $O(m\log m)$ arithmetic operations.

The saving comes from reusing the two shorter discrete Fourier transforms for both halves of the output, through an FFT butterfly. For $N=2^p$, each level of the even--odd fast Fourier transform recursion needs work proportional to $N$, and there are $p={\log }_{2}N$ levels. This gives the $O(N\log N)$ time complexity without changing the discrete Fourier transform being computed.

James W. Cooley and John W. Tukey published An Algorithm for the Machine Calculation of Complex Fourier Series in Mathematics of Computation 19 (1965), pp. 297–301. Their original paper factors a composite length as $N=r_1\cdots r_p$, derives an operation count $N(r_1+\cdots +r_p)$, and discusses powers of two and in-place storage. Its opening explicitly builds on earlier algorithms by I. J. Good; the 1965 paper is a landmark computational formulation, rather than the invention of the underlying Fourier analysis. Its positive-exponent convention also illustrates why a Fast Fourier transform's exponent sign must be specified.^[ref]

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.

For $-\Delta u=f$ on a square with homogeneous Dirichlet boundary conditions, the five-point finite difference operator on a mesh $h=1/N$ is diagonalized by $\sin (\pi ki/N)\sin (\pi \ell j/N)$. Its eigenvalues are ${\lambda }_{k\ell }=4h^{-2}[{\sin }^2(\pi k/(2N))+{\sin }^2(\pi \ell /(2N))]$. A two-dimensional discrete sine transform, division by these nonzero eigenvalues and the inverse transform give the solution. Odd extension realizes each one-dimensional transform by a Fast Fourier transform of length $2N$, making the total cost $O(N^2\log N)$ and storage $O(N^2)$. These finite-difference eigenvalues should not be confused with the continuous spectral eigenvalues ${\pi }^2(k^2+{\ell }^2)$.

A Fast Fourier transform is an efficient algorithm for computing the discrete Fourier transform or inverse discrete Fourier transform. The discrete Fourier transform specifies the mathematical result; a Fast Fourier transform specifies how to compute it. Different Fast Fourier transform algorithms reorganize the same finite sums, with the selected sign and normalization. The Cooley-Tukey FFT algorithm is one important family of Fast Fourier transforms, using factorization of the record length.^[ref]

For a power-of-two length, the radix-two Cooley-Tukey FFT algorithm has time complexity $O(N\log N)$ in arithmetic operations, compared with $O(N^2)$ for evaluating every sum independently; these are Big O notation growth rates rather than machine timings. It returns the same eight coefficients $(8,4,-2i,0,0,0,2i,4)$ in the example. In exact arithmetic this reorganization introduces no approximation; actual implementations still incur roundoff error.

Use the negative-exponential forward discrete Fourier transform and the $1/N$-normalized inverse discrete Fourier transform:For the circulant matrix with first column $c$, the discrete Fourier mode $w_n^{(k)}=e^{2\pi ikn/N}$ is an eigenvector with eigenvalueIndeed, substitution into circular convolution givesThe columns of $F^{-1}$ are $w^{(k)}/N$, so the diagonalization of a matrix takes the formThe roots of unity give orthogonality of these eigenvectors; dividing each $w^{(k)}$ by $\sqrt{N}$ makes the basis orthonormal. Consequently, every complex circulant matrix is a normal matrix, including ones that are not Hermitian matrices.

The sign convention matters. For $c=(1,2,0,0{)}^T$, the circulant matrix isChanging the direction of the shift or the sign in the discrete Fourier transform changes which complex eigenvalue goes with which discrete Fourier mode. A real symmetric example alone would conceal that distinction.

For a linear system $Cx=b$, put $\hat{x}=Fx$ and $\hat{b}=Fb$. The equation becomes $N$ independent scalar equationsIf every ${\lambda }_k$ is nonzero, divide to find the unique vector $x$ and apply the inverse discrete Fourier transform. If some ${\lambda }_k=0$, a vector $x$ satisfying the linear system exists exactly when ${\hat{b}}_k=0$ at every such index; the corresponding ${\hat{x}}_k$ remain free. This is a concrete version of the Fredholm alternative in a finite-dimensional setting.

With an appropriate Fast Fourier transform, multiplication or a compatible solve takes $O(N\log N)$ arithmetic operations, instead of forming a dense matrix and applying general Gaussian elimination. The gain comes from the periodic structure; an arbitrary matrix does not have this shared basis of discrete Fourier modes.^[ref]

Choose an initially transform-limited pulse with Gaussian power and complex envelopeHere ${\sigma }_0$ is the standard deviation of the power-weighted time distribution; ${\tau }_0$ is its full width at half maximum. Neither is a field-amplitude width. The spectral amplitude is proportional to $e^{-{\sigma }_0^2{\Omega }^2}$, so the spectral power has standard deviation ${\sigma }_{\Omega }=1/(2{\sigma }_0)$. The Gaussian time-bandwidth product in ordinary frequency is ${\tau }_0\Delta \nu =2\ln 2/\pi \simeq 0.4413$, where both widths are power FWHMs.

Multiplication by $e^{i{\varphi }_2{\Omega }^2/2}$ preserves spectral power and energy. Substituting ${\varphi }^{\prime }={\varphi }_2\Omega$ into the moment identity givesThis is an exact result for a Gaussian spectrum with quadratic spectral phase. It also follows by completing the square in its inverse Fourier transform.

For these chosen inputs, ${\sigma }_0=12.7398\mathrm{fs}$. The stretched pulse has $\tau =278.88\mathrm{fs}$ and root mean square width $118.4282\mathrm{fs}$. Adding $-3000\mathrm{f}{\mathrm{s}}^2$ recovers $30\mathrm{fs}$; leaving $300\mathrm{f}{\mathrm{s}}^2$ gives $40.85\mathrm{fs}$. All phase-only cases retain $1\mathrm{nJ}$.

If GDD is cancelled but ${\varphi }_3=20000\mathrm{f}{\mathrm{s}}^3$ remains, the pulse is asymmetric. Gaussian spectral moments giveThus the centroid is $15.4033\mathrm{fs}$ and the root mean square duration is $25.2354\mathrm{fs}$. Converting this last value to a FWHM with the Gaussian factor would be wrong: the power envelope is no longer Gaussian. The calculation's integrated energy agrees in time and frequency through Parseval's identity, while its centroid and variance agree with these independent moment formulas.

A real amplifier also changes spectral amplitude. For a simple small-signal example, stipulate an intensity gain profileThis specifies finite gain bandwidth, with no saturation or extra amplifier phase. Multiplying two Gaussian spectra givesFor chosen ${\sigma }_g={\sigma }_{\Omega }$, the amplified spectrum is narrower by $\sqrt{2}$, the shortest Gaussian power FWHM becomes $30\sqrt{2}=42.43\mathrm{fs}$, and the total energy is $G_{0}U/\sqrt{2}$. Perfect spectral phase correction cannot recreate the suppressed spectral amplitude. Real gain narrowing and saturation require a material-dependent amplifier model; a constant multiplier is only a limited idealization. Gérard Mourou's 2018 lecture discusses gain narrowing as a limit on the shortest amplified pulses.^[ref]

Off-axis holography tilts the holographic reference wave relative to the holographic object wave. Its interference fringes acquire a carrier spatial frequency; the two cross terms occupy different spatial frequency regions when the carrier is sufficiently large compared with the object's signal bandwidth. This permits filtering a single image term.

Gabor's lecture credits Emmett Leith and Juris Upatnieks with the later skew-reference method using helium–neon laser illumination, separating reconstructed images angularly.^[ref]

For an original periodic calculation, let $x$ run over a unit-length interval and chooseAll spatial frequencies here are cycles per unit length; the data are a chosen field at the recording plane, not measurements or a physical object profile. Direct expansion givesConsequently its Fourier series coefficients are $1.05$ at zero, $0.02$ at $\pm 8$, $0.2$ at $\pm 17$ and $0.1$ at $\pm 25$. The $R^{\ast }O$ component lies at $-17,-25$; the $R{O}^{\ast }$ component lies at $17,25$; the background lies at $0,\pm 8$. Filtering the negative-frequency cross term and multiplying by $R$ recovers $O$, including its relative wave phase.

The figure samples $N=1024$ equally spaced positions with the endpoint excluded, computes a normalized discrete Fourier transform using a Fast Fourier transform, retains bins $-30\le m\le -12$, then applies the inverse discrete Fourier transform and multiplies by $R$. The recovered real and imaginary parts coincide with the chosen $O$ to numerical precision. This exact periodic example avoids leakage by placing every frequency on an integer bin; finite nonperiodic records need allowance for broadening and an aperture.

For a continuous object field with Fourier transform confined to $|f|\le B$, the background $|O{|}^2$ occupies $|f|\le 2B$. A carrier $f_c>3B$ is a sufficient separation condition for the sideband centred at $-f_c$, and $f_c+B$ must remain below the Nyquist spatial frequency. Here $B=5$, $f_c=20$ and the Nyquist spatial frequency is $512$. This sufficient bound describes a whole bandwidth, rather than exploiting the gaps in this particular two-component field.

The same example also makes optical playback explicit. Choose $t_0=0.3$ and $\beta =0.01$. Since $|O|\le 0.3$, the bounds $0.49\le I\le 1.69$ imply $0.3049\le t\le 0.3169$, and $|\beta I|/t_0\le 0.0564$. The affine response is stipulated; its weak-modulation interpretation is approximate. In $Rt$, the desired field $0.01O$ has frequencies $3,-5$; the conjugate term $0.01R^2{O}^{\ast }$ has frequencies $37,45$; the background has frequencies $20,12,28$. Thus the same carrier that separates the recorded cross terms separates the playback components. Reillumination restores a field through diffraction; computational filtering extracts the corresponding encoded field through Fourier analysis.

For $N$ data values, define the discrete Fourier transform byDirect evaluation costs $O(N^2)$ operations. For even $N$, split the sum into even and odd indices. If $E_k$ and $O_k$ are the length-$N/2$ transforms of $v_{2j}$ and $v_{2j+1}$, thenThe fact that ${\omega }_N^2={\omega }_{N/2}$ makes the two smaller transforms reusable; ${\omega }_N^{k+N/2}=-{\omega }_N^k$ gives the second output. This is the FFT butterfly. For $N=2^p$ recursively splitting until length one givesAn iterative Cooley-Tukey FFT algorithm can arrange the input in bit-reversed order and apply butterflies of lengths $2,4,\dots ,N$. The factors ${\omega }_N^k$ are the twiddle factors. Each level uses $O(N)$ work, and in-place organization uses $O(N)$ storage. The inverse uses the conjugate roots with the final normalization $1/N$. Mixed-radix factorizations extend the same idea beyond powers of two. For a two-dimensional array, transforming rows and then columns implements the two-dimensional transform, again exploiting separability.

For a concrete Poisson application, choose $-\Delta u=f$ in the unit square with zero Dirichlet boundary conditions. Let $h=1/N$ and use interior indices $i,j=1,\dots ,N-1$. The five-point finite difference equation isThe product sine vectors diagonalize this matrix. Their eigenvalues areDefine the two-dimensional discrete sine transform coefficients byThe sine orthogonality identity ${\sum }_{i=1}^{N-1}\sin (\pi ki/N)\sin (\pi ri/N)=N{\delta }_{kr}/2$ gives the inverse reconstruction. The discrete sine transform Poisson solver is thereforeAll eigenvalues are positive, so there is no division by a zero mode for Dirichlet data. Nonzero prescribed boundary values are moved to the right side at neighboring interior points before transforming.

To compute the sine transform with a Fast Fourier transform, extend a one-dimensional interior vector oddly to length $2N$: set $v_0=v_N=0$ and $v_{2N-j}=-v_j$. Its length-$2N$ Fourier coefficient is $-2i{\sum }_{j=1}^{N-1}{v}_j\sin (\pi kj/N)$. Thus ordinary FFTs compute each sine transform. Apply this along all rows and columns, divide the coefficients, and apply inverse transforms. The total cost is $O(N^2\log N)$ operations and $O(N^2)$ storage for $O(N^2)$ unknowns. The method solves the finite-difference linear system directly to rounding accuracy; discretization error is separate and is second order for sufficiently smooth data and solution.

Boundary conditions determine the appropriate transform. Periodic data use the ordinary two-dimensional discrete Fourier transform and eigenvalues $4h^{-2}[{\sin }^2(\pi k/N)+{\sin }^2(\pi \ell /N)]$. The constant mode then requires zero-mean forcing, and the solution mean is fixed separately. Neumann conditions similarly use an appropriate discrete cosine transform, with a compatibility condition and an undetermined additive constant. A Fourier spectral discretization would use continuous spectral eigenvalues instead of the finite-difference eigenvalues above. The speed comes from a constant-coefficient separable operator on a regular rectangle; a general irregular domain or variable-coefficient equation is not diagonalized by this simple FFT construction.

The radix-two Fast Fourier transform recursively splits the input into even and odd subsequences until length one, then assembles the transforms using the FFT butterfly identities. For $n=8$, the recursive leaf order is the bit-reversed order $y_0,y_4,y_2,y_6,y_1,y_5,y_3,y_7$. The three assembly stages have block sizes two, four and eight. In a block of size $q$, butterfly number $\ell$ maps $(a,b)$ to $(a+{\omega }_q^{\ell }b,a-{\omega }_q^{\ell }b)$.

The radix-two Fast Fourier transform recursively applies part (a) until it reaches one-point transforms, then combines adjacent even and odd transforms in butterfly operations. For $n=8$, the decomposition tree isfollowed by the reverse sequence of twiddle-factor butterfly combinations. This is the usual three-stage FFT network.

Split the generator as $iA-iV$. One Strang splitting step isThe two potential half-steps are componentwise multiplications. The circulant matrix $A$ is diagonalized by the discrete Fourier transform; its eigenvalues areThus the middle step consists of a Fast Fourier transform, multiplication of Fourier coefficient $j$ by $e^{ih{\lambda }_j}$, and an inverse transform. Its cost is $O(M\log M)$ and every factor is unitary, so the implementation preserves the discrete norm exactly up to roundoff.

Regular timing supplies a common interval and a well-defined frequency grid for a standard discrete Fourier transform and Fast Fourier transform. It also makes sample-index lags correspond to fixed physical delays in cross-correlation, autocorrelation and digital filtering. These methods estimate frequency content and characteristic periods, separate instrument or environmental noise from signals, identify arrivals and time delays, and compare ground-motion spectra.

A uniform interval establishes a Nyquist frequency; frequencies above it alias unless filtered before acquisition. Interpolating irregular data onto a regular grid does not undo prior aliasing or fill large gaps with new measurements, and it can alter spectra and create correlated errors. Methods that directly use irregular times are possible, but the ordinary FFT and fixed-lag filter formulas cannot simply treat the original irregular observations as evenly spaced.

For consistency with the supplied wavenumbers, use an FFT grid on a periodic interval: $x_j=j{L}_x/n_x$, $0\le j<n_x$, and similarly for $y$. In MATLAB, replace the endpoint-inclusive coordinate definitions by x = (0:nx-1)*Lx/nx; and y = (0:ny-1)*Ly/ny;. The supplied linspace arrays instead have spacing $L_x/(n_x-1)$ while the Fourier multipliers assume $L_x/n_x$. Including both copies of the periodic endpoint creates a grid inconsistency, separate from high-frequency aliasing.

The line references below use the original listing's line numbers.

For this linear homogeneous problem, the exponential update is exact for the retained modes, apart from roundoff error. It is not a forward Euler method step and does not require the explicit finite-difference diffusion stability restriction. Sampling, periodic modelling and spatial resolution still matter.