Defines the Fourier transform of an integrable function on the line in the handbook's convention, records the six conventions used by the cited books and the dictionary between them, proves the operational rules for translation, modulation, dilation, reflection, conjugation, differentiation and multiplication by the variable, and tabulates the standard transform pairs.
Let f:R→C be integrable, written f∈L1(R), meaning that f is measurable and ∥f∥1=∫−∞∞∣f(x)∣dx is finite. The Fourier transform of f is the function
f^(k)=F[f](k)=∫−∞∞f(x)e−ikxdx,k∈R,
and the associated inverse transform of a function g∈L1(R) is
F−1[g](x)=2π1∫−∞∞g(k)eikxdk.
That the second really inverts the first is the inversion theorem; nothing below depends on it. The integral defining f^ converges absolutely for every k, since ∣f(x)e−ikx∣=∣f(x)∣, and three properties follow at once:
Boundedness.∣f^(k)∣≤∥f∥1 for every k, so f^ is a bounded function with supk∣f^(k)∣≤∥f∥1.
Continuity.f^ is continuous on R, by dominated convergence (Dominated convergence) with the fixed integrable majorant ∣f∣.
Decay.f^(k)→0 as ∣k∣→∞. This is the Riemann-Lebesgue lemma; Folland (Folland (1999), Thm 8.22(f), p. 250) states it as F(L1)⊂C0, the continuous functions vanishing at infinity. Its Fourier-series counterpart is proved in Mean-square convergence and Parseval's identity.
Nothing more is true in general. Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §1.2, p. 134) emphasise that integrability of f guarantees no rate of decay for f^, so the inverse integral above need not converge; this is the reason the inversion theorem carries an extra hypothesis and the reason the Schwartz class is introduced there. Two convenient sufficient conditions for membership of L1(R) are used throughout: f continuous with ∣f(x)∣≤A/(1+x2), which Stein and Shakarchi call moderate decrease, and f piecewise continuous with compact support.
The transform is linear: F[αf+βg]=αf^+βg^ for α,β∈C and f,g∈L1. It is one instance of the general integral transform g(x)=∫abf(t)K(x,t)dt that Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.1, Eq. (20.1), p. 963) take as the organising idea of their chapter, with kernel K(k,x)=e−ikx and (a,b)=(−∞,∞); the Laplace transform and the Hankel and Mellin transforms are others.
Where the definition comes from. Take the Fourier series of Fourier series on (−l,l), substitute the coefficient integral into the series, and write kn=nπ/l, so that successive kn are spaced Δk=π/l apart:
f(x)=2π1n=−∞∑∞(∫−llf(y)e−iknydy)eiknxΔk.
The sum is a Riemann sum in k. Letting l→∞ turns Δk into dk, the inner integral into f^(k), and the sum into 2π1∫f^(k)eikxdk. Strauss (Strauss (2008), §12.3, pp. 343–344) warns that this limit is delicate, because the interval grows while the terms change; Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1, p. 435) and Hassani (Hassani (2013), §9.2, pp. 277–278) run the same argument and both state that the existence of the limit is assumed, not proved, at that point. The argument is a derivation of the form of the pair, not of its validity.
Only the product of the two prefactors is fixed, at 1/2π; the split between the forward and the inverse transform, and the sign in the exponent, are free. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1, p. 435) say this explicitly, and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.2, p. 965) add the practical warning that formulas taken from different books differ by powers of 2π and must be checked before being combined. Writing f^(k)=∫f(x)e−ikxdx for our transform, the library divides as follows.
Convention
Forward
Inverse
Books
In terms of f^
Ours
∫f(x)e−ikxdx
2π1∫F(k)eikxdk
Strauss §12.3
f^(k)
2π in the exponent
∫f(x)e−2πixξdx
∫f^(ξ)e2πixξdξ
Stein and Shakarchi Ch. 5; Folland §8.3
f^(2πξ)
Symmetric, physics sign
2π1∫f(x)e−ikxdx
2π1∫f~(k)eikxdk
Riley, Hobson and Bence §13.1; Hassani §9.2; Reed and Simon §IX.1
2π1f^(k)
Symmetric, opposite sign
2π1∫f(t)e+iωtdt
2π1∫g(ω)e−iωtdω
Arfken, Weber and Harris §20.2
2π1f^(−ω)
Asymmetric, opposite sign
∫f(x)e+iωxdx
2π1∫f∧(ω)e−iωxdω
Stakgold and Holst §2.4
f^(−ω)
Symmetric, n dimensions
(2π)−n/2∫u(x)e−ix⋅ydnx
(2π)−n/2∫u^(y)eix⋅ydny
Stakgold and Holst §10.3
(2π)−n/2f^(y) at n=1
Two remarks on the table. First, Stakgold and Holst use the fifth line in Chapter 2 and the sixth in Chapter 10, so a formula of theirs must be matched to its chapter; we have followed each where we cite it, and the two are not interchangeable. Second, the entries in the last column are the whole dictionary: a sign change in the exponent replaces k by −k and nothing else, and a prefactor c on the forward transform multiplies the result by c. Three consequences are worth memorising, since they are where conventions bite.
In Arfken's opposite-sign convention the derivative rule reads F[f′](ω)=−iωg(ω), their Eq. (20.55): the sign of the exponent and the sign in the derivative rule always travel together. The symmetric convention is the only one in which the forward and inverse transforms differ solely by the sign of the exponent, which is why books that stress the duality of x and k prefer it; the 2π-in-the-exponent convention is the only one with no stray constants anywhere, which is why books that stress convolution and L2 prefer it. We follow Strauss because it keeps the exponent clean and puts every factor of 2π in one place, on the inverse transform, where a reader can see it.
Let f∈L1(R), let a∈R and let c∈R∖{0}. Then each function on the left is in L1(R) and
#
Function of x
Its transform at k
Extra hypothesis
1
f(x−a) (translation)
e−ikaf^(k)
none
2
eiaxf(x) (modulation)
f^(k−a)
none
3
f(cx) (dilation)
∣c∣1f^(ck)
none
4
f(−x) (reflection)
f^(−k)
none
5
f(x) (conjugation)
f^(−k)
none
6
f′(x) (differentiation)
ikf^(k)
f absolutely continuous, f′∈L1
7
xf(x) (multiplication)
idkdf^(k)
xf(x)∈L1
Rules 6 and 7 iterate: if f is m times differentiable with f,f′,…,f(m) all in L1 then f(m)(k)=(ik)mf^(k), and if xmf(x)∈L1 then f^ is m times continuously differentiable with xmf(k)=imdmf^/dkm. Two corollaries are used constantly:
Smoothness of f forces decay of f^. If f,…,f(m)∈L1 then ∣k∣m∣f^(k)∣=∣f(m)(k)∣≤∥f(m)∥1, so f^(k)=O(∣k∣−m). Decay of f likewise forces smoothness of f^. This is the continuous form of the smoothness-decay correspondence.
Symmetry. If f is real then f^(−k)=f^(k), by rule 5; if f is real and even, f^ is real and even; if f is real and odd, f^ is purely imaginary and odd. In the even and odd cases the transform collapses to the cosine and sine transformsf^(k)=2∫0∞fcoskxdx and f^(k)=−2i∫0∞fsinkxdx, which Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.2, (20.27)–(20.30), p. 971) normalise as separate transforms with factor 2/π and treat as the natural tools on the half-line.
All in our convention f^(k)=∫f(x)e−ikxdx, with a>0. Each is an ordinary convergent integral; the entries for δ, for the constant function and for the Heaviside step belong to Fourier transforms of distributions and are not repeated here.
f(x)
f^(k)
Check
χ[−a,a](x) (square pulse)
k2sinak
value 2a at k=0
(1−∣x∣/a)+ (triangle)
ak24sin2(ak/2)
value a at k=0
e−a∣x∣
a2+k22a
real and even
a2+x21
aπe−a∣k∣
dual of the previous line
e−ax2
aπe−k2/4a
a=21 gives 2πe−k2/2
H(x)e−ax
a+ik1
modulus (a2+k2)−1/2
The Gaussian line is the only self-reciprocal family: at a=21 the pair is e−x2/2↔2πe−k2/2, the same shape in both variables. The square pulse illustrates the opposite extreme. Both show the reciprocity Hassani (Hassani (2013), §9.2, pp. 279–280) draws out: the pulse of half-width a has a transform whose first zeros are at k=±π/a, so halving the width of f doubles the width of f^, and in the Gaussian line a↦λa scales the two widths by λ−1/2 and λ1/2. The uncertainty principle makes this exact.
Idea of the proofOriginal
Every rule in the table is the same two-line calculation: substitute the new function into the defining integral, then move the change from f onto the exponential and read off what the exponential has become. A translation of f turns into a phase, because e−ik(y+a)=e−ikae−iky splits; a dilation turns into a dilation the other way, because the substitution u=cx rescales both the exponent and dx; and differentiation turns into multiplication by ik, because an integration by parts transfers the derivative onto e−ikx, which reproduces itself times −ik, with a second sign change from moving the derivative across. The whole content of the transform is that e−ikx is an eigenfunction of d/dx.
Proof(Operational rules)AdaptedAdapted from Riley, Hobson & Bence (2006), §13.1.5, the integration-by-parts proof of (13.28) and the modulation calculation (13.33)–(13.34), pp. 444–445; Folland (1999), Thm 8.22(a),(b),(d),(e), pp. 249–250; Reed & Simon (1980), §IX.1, the Lemma on p. 319
Throughout, every integral converges absolutely because ∣eiθ∣=1 and the stated hypothesis puts the relevant function in L1.
Rules 1, 4 and 5: a substitution or a conjugation. For rule 1 substitute y=x−a, so dx=dy and e−ikx=e−ikae−iky:
∫f(x−a)e−ikxdx=e−ika∫f(y)e−ikydy=e−ikaf^(k).
Rule 4 is the substitution y=−x, which reverses the limits and supplies a second sign from dx=−dy, leaving ∫f(y)eikydy=f^(−k). For rule 5, conjugation passes under the integral sign:
∫f(x)e−ikxdx=∫f(x)eikxdx=f^(−k).
Rule 2: no substitution at all. The two exponentials combine, eiaxe−ikx=e−i(k−a)x, so the integral is by definition f^(k−a). This is the one rule that needs nothing; it is the dual of rule 1, and it is the calculation Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.5, (13.33)–(13.34), p. 445) perform to show that modulating a signal onto a carrier frequency ωc shifts its whole spectrum by ωc.
Rule 3: the substitution u=cx. Take c>0 first. Then x=u/c, dx=du/c, and the exponent becomes −i(k/c)u:
∫f(cx)e−ikxdx=c1∫f(u)e−i(k/c)udu=c1f^(ck).
For c<0 the same substitution reverses the limits of integration, which contributes one further minus sign, and the two combine into 1/∣c∣. The factor 1/∣c∣ is forced by the case k=0: both sides must then equal ∫f(cx)dx, which is ∣c∣−1∫f(u)du=∣c∣−1f^(0).
Rule 6: integration by parts. Let f be absolutely continuous with f′∈L1. Integrating by parts on [−N,N],
the +ik arising because dxde−ikx=−ike−ikx and the parts formula supplies a further minus sign. For the boundary term: since f′∈L1, f(N)=f(0)+∫0Nf′ has a limit L+ as N→∞, and L+ must be 0, because a non-zero limit would make ∫∞∣f∣ diverge, contradicting f∈L1. The same argument at −∞ gives L−=0. Hence the boundary term vanishes in the limit and f′(k)=ikf^(k). This is the argument Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.5, p. 444) give, with the vanishing of f at ±∞ justified from the finiteness of ∫∣f∣; Folland (Folland (1999), Thm 8.22(e), p. 250) states it with the explicit hypothesis that f and its derivatives up to order k−1 lie in C0, which is what the limit argument establishes here.
Rule 7: differentiation under the integral sign. Suppose xf(x)∈L1. For h=0,
hf^(k+h)−f^(k)=∫f(x)e−ikxhe−ihx−1dx.
The bracket tends pointwise to −ix as h→0, and ∣e−ihx−1∣≤∣hx∣ by the mean value theorem applied to t↦e−itx, so the integrand is dominated by ∣xf(x)∣∈L1 uniformly in h. Dominated convergence (Dominated convergence) gives
dkdf^(k)=∫(−ix)f(x)e−ikxdx=−ixf(k),
and multiplying by i gives rule 7. Continuity of df^/dk follows by applying the continuity statement of the definition to xf∈L1. Iterating rules 6 and 7 gives the m-fold versions; Reed and Simon (Reed & Simon (1980), §IX.1, Lemma, p. 319) package both directions as one identity on multi-indices, which is the form in which the two rules are plainly dual. □
Remark(Where the hypotheses are used)Original
Integrability is used for everything, and it is all that rules 1–5 use. Those five are substitutions in an absolutely convergent integral; they would fail for a merely conditionally convergent integral, where a change of variable can change the value.
Rule 6 needs f′∈L1, not merely that f′ exists. The hypothesis enters twice: to make f′ defined at all, and to force f(±∞)=0 so that the boundary term dies. Without the second, the rule acquires a boundary term. Take f=H(x)e−x, which is in L1 with f^(k)=1/(1+ik) but is not absolutely continuous at 0: its pointwise derivative off the origin is −f, whose transform is −1/(1+ik), whereas ikf^(k)=ik/(1+ik). The two differ by exactly 1, the size of the jump, which is the delta that the pointwise derivative omits (The Heaviside step and related functions).
Rule 7 needs xf(x)∈L1, which is strictly stronger than f∈L1. For f(x)=(1+x2)−1 we have f∈L1 but xf∈/L1, and indeed f^(k)=πe−∣k∣ is not differentiable at k=0. The failure is exactly at the one point the hypothesis was needed for.
Nothing here needs f^ to be integrable, and in general it is not. The square pulse has f^(k)=2sinak/k, which decays like 1/∣k∣ and so is not in L1. Riemann-Lebesgue gives decay but never a rate; a rate costs smoothness of f, by the first corollary to the rules.
The Riemann-Lebesgue lemma is not reversible.F maps L1 into C0 but not onto it, as Reed and Simon (Reed & Simon (1980), §IX.2, Thm IX.7 and the remark following it, p. 327) note; there are continuous functions vanishing at infinity that are the transform of no integrable function.
At a jump the transform cannot see which value f takes. Changing f at finitely many points changes no integral, so f^ is unchanged; correspondingly the inversion integral returns the midpoint 21[f(x+)+f(x−)] rather than any particular value, as Hassani (Hassani (2013), §9.2, Eq. (9.22), p. 280) records.
Method(Computing a Fourier transform)Original
Fix the convention before writing anything. Decide where the 2π sits and which sign the exponent carries, and if a formula is being imported, convert it using the last column of the conventions table. Almost every discrepancy between a textbook answer and a computed one is a factor of 2π or a sign of k.
Check integrability. If ∫∣f∣=∞ the transform is not defined by the integral; move to Fourier transforms of distributions, where 1, H and δ acquire transforms.
Exploit symmetry first. If f is even, write f^(k)=2∫0∞fcoskxdx; if odd, f^(k)=−2i∫0∞fsinkxdx. This halves the work and tells you in advance whether the answer is real or imaginary.
Try to reach a tabulated pair by the operational rules. Translation, modulation, dilation and the two differentiation rules generate a large family from each entry of the table; this is almost always faster than integrating.
Otherwise integrate directly. Split the range at any point where the formula for f changes (typically at 0 for anything involving ∣x∣ or H). For a Gaussian exponent, complete the square and quote Gaussian integrals; the resulting shift of the contour into the complex plane is justified by Cauchy's theorem, as Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Example 20.2.2 and Fig. 20.3, pp. 968–969) set out. For a rational f, close the contour in the half-plane where e−ikx decays — the lower one for k>0, the upper one for k<0 — and use The residue theorem.
Check three things. The value at k=0 must be ∫f. If f is real, f^(−k) must equal f^(k). And f^ must tend to 0 at ±∞ at the rate the smoothness of f predicts: a jump in f gives 1/∣k∣, a kink gives 1/k2.
Example(A modulated, shifted exponential)Original
Compute the transform of
f(x)=e−2∣x−3∣cos(5(x−3)),
a damped oscillation of decay rate 2 and angular frequency 5, centred at x=3. It is continuous and ∣f(x)∣≤e−2∣x−3∣, so f∈L1(R) with ∥f∥1≤1, and the transform is defined.
Step 1: the base pair. From the table, g(x)=e−2∣x∣ has g^(k)=4/(4+k2). Derived directly: g is even, so g^(k)=2∫0∞e−2xcoskxdx, and integrating by parts twice (or taking the real part of ∫0∞e−(2+ik)xdx=1/(2+ik)) gives 2Re2+ik1=4+k24.
Step 2: the modulation. Write cos5y=21(e5iy+e−5iy), so that g(y)cos5y=21e5iyg(y)+21e−5iyg(y). Rule 2 applied to each term, with a=5 and a=−5, gives for h(y)=e−2∣y∣cos5y
Two Lorentzian peaks of half-width 2 sit at k=±5: the oscillation has put the energy at the carrier frequency, and the damping has given each line a width.
Step 3: the translation.f(x)=h(x−3), so rule 1 with a=3 multiplies by the phase e−3ik:
f^(k)=e−3ik[4+(k−5)22+4+(k+5)22].
The shift changes only the phase, never ∣f^∣ — the physical reading is that moving a pulse in time does not change its power spectrum.
Step 4: a number. At k=1 the bracket is 4+162+4+362=101+201=203=0.15, so
Numerical integration of ∫−∞∞e−2∣x−3∣cos(5(x−3))e−ixdx, split at x=3, returns −0.1484988745−0.0211680012i, agreeing to ten decimal places.
Step 5: the checks of the method.f^(0)=292+292=294=0.137931, which must equal ∫f; direct integration of ∫e−2∣y∣cos5ydy=2⋅4+252 confirms it. And f is real, so f^(−k)=f^(k) must hold: at k=−1 the bracket is unchanged (it is even in k) and the phase conjugates, as required.
ApplicationOriginal
Spectroscopy and line shapes. The example above is the standard model of a radiating atom: an oscillation at the transition frequency ω0, damped at rate γ by spontaneous emission. Its transform is a pair of Lorentzians of half-width γ centred at ±ω0, and the measured line width is therefore the inverse lifetime. Nothing but the operational rules is needed to get there, and the same calculation with a Gaussian damping — Doppler broadening instead of radiative damping — gives a Gaussian line, which is how the two broadening mechanisms are told apart in a spectrum.
Fraunhofer diffraction. For a slit of half-width Y illuminated normally by light of wave number K, the amplitude observed at angle θ is, in the far field, the Fourier transform of the aperture function f(y) evaluated at the spatial frequency k=Ksinθ; Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.2, pp. 437–439) derive this from the superposition of Huygens wavelets. A uniform slit is the square pulse of the table, so the amplitude is 2sin(kY)/k and the intensity the familiar sinc2 pattern, with first zeros at sinθ=±π/(KY). Rule 1 says that displacing the slit changes only a phase, so the intensity pattern is independent of where the aperture sits; rule 3 says that a narrower slit diffracts more widely.
Solving differential equations. Rule 6 turns d/dx into multiplication by ik, so a linear constant-coefficient differential equation on the line becomes an algebraic equation for u^, which is why the transform is the standard tool for the heat, wave and Poisson equations on unbounded domains; Hassani (Hassani (2013), §9.2.1, pp. 284–285) works this through for a general second-order ordinary differential equation and obtains y~(k)=f~(k)/(−C2k2+iC1k+C0). The cost, as Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.3, pp. 984–985) note, is that the method presupposes the solution is unique and decays at infinity, which is really a boundary condition in disguise.
Momentum space in quantum mechanics. For a wave function ψ(x), ψ^(k) is the amplitude for the particle to be found with momentum ℏk. Rule 6 is then the statement that the momentum operator −iℏd/dx acts as multiplication by ℏk in the transform variable, and rule 7 that position acts as id/dk; the two rules together are the canonical commutation relation. Rule 3, read with the Gaussian line of the table, is already the statement that a narrow wave packet has a broad momentum distribution, made quantitative in The uncertainty principle.
Proves that a rapidly decreasing function is recovered from its Fourier transform by the inverse integral, using the Gaussian as an approximate identity and the multiplication formula, shows that the transform is a bijection of the Schwartz class whose square is reflection times a constant, and records the weaker versions that hold when the transform is merely integrable or the function merely piecewise smooth.
The Schwartz spaceS(R) consists of the functions f:R→C that are infinitely differentiable and satisfy
x∈Rsupxmf(n)(x)<∞for all integers m,n≥0.
Equivalently: f and every derivative decay faster than any power of 1/∣x∣. Such f are called rapidly decreasing. The Gaussians e−ax2 with a>0 are in S, as is any smooth compactly supported function; e−∣x∣ is not (it is not differentiable at 0) and neither is (1+x2)−1 (it decays too slowly).
Three facts are used below and all follow directly from the definition.
S⊂L1. Taking n=0 with m=0 and with m=2 gives two bounds, ∣f(x)∣≤C0 everywhere and ∣f(x)∣≤C2/x2 for x=0. Put C=2max(C0,C2): for ∣x∣≤1 the first bound gives ∣f(x)∣≤C0≤C/2≤C/(1+x2), and for ∣x∣≥1 the second gives ∣f(x)∣≤C2/x2≤C/(2x2)≤C/(1+x2). So ∣f(x)∣≤C/(1+x2), which is integrable. Every f∈S therefore has a Fourier transform in the sense of The Fourier transform.
S is stable under differentiation and under multiplication by a polynomial, directly from the definition, hence under any operator xm(d/dx)n.
Every f∈S is bounded and uniformly continuous, since f and f′ are bounded.
The two Gaussian kernels used in the proof are, for ε>0,
Gε(k)=e−εk2/2,Kε(u)=2πε1e−u2/(2ε).
Gε is a convergence factor that tends to 1 pointwise as ε→0+; Kε is a probability density of variance ε that concentrates at the origin. Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, Thm 1.6, p. 139) call a family with the properties of Kε a family of good kernels.
Theorem(Fourier inversion on the Schwartz class)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, Thm 1.3, p. 137, Thm 1.9 and Cor. 1.10, pp. 141–142; Reed & Simon (1980), §IX.1, Thm IX.1, p. 320; Folland (1999), §8.3, Cor. 8.28, p. 252
Let f∈S(R) and let f^(k)=∫f(x)e−ikxdx. Then
The class is preserved.f^∈S(R). In particular f^∈L1, so the integral in part 2 converges absolutely.
Inversion. For every x∈R,
f(x)=2π1∫−∞∞f^(k)eikxdk.
Bijectivity.F:S→S is a linear bijection, with inverse F−1[g](x)=2π1∫g(k)eikxdk. Writing (Rf)(x)=f(−x) for reflection,
F2=2πR,F4=(2π)2I,
so up to the constant the transform has order four, and its only possible eigenvalues on S are the fourth roots of (2π)2, namely ±2π and ±i2π.
Theorem(Two weaker versions)AdaptedAdapted from Folland (1999), §8.3, Thm 8.26 and Cor. 8.27, pp. 251–252; Stakgold & Holst (2011), §2.3, p. 130 and §2.4, (2.4.5) and the Fourier Integral Theorem (2.4.6), pp. 146–147; Stein & Shakarchi (2003), Ch. 5, §1.7, p. 144; Hassani (2013), §9.2, (9.22), p. 280
The integrable version. If f∈L1(R) and f^∈L1(R), then some continuous f0 satisfies f=f0 off a set of measure zero, and for every x
f0(x)=2π1∫−∞∞f^(k)eikxdk.
Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §1.7, p. 144) record the same extension for functions of moderate decrease whose transforms also have moderate decrease, and remark that it needs no new ideas.
Uniqueness. In particular, if f∈L1 and f^≡0 then f=0 almost everywhere: two integrable functions with the same transform are equal almost everywhere.
The piecewise-smooth version.Stakgold and Holst (Stakgold & Holst (2011), §2.4, Fourier Integral Theorem (2.4.6), p. 146) ask for no smoothness beyond piecewise regularity, and for no decay of f^ at all. A function obeys their Dirichlet conditions (Stakgold and Holst (Stakgold & Holst (2011), §2.3, p. 130)) when it is piecewise continuous and the interval admits a finite subdivision inside each piece of which f has a continuous derivative, with one-sided derivatives at the subdivision points; this is the same hypothesis under which the Fourier series of a periodic function converges, treated in Pointwise convergence of Fourier series. If f∈L1(−∞,∞) and f satisfies Dirichlet conditions, then for every x
2f(x+)+f(x−)=R→∞lim2π1∫−RRf^(k)eikxdk,
the limit being over symmetric intervals, which they stress is how the integral must be read. Their statement is written with the kernel e+iωξ and the inverse kernel e−iωx; substituting ω=−k turns it into the display above, with no other change. Hassani (Hassani (2013), §9.2, Eq. (9.22), p. 280) states the same midpoint rule for piecewise continuous f and f^ in the symmetric convention.
Idea of the proofOriginal
We would like to substitute the definition of f^ into the inversion integral and exchange the two integrations, whereupon 2π1∫eik(x−y)dk would act as a delta function and pick out f(x). That is exactly the calculation Hassani (Hassani (2013), §9.2, pp. 280–281) performs, and it is the right picture; but the exchange is illegitimate, because the double integrand has modulus ∣f(y)∣, which is not integrable over the whole plane. Folland (Folland (1999), §8.3, p. 251) names this as precisely the obstruction.
The repair is to insert a convergence factor e−εk2/2, which makes the double integral absolutely convergent so that Fubini applies, and then to let ε→0+. The Gaussian is the right factor for two reasons: it is its own transform up to scaling, so inserting it in the k integral produces a Gaussian in the y integral, and that Gaussian is positive with total mass one, so as ε→0+ it concentrates at y=x and reproduces f(x). In short, the unproved delta function of the formal calculation is replaced by an honest, explicitly computable approximation to it.
Let f∈S and fix integers m,n≥0. By item 2 of the definition, xnf(x) lies in S, hence in L1, so rule 7 of The Fourier transform may be applied n times:
dkndnf^(k)=(−ix)nf(k).
Put g(x)=(−ix)nf(x), again in S, so g and all its derivatives are in L1. Rule 6 applied m times gives g(m)(k)=(ik)mg^(k), that is
kmdkndnf^(k)=(−i)mdxmdm[(−ix)nf(x)](k).
The function being transformed on the right is in S, hence in L1, so by the boundedness statement of The Fourier transform the right-hand side is bounded by dxmdm[xnf]1<∞, uniformly in k. As m and n were arbitrary, f^∈S. The only two inputs are the two operational rules, and they enter in opposite directions: rule 7 converts derivatives in k into powers of x, rule 6 converts powers of k into derivatives in x. □
Proof(Part 2: inversion)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, Thm 1.4, pp. 138–139, Thm 1.6 and Cor. 1.7, pp. 139–140, Prop. 1.8, pp. 140–141, and Thm 1.9, p. 141; Folland (1999), §8.3, Lemma 8.25 and Thm 8.26, pp. 251–252
Step 1: the Gaussian transform. For a>0, completing the square as in Gaussian integrals and shifting the contour gives
∫−∞∞e−ax2e−ikxdx=aπe−k2/(4a).
Specialising to a=ε/2 and writing the transform variable as u,
Gε(u)=∫e−εk2/2e−ikudk=ε2πe−u2/(2ε)=2πKε(u).
This is the one place a specific function is computed, and it is why the Gaussian is used: it is the only elementary family that reproduces itself. (Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, Thm 1.4, p. 138) obtain it instead from the differential equation F′(ξ)=−2πξF(ξ) satisfied by the transform of e−πx2, which avoids the contour shift.)
Step 2: Kε is an approximate identity. Three properties, all immediate from step 1 and the substitution u=εv:
Kε≥0;
∫Kε(u)du=2π1∫e−v2/2dv=1, by the substitution u=εv and the normalisation of Gaussian integrals;
for each fixed η>0, ∫∣u∣>ηKε(u)du=2π1∫∣v∣>η/εe−v2/2dv⟶0 as ε→0+, since the region of integration escapes to infinity.
Claim. If h is bounded and continuous at x, then ∫h(y)Kε(y−x)dy→h(x) as ε→0+.
Proof of the claim. Let B=sup∣h∣ and let σ>0. Using the second property to write h(x)=∫h(x)Kε(y−x)dy,
where positivity of Kε was used to take the modulus inside. Continuity of h at x lets us choose η so that the first integrand is at most σ on its range, making the first term at most σ; the third property then makes the second term at most σ for all small enough ε. As σ was arbitrary, the claim follows. This is the argument Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, Cor. 1.7, p. 140) give for uniform convergence when h is uniformly continuous; we need it only at a single point.
Step 3: the multiplication formula. For p,q∈L1(R),
∫p(y)q^(y)dy=∫p^(k)q(k)dk.
Both sides are the double integral ∬p(y)q(k)e−ikydkdy, taken in the two possible orders. The double integrand has modulus ∣p(y)∣∣q(k)∣, whose integral over R2 is ∥p∥1∥q∥1<∞, so Tonelli's theorem licenses Fubini's theorem (Fubini and Tonelli theorems) and the two orders agree. This is Folland's (Folland (1999), Lemma 8.25, p. 251) lemma; note that it needs nothing beyond integrability, because the Gaussian is not yet present.
Step 4: putting them together. Fix x∈R and ε>0, and apply step 3 with p=f and q(k)=eikxGε(k), both in L1. By rule 2 of The Fourier transform (modulation, with the roles of the variables exchanged) and step 1,
q^(y)=∫eikxe−εk2/2e−ikydk=Gε(y−x)=2πKε(y−x).
The multiplication formula therefore reads
2π∫f(y)Kε(y−x)dy=∫f^(k)eikxe−εk2/2dk.
Now let ε→0+ on both sides.
On the left, f is bounded and continuous (item 3 of the definition), so the claim of step 2 gives 2πf(x).
On the right, the integrand tends pointwise to f^(k)eikx and is dominated in modulus by ∣f^(k)∣, which is integrable by part 1. Dominated convergence (Dominated convergence) gives ∫f^(k)eikxdk.
Equating the two limits and dividing by 2π proves part 2. □
Proof(Part 3: bijectivity and the order-four identity)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, §1.5, pp. 141–142, Cor. 1.10; Folland (1999), §8.3, Cor. 8.28, p. 252
By part 1, F maps S into S, and the same proof applied to the inverse integral shows that F−1 does too, since the two differ only by the sign of the exponent and the constant. Part 2 says F−1F=I on S, which makes F injective. Comparing the two integrals gives F−1[g](x)=2π1g^(−x), that is F−1=2π1RF; since R is a bijection of S and F is injective, F−1 is injective too. Now take g∈S and apply part 2 to F−1[g]: it gives F−1[FF−1[g]]=F−1[g], and injectivity of F−1 cancels the outer F−1, leaving FF−1[g]=g. So F is surjective. Hence F is a bijection of S with the stated inverse.
For the order-four identity, apply F to f^ and compare with the inversion formula read at −x:
So F2=2πR, and since R2=I we get F4=(2π)2I. If Ff=λf with f=0 then λ4=(2π)2, giving the four listed values. (In the symmetric convention of Reed and Simon (Reed & Simon (1980), §IX.1, p. 318) the constant disappears and F4=I exactly; this is one reason that convention is preferred when the transform is studied as an operator.) □
Remark(Where the hypotheses are used)Original
Integrability of f^ is used once, in the dominated convergence step, and it is exactly what can fail. Take f=χ[−1,1], which is in L1 but not continuous and not in S. Its transform f^(k)=2sink/k decays like 1/∣k∣ and is not integrable, so part 2 and part 1 of the weaker statement both fail to apply. The symmetric limit of the second weaker statement still works and returns 1 on (−1,1), 0 outside, and 21 at x=±1. This is why Stakgold and Holst (Stakgold & Holst (2011), §2.4, (2.4.5), p. 146) insist that the inversion integral be read as limR→∞∫−RR and not as an absolutely convergent integral.
Continuity of f at x is used once, in the claim of step 2, and it is the only place a property of f at the single point x enters. Drop it and the approximate identity returns the average of the two one-sided limits instead, because Kε is even and so splits its mass equally on the two sides — which is exactly the midpoint value of the second weaker statement.
Boundedness of f is used in the same claim, to control the part of the integral outside ∣y−x∣≤η. For the L1 statement, where f need not be bounded, Folland (Folland (1999), Thm 8.26, pp. 251–252) replaces this step by convergence of f∗Kε→f in the L1 norm, which is why his conclusion is only almost everywhere and names a continuous representative rather than f itself.
The Gaussian is a convenience, not the mechanism. Any positive, integrable, explicitly transformable approximate identity would do. Reed and Simon (Reed & Simon (1980), §IX.1, Thm IX.1, pp. 320–321) prove the same theorem with no Gaussian at all, by expanding a compactly supported f in a Fourier series on a large cube and observing that the series is a Riemann sum for the inversion integral; they attribute the idea to Fourier. Stakgold and Holst (Stakgold & Holst (2011), §2.4, p. 146) use the sinc kernel (R/π)sincRy instead, which is not positive, and pay for it with a more delicate limit.
Part 1 is what makes part 2 unconditional on S. Without it we would have to assume f^∈L1 separately, as the L1 statement does. The chain is: rapid decrease of f gives smoothness of f^, smoothness of f gives rapid decrease of f^, and S is the class where both hold at once. This is the reason Schwartz introduced it, and the reason Fourier transforms of distributions is built on it.
Uniqueness is a corollary, not an input. The L1 statement gives: if f^=g^ with f,g∈L1 then f=g almost everywhere. "Almost everywhere" cannot be improved, because changing f at one point changes no integral.
Method(Inverting a Fourier transform)Original
Decide which theorem applies. If the function at hand is rapidly decreasing, part 2 gives the inverse at every point. If only f,f^∈L1, the inverse integral gives a continuous representative. If f^ decays like 1/∣k∣ or slower, the integral is not absolutely convergent: read it as a symmetric limit and expect the midpoint value at any jump. If f^ does not decay at all, the inverse lives in Fourier transforms of distributions.
Look the pair up backwards. Most inversions in practice are a table entry of The Fourier transform read right to left, combined with the operational rules: a factor e−ika is a translation by a, a factor f^(k−a) is a modulation, and a polynomial factor (ik)m is a derivative.
Use duality rather than integrating again. By part 3, F2=2πR, so a known pair f↔g immediately yields the pair g↔2πf(−⋅). This converts the exponential-to-Lorentzian pair into the Lorentzian-to-exponential pair at no cost.
Differentiate a known pair with respect to a parameter. Suppose fa↔ga for a in an interval, and that on each compact subinterval the difference quotients [fa′(x)−fa(x)]/(a′−a) are bounded in modulus by one fixed integrable function of x. Then dominated convergence (Dominated convergence) lets the a-derivative pass under the integral sign, and ∂afa↔∂aga. This generates whole families from one entry, and is the cheapest route to transforms of repeated poles.
Otherwise use residues. Write eikx and complete the contour through the upper half of the k-plane when x>0, through the lower half when x<0, since ∣eikx∣=e−xImk; see The residue theorem. Poles of f^ in the upper half-plane produce the behaviour for x>0, which is how causality is read off a response function.
Check. Evaluate the recovered f at one convenient point against an independent computation of the same number; confirm ∫∣f∣<∞ if the theorem used requires it; and confirm that real f came from a f^ satisfying f^(−k)=f^(k).
Example(A fourth-power Lorentzian, transformed and recovered)Original
Take
f(x)=(1+x2)21,
which is continuous, positive and O(x−4), so f∈L1 with ∫f=π/2. It is not in S: it is not rapidly decreasing. So part 2 of the main statement does not apply, and we shall need the L1 version.
Step 1: the transform, by differentiating a parameter. Start from the table entry e−a∣x∣(k)=2a/(a2+k2) for a>0. Neither function is in S, so part 2 of the main statement is unavailable; but e−a∣x∣ is continuous and in L1, and its transform 2a/(a2+k2) is in L1 as well, so part 1 of the weaker statement applies and gives, at every x,
e−a∣x∣=2π1∫−∞∞a2+k22aeikxdk.
Read this backwards: put x=−k′ and rename the integration variable k as x, which turns eikx into e−ixk′ and the right-hand side into 2π1 times a transform in our own convention. The result is the dual pair
Both sides are differentiable in a on a>0. Fix 0<a0≤a≤a1. The map a↦(a2+x2)−1 is real-valued with ∣∂a(a2+x2)−1∣=2a/(a2+x2)2≤2a1/(a02+x2)2, so by the mean value theorem every difference quotient of the x-integrand in a is bounded in modulus by the single integrable function 2a1/(a02+x2)2, independent of a. Dominated convergence (Dominated convergence) therefore lets the a-derivative pass under the integral sign. Since ∂a(a2+x2)−1=−2a(a2+x2)−2,
Sanity checks: f^(0)=π/2=∫f, correct; f^ is real and even, as it must be for real even f; and f^ has a kink at k=0, matching the merely algebraic decay of f. Numerically the closed form gives f^(2.3)=0.5197045373799962, and oscillatory quadrature of ∫f(x)e−2.3ixdx returns the same digits (the integrand must be integrated period by period; a naive adaptive rule loses accuracy in the seventh digit here).
Step 2: the inversion theorem applies.f^ decays exponentially, so f^∈L1; and f∈L1. Part 1 of the weaker statement therefore applies, and since f is already continuous, its continuous representative is itself: the inversion formula must return f(x) at every x.
Step 3: carry out the inversion by hand. Because f^ is even, eikx may be replaced by coskx and the range halved:
Step 4: a number. At x=2 the inversion integral must equal f(2)=1/25=0.04 exactly. Numerical evaluation of 2π1∫−∞∞2π(1+∣k∣)e−∣k∣e2ikdk returns 0.0400000000 with zero imaginary part, to the precision of the quadrature.
What the example shows. The Schwartz-class theorem would not have covered this f, but the L1 version does, and the price was checking two integrability conditions rather than one. Had we instead taken f=χ[−1,1], whose transform decays only like 1/∣k∣, neither version would apply and only the symmetric-limit statement would be available.
ApplicationOriginal
The transform method for differential equations. Inversion is what makes the scheme of The Fourier transform a method rather than a one-way map: transform the equation, solve the resulting algebraic problem for u^(k), invert. For the heat equation ut=κuxx on the line with u(x,0)=ϕ(x), the transformed equation is u^t=−κk2u^, so u^(k,t)=ϕ^(k)e−κk2t; inverting the Gaussian factor by step 1 of the proof above gives the heat kernel (4πκt)−1/2e−x2/4κt and, with The convolution theorem, the solution as its convolution with ϕ. See The heat kernel. Note how the hypothesis is met for free: e−κk2t makes u^ rapidly decreasing for every t>0 however rough ϕ was, which is the Fourier-side statement that the heat equation smooths instantly.
Completeness of the plane waves. Read as an identity between operators, inversion says 2π1∫eik(x−y)dk=δ(x−y): the plane waves eikx form a complete set, as the discrete exponentials do for Fourier series. Hassani (Hassani (2013), §9.2, Eqs. (9.26)–(9.28), pp. 281–282) writes this as the resolution of the identity 1=∫∣k⟩⟨k∣dk and uses it to pass between the position and momentum representations of a quantum state; the theorem proved here is what licenses that bookkeeping for wave functions in S, and Fourier transforms of distributions is what licenses it for the plane waves themselves, which are not in S.
Reconstruction from measured spectra. An interferometer or a diffraction experiment measures f^, not f: the Fourier transform spectrometer records the interferogram and the crystallographer records the structure factors. Inversion is the step that turns the measurement into the object, and the hypothesis f^∈L1 is the reason the data must be apodised — truncating the measurement at ∣k∣≤R amounts to multiplying f^ by a square pulse, whose inverse transform is a sinc, so the reconstruction is convolved with a sinc and acquires the ringing that a smooth apodisation window is designed to suppress.
Why the order-four identity matters.F2=2πR says that transforming twice returns the original function reflected. In optics this is the statement that two successive Fourier lenses produce an inverted image; in quantum mechanics it is the statement that the transform generates a quarter-turn in the position-momentum phase plane, whose eigenfunctions are the Hermite functions hn(x)=Hn(x)e−x2/2, with eigenvalues (−i)n2π in the convention used here of Hermite polynomials — the stationary states of the harmonic oscillator, which are therefore exactly the functions that look the same in position and in momentum space.
Defines the convolution and the cross-correlation of two functions on the line, proves that the transform turns convolution into pointwise multiplication and correlation into the product of one transform with the conjugate of the other, and shows how the two results reduce a linear translation-invariant system, a deconvolution and a power spectrum to algebra in the transform variable.
whenever the integral converges, and the cross-correlation is
(f⋆g)(z)=∫−∞∞f(x)g(x+z)dx.
The special case g=f gives the autocorrelationa(z)=(f⋆f)(z)=∫f(x)f(x+z)dx.
The two integrals look alike and mean different things. In a convolution the two functions run in opposite directions, so that f∗g at z collects the contributions of f at x spread by g through the displacement z−x; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.4, Fig. 20.7, p. 986) call this Faltung, folding, because g(z−x) is the graph of g reflected and then slid. In a correlation both run the same way, and f⋆g at z measures how much f and g resemble each other when one is displaced by z; Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.8, p. 449) stress that despite the formal similarity the two are used for entirely different purposes.
Where the convolution comes from.Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.7, pp. 446–447) derive it from measurement. Let f(x)dx be the probability that the true value lies in dx, and let the instrument have resolution functiong, meaning that a true value x is recorded somewhere in [z,z+dz] with probability g(z−x)dz. The joint probability of the two events is f(x)g(z−x)dxdz, and summing over all x that could produce a reading near z gives the observed distribution h=f∗g. Two consequences are immediate from the formula: a perfect instrument, g=δ, gives h=f; and convolving with a sum of delta functions leaves a copy of g at each spike.
Our convention throughout is that of The Fourier transform, f^(k)=∫f(x)e−ikxdx, in which the convolution carries no prefactor. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.4, Eq. (20.68), p. 986) and Reed and Simon (Reed & Simon (1980), §IX.1, p. 323) work in the symmetric convention and must therefore insert 1/2π somewhere: Arfken puts it in the definition of f∗g, Reed and Simon leave the definition bare and pay (2π)n/2 in the theorem. The three choices describe the same mathematics.
It is defined, and it does not leave L1.(f∗g)(z) exists for almost every z, and
∥f∗g∥1≤∥f∥1∥g∥1.
This is the case p=1 of Young's inequality (Folland (1999), Prop. 8.7, p. 240), which gives ∥f∗g∥p≤∥f∥1∥g∥p for every p∈[1,∞].
It is a commutative, associative, distributive product.f∗g=g∗f; (f∗g)∗h=f∗(g∗h); f∗(αg+βh)=α(f∗g)+β(f∗h). It commutes with translations, and Folland (Folland (1999), Prop. 8.6(d), p. 240) adds that supp(f∗g) is contained in the closure of {x+y:x∈suppf,y∈suppg}. There is no identity element in L1: the identity would have to be δ, which is not a function (Test functions and distributions).
Smoothing.f∗g is at least as smooth as the smoother factor: if g is m times differentiable with g,…,g(m) bounded, then f∗g is m times differentiable and (f∗g)(m)=f∗g(m). In particular, if f,g∈S(R) then f∗g∈S(R).
Correlation in terms of convolution. With f~(x)=f(−x) one has f⋆g=f~∗g. Correlation is not commutative; instead
(f⋆g)(z)=(g⋆f)(−z).
Convolution becomes multiplication.f∗g(k)=f^(k)g^(k)for every k∈R.
Correlation becomes a conjugated product (Wiener-Kinchin).f⋆g(k)=f^(k)g^(k),and in particularf⋆f(k)=f^(k)2.
The transform of the autocorrelation is the energy spectrum∣f^∣2, which is real and non-negative whatever f is.
The dual statement: multiplication becomes convolution. If moreover f,g∈S(R), then
fg(k)=2π1(f^∗g^)(k).
The lone factor 2π1 in part 3 is the same 2π that sits on the inverse transform in our convention; parts 1 and 2 are free of constants. In the symmetric convention the constants redistribute, and Reed and Simon (Reed & Simon (1980), §IX.1, Thm IX.3(b), p. 323) state the pair as f∗g=(2π)n/2f^g^ and fg=(2π)−n/2f^∗g^, while Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.7, Eqs. (13.38)–(13.39), p. 448) write 2πf^g^ and (2π)−1/2f^∗g^ in one dimension. Setting z=0 in part 2 and inverting produces the multiplication theorem and Parseval's relation, which are the subject of Plancherel's theorem.
Idea of the proofOriginal
Everything rests on one algebraic accident: e−ikz factorises along the same split as the convolution. Writing z=(z−x)+x gives e−ikz=e−ik(z−x)e−ikx, so the double integral ∬f(x)g(z−x)e−ikzdxdz separates into a product of two single integrals the moment the exponential is split that way and u=z−x is substituted. This is the entire proof; the only real work is checking that the double integral converges absolutely, so that the order of integration may be exchanged.
Correlation works the same way with one sign changed: z=(x+z)−x gives e−ikz=e−ik(x+z)eikx, and the factor e+ikx against f(x) assembles into f^(k) rather than f^(k). The conjugate in the Wiener-Kinchin theorem is therefore not decoration: it records that correlation slides the two functions the same way while convolution slides them opposite ways.
the inner integral being ∥g∥1 for each x by translation invariance. For fixed z, (f∗g)(z) is an integral in x, so the slice needed is the one in that variable, and it is the other order of the same Tonelli computation: ∫(∫∣f(x)∣∣g(z−x)∣dx)dz is the same finite number ∥f∥1∥g∥1, so the inner integral — which is exactly the integral defining (f∗g)(z), with the modulus taken — is finite for almost every z. Hence (f∗g)(z) converges absolutely for almost every z, and the displayed computation bounds ∫∣f∗g∣ by ∥f∥1∥g∥1.
Part 2. For commutativity, which is Folland's (Folland (1999), Prop. 8.6(a), p. 240) Proposition 8.6(a), substitute u=z−x, so x=z−u and dx=−du with the limits reversed:
(f∗g)(z)=∫f(x)g(z−x)dx=∫f(z−u)g(u)du=(g∗f)(z).
Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, p. 143) observe that this substitution is a reflection followed by a translation, both of which preserve the integral. Associativity and distributivity follow the same way, by Fubini and a substitution respectively; Reed and Simon (Reed & Simon (1980), §IX.1, Thm IX.3(c), p. 324) instead deduce both from the transform identity of the next proof together with injectivity of the transform, which is shorter but uses more.
For the absence of an identity: if e∈L1 satisfied e∗f=f for all f∈L1, then transforming with the theorem below gives e^(k)f^(k)=f^(k) for all f and all k. Fix k and take f(x)=e−x2, whose transform πe−k2/4 is nowhere zero; dividing by it gives e^(k)=1, and since k was arbitrary, e^≡1. But the Riemann-Lebesgue lemma of The Fourier transform says e^(k)→0 as ∣k∣→∞, a contradiction.
Part 3. Take m=1 first, so g is differentiable with g and g′ bounded. The z-dependence of (f∗g)(z)=∫f(x)g(z−x)dx sits entirely in g, so for t=0
t(f∗g)(z+t)−(f∗g)(z)=∫f(x)tg(z+t−x)−g(z−x)dx.
For each fixed x the quotient inside tends to g′(z−x) as t→0, by definition of the derivative. For the majorant, apply the mean value theorem to the real and the imaginary part of g separately — it is false for a complex-valued function taken whole — which gives ∣Reg(b)−Reg(a)∣≤∣b−a∣sup∣g′∣ and the same for Im, hence ∣g(b)−g(a)∣≤2∣b−a∣sup∣g′∣. The integrand is therefore dominated in modulus by 2∣f(x)∣sup∣g′∣, which is integrable because f∈L1, and which does not depend on t. Dominated convergence (Dominated convergence) applied along an arbitrary sequence tn→0 gives that the limit exists and equals ∫f(x)g′(z−x)dx; as the sequence was arbitrary, (f∗g)′(z) exists and equals (f∗g′)(z). Iterating gives the mth derivative. Since convolution is commutative, the derivative may be placed on whichever factor is smoother — which is the precise sense in which convolution smooths. For the Schwartz statement, Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, Prop. 1.11(i), p. 142) combine this with the bound supx∣x∣ℓ∣g(x−y)∣≤Aℓ(1+∣y∣)ℓ, which gives supx∣xℓ(f∗g)(x)∣≤Aℓ∫∣f(y)∣(1+∣y∣)ℓdy<∞; rapid decrease of f makes every such integral finite, so f∗g and all its derivatives decay faster than any power.
Part 4. Substituting u=−x in the definition of f~∗g,
which is Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §13.1.8, Eq. (13.41), p. 449) relation (13.41). □
Proof(The convolution theorem)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.4, the derivation of (20.70), pp. 986–987; Folland (1999), Thm 8.22(c), p. 250; Riley, Hobson & Bence (2006), §13.1.7, the derivation of (13.38), p. 448 and §13.1.8, the derivation of (13.42), p. 450
Part 1. By part 1 of the previous statement the double integrand f(x)g(z−x)e−ikz has modulus ∣f(x)∣∣g(z−x)∣, whose integral over R2 is ∥f∥1∥g∥1<∞. Fubini's theorem (Fubini and Tonelli theorems) therefore permits the order of integration to be exchanged, and this is the only analytic step:
In the inner integral substitute u=z−x, so dz=du and e−ikz=e−ikue−ikx; the factor e−ikx does not involve u and comes out:
∫g(z−x)e−ikzdz=e−ikx∫g(u)e−ikudu=e−ikxg^(k).
Since g^(k) no longer depends on x it leaves the outer integral as well, and what remains is g^(k)∫f(x)e−ikxdx=f^(k)g^(k). Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.4, p. 987) describe the same manoeuvre as splitting the exponential into two factors and then changing variable so that each bracket contains only one of them.
Part 2. Identical, with the two changes forced by the definition of ⋆. The double integrand is now f(x)g(x+z)e−ikz, of modulus ∣f(x)∣∣g(x+z)∣, again with double integral ∥f∥1∥g∥1, so Fubini again applies. Substituting u=x+z in the inner integral gives e−ikz=e−ikue+ikx, and
f⋆g(k)=∫f(x)eikxdx∫g(u)e−ikudu=f^(k)g^(k),
the first factor being f^(k) because conjugating f^(k)=∫f(x)e−ikxdx turns the exponent's sign around. Setting g=f gives ∣f^(k)∣2. This is the calculation Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.8, Eq. (13.42), p. 450) present as the Wiener-Kinchin theorem.
Part 3. Let f,g∈S. By Fourier inversion, f^,g^∈S⊂L1, so part 1 may be applied to the pair f^,g^:
using F2=2πR from part 3 of the inversion theorem. Applying F once more to both sides and using F2=2πR again on the left,
2π(f^∗g^)(−k)=(2π)2F[(fg)(−⋅)](k)=(2π)2fg(−k),
the last step by the reflection rule of The Fourier transform. Replacing −k by k and dividing by (2π)2 gives part 3. (A direct proof is also possible, but it needs the inversion theorem anyway to make sense of f^∗g^, so nothing is saved.) □
Remark(Where the hypotheses are used)Original
Integrability of both factors is used exactly once, to license Fubini. It is not a technicality, and what it buys is existence of the convolution at almost every point, never at every point. For f=g=∣x∣−1/2χ[−1,1], which is in L1, the two singularities of ∣x∣−1/2∣z−x∣−1/2 sit at x=0 and x=z and are separately integrable, so (f∗g)(z) converges for every z=0 and diverges only at z=0 — a single point, exactly as part 1 of the algebra statement allows. For f=g=∣x∣−1χ[−1,1], which is not in L1, the two singularities are no longer integrable and (f∗g)(z) diverges for every z with ∣z∣≤1, a set of positive measure. That is the real contrast: without finiteness of the double integral the exceptional set stops being null, and then part 1 has nothing to say.
Part 1 needs nothing about smoothness or decay beyond L1, which is why it is the version Folland (Folland (1999), Thm 8.22(c), p. 249) states. Part 3 is genuinely stronger: a product of two L1 functions need not be in L1 (take f=g=∣x∣−1/2χ[−1,1] again), so fg may not exist, and the right-hand side f^∗g^ may not either, since f^ and g^ are merely bounded and vanishing at infinity. Restricting to S removes both problems at once. Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §1.7, p. 144) note that the pair of results extends to functions of moderate decrease, whose convolution is again of moderate decrease.
The theorem is an identity of functions, not of values at a point. Both sides of part 1 are continuous in k, so the identity does hold everywhere; but f∗g itself is only defined almost everywhere when f,g are merely integrable, and changing it on a null set changes neither side.
The conjugate in part 2 cannot be dropped. If it is, the result is a statement about ∫f(x)g(x+z)dx, which equals (f~r∗g)(z) with f~r(x)=f(−x) and therefore has transform f^(−k)g^(k). For real f the two agree, since f^(−k)=f^(k), which is why the distinction is invisible in most of the engineering literature.
Non-commutativity of correlation is not an oversight. Part 4 of the algebra statement says that exchanging the arguments reflects and conjugates. The autocorrelation is therefore automatically Hermitian, a(−z)=a(z), which is exactly the condition for a^ to be real — consistent with a^=∣f^∣2.
Deconvolution is not the inverse operation in any stable sense. Part 1 says f^=h^/g^, but g^(k)→0, so the division amplifies whatever noise sits in h^ at large ∣k∣ without bound. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.7, p. 449) state the recipe and immediately warn that experimental and statistical uncertainty, and data over a limited range, make it imprecise in practice. This is the standard example of an ill-posed inverse problem.
Recognise the convolution. Any integral of the form ∫K(z−x)u(x)dx over the whole line is a convolution. If the kernel depends on z and x separately and not only through z−x, the problem is not translation-invariant and the method does not apply.
Transform. Part 1 turns the equation into an algebraic one in k. An integral equation u−λK∗u=f becomes u^(1−λK^)=f^; a differential equation becomes a polynomial relation by the derivative rule; and a boundary value problem on a strip or half-plane becomes an ordinary differential equation in the remaining variable.
Solve. Divide: u^(k)=f^(k)/(1−λK^(k)), or whatever the algebra gives. Record where the denominator vanishes; those k are where the method fails and are usually the spectrum of the problem.
Invert, reading products as convolutions. If the answer factorises as u^=G^f^ with G^ a recognisable transform, then u=G∗f without any further integration: this is how Stakgold and Holst (Stakgold & Holst (2011), §6.2, (6.2.31), p. 367 and §7.3, (7.3.49), p. 453) finish both an integral equation and a Dirichlet problem in a strip. The function G is the Green's function or impulse response of the problem; see Causal Green's functions for initial-value problems.
If only a number is wanted, do not invert at all. Evaluating part 1 at k=0 gives ∫(f∗g)=(∫f)(∫g), and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.4, Eq. (20.71), p. 987) note more generally that the inverse transform of the product is often a better-behaved integral than the original, because both factors are then evaluated at the same point.
To deconvolve, regularise. Form f^=h^/g^ only where ∣g^(k)∣ is comfortably above the noise floor, and suppress the rest — for instance by the Wiener filter f^=h^g^/(∣g^∣2+σ2), which reduces to plain division where the signal dominates and to zero where it does not.
Check. The area rule of step 5 is the cheapest test; a second is that the convolution of two non-negative functions is non-negative, and a third is the support rule of part 2 of the algebra statement, so the support of the answer must lie in the sum of the two supports.
Example(Two decaying pulses, convolved two ways)Original
Let
f(t)=H(t)e−t,g(t)=H(t)e−3t,
with H the Heaviside step (The Heaviside step and related functions). Both are in L1(R), with ∥f∥1=1 and ∥g∥1=31. Physically, f is a pulse of unit area with relaxation time 1 driving a system whose impulse response g relaxes three times faster.
Route A: the integral directly. Both factors vanish for negative argument, so in (f∗g)(t)=∫f(s)g(t−s)ds the integrand is non-zero only for 0≤s≤t, and the convolution vanishes for t<0. For t≥0,
Route B: through the transform. From the table of The Fourier transform, f^(k)=1/(1+ik) and g^(k)=1/(3+ik), so by part 1 of the theorem
f∗g(k)=(1+ik)(3+ik)1=21[1+ik1−3+ik1],
splitting by partial fractions in the variable ik (Partial fractions for inverse transforms). Each term is a tabulated transform, so reading backwards gives (f∗g)(t)=21(e−t−e−3t)H(t), the same answer, with no integration performed. Route B is the one that generalises: the partial-fraction step works for any product of rational transforms.
A number. At t=ln2, where e−t=21 and e−3t=81,
(f∗g)(ln2)=221−81=163=0.1875.
Direct numerical quadrature of ∫0ln2e−se−3(ln2−s)ds gives 0.1875000000, and numerical evaluation of f∗g at k=1.3 from the closed form returns 0.0455555517−0.1808311976i, against 1/((1+1.3i)(3+1.3i))=0.0455555517−0.1808311976i.
The checks of the method. The area rule gives ∫(f∗g)=1⋅31=31, and indeed ∫0∞21(e−t−e−3t)dt=21(1−31)=31. The support rule gives supp(f∗g)⊆[0,∞)+[0,∞)=[0,∞), as found. And the peak of f∗g, at t=21ln3≈0.549, lies to the right of the peaks of both factors at t=0: the convolution has delayed and broadened the pulse, which is what a resolution function does.
An autocorrelation, for contrast. Take instead p=χ[−a,a]. Its autocorrelation is the overlap of the pulse with its own translate, ap(z)=(2a−∣z∣)+, a triangle of height 2a and half-width 2a. Part 2 predicts a^p=∣p^∣2=4sin2(ak)/k2, and the triangle entry of the table, applied with half-width 2a and height 2a, gives 2a⋅2ak24sin2(ak)=k24sin2(ak) — the same function. Notice that the autocorrelation is even and its transform is non-negative, as part 2 requires, whereas p^ itself takes both signs.
ApplicationOriginal
Linear time-invariant systems. A system that is linear and commutes with time translation has output y=h∗u for a fixed impulse response h, since u=∫u(s)δ(⋅−s)ds and linearity plus translation invariance send δ(⋅−s) to h(⋅−s). The theorem then says y^=h^u^: in the frequency domain the system acts by multiplication, one frequency at a time, with no mixing. h^(ω) is the transfer function, its modulus the gain and its argument the phase shift, and this is the whole content of circuit analysis by impedances, of filter design, and of the Bode plot. Causality — that h vanishes for negative time — is what forces the real and imaginary parts of h^ to be related, by the Kramers-Kronig relations of Dispersion relations.
Instrumental resolution and deconvolution. The derivation in the definition is the general model of a measurement: every spectrometer, telescope and detector returns the true distribution convolved with its point-spread or resolution function. Two consequences follow immediately from the theorem. First, resolutions compose by convolution and therefore widths add in quadrature for Gaussians, since the product of two Gaussians in k is a Gaussian. Second, information at spatial frequencies where g^ is negligible is lost and cannot be recovered by any amount of processing — the diffraction limit of an optical instrument is the statement that g^ is supported in ∣k∣≤2π/λ times a numerical aperture factor.
Power spectra of noise. For a stationary random signal the autocorrelation a(z)=⟨f(t)f(t+z)⟩ is the natural object, and part 2 identifies its transform as the power spectral density. It gives two routes to the same quantity: transform an estimate of the autocorrelation, or square the modulus of a single transform of the data. The second is the periodogram of The periodogram and spectral estimation, and part 2 is what says the two estimate the same thing; they differ only in how the estimate is smoothed, which is why the raw periodogram is noisy and why windowing and averaging are applied to it. The same identity run backwards generates coloured noise with a prescribed spectrum: draw white noise, multiply its transform by a^, and invert.
Green's functions and potentials. For a constant-coefficient linear differential operator L on the line, Lu=f transforms to P(ik)u^=f^ with P the symbol of L, so u^=f^/P(ik) and u=G∗f with G^=1/P(ik): the Green's function is the inverse transform of the reciprocal symbol, and the solution is the source convolved with it (Causal Green's functions for initial-value problems). Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Example 20.4.1, (20.79)–(20.80), pp. 988–989) run exactly this in three dimensions for electrostatics, using the R3 form of part 1 and the three-dimensional pair of Fourier transform in several dimensions: the Coulomb integral ψ=4π1∫ρ(r′)/∣r−r′∣d3r′ is a convolution of ρ with 1/r, whose transform is 4π/k2, so ψ^(k)=ρ^(k)/k2 — Poisson's equation in one line.
Matched filtering. To detect a known waveform s buried in noise, correlate the data with s rather than convolve: part 2 says the correlation has transform s^d^, which weights each frequency by the amplitude the signal has there, and this is the filter that maximises the signal-to-noise ratio at the detection instant for white noise. Gravitational-wave searches and radar ranging are the same calculation.
Proves that the Fourier transform preserves the L2 inner product up to the factor 2π fixed by the convention, by transforming the cross-correlation and evaluating the inversion formula at the origin, then extends the transform from the Schwartz class to a unitary map of L2 onto itself and uses the resulting Parseval relation to evaluate integrals in closed form.
Let L2(R) be the space of measurable f with ∫∣f∣2<∞, equipped with
(f,g)=∫−∞∞f(x)g(x)dx,∥f∥2=(f,f)1/2,
linear in the first slot, as in L^p spaces. It is a Hilbert space (Hilbert spaces), and the Schwartz class S(R) of Fourier inversion is a dense subspace of it: the continuous functions of compact support are dense in Lp for p<∞ (L^p spaces), and smoothing one of them by convolution with a bump puts it in Cc∞⊂S without moving it far in ∥⋅∥2, which is Folland's (Folland (1999), §8.2, Prop. 8.17, p. 245) Proposition 8.17.
Neither of L1(R) and L2(R) contains the other, which is the one structural fact the theorem has to accommodate.
f(x)=(1+∣x∣)−1 is in L2 but not in L1, so the defining integral ∫f(x)e−ikxdx of The Fourier transform need not converge for an L2 function.
f(x)=∣x∣−1/2χ[−1,1](x) is in L1 but not in L2, so an integrable function need not have a square-integrable transform.
The overlap L1∩L2 contains S and is dense in L2. On it the transform is given by the integral, and the theorem below says it is an isometry there up to a constant; the transform of a general L2 function is then defined by continuous extension. Reed and Simon (Reed & Simon (1980), §IX.2, p. 328) give the resulting integral representation as a limit in the L2 norm. They write it in their symmetric convention, as l.i.m.R→∞(2π)−n/2∫∣x∣≤Re−ik⋅xf(x)dx; taking n=1, so that their prefactor is (2π)−1/2, and then dropping that prefactor, which our F does not carry, it becomes
f^(k)=R→∞l.i.m.∫−RRf(x)e−ikxdx,
where l.i.m. ("limit in the mean") means convergence of the truncated integrals in ∥⋅∥2, not at any particular k. The rescaling is the only change: our f^ is 2π times theirs, which is also why their unitary operator is F itself and ours is (2π)−1/2F.
The extension to L2. Both identities hold for every f,g∈L1∩L2, and the operator
U=2π1F
extends uniquely from L1∩L2 to a unitary operator of L2(R)ontoL2(R), with U−1=U† given by the correspondingly normalised inverse transform. In particular ∥f∥22=2π1∥f^∥22 for every f∈L2, with the transform read as the limit in the mean of the definition above.
Parts 1 and 2 are stated in exactly this form, with the 1/2π in the same place, by Strauss (Strauss (2008), §12.3, (14)–(15), p. 346), who adds that if either integral is finite so is the other. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.4, (20.75)–(20.77), p. 988) write part 1 as ⟨f∣g⟩=⟨Ff∣Fg⟩ and read off F†=F−1, which is part 3 in operator language. The constant depends on the convention, and only on it: it is 1 in the symmetric conventions of Reed and Simon (Reed & Simon (1980), §IX.2, Thm IX.6, p. 327) and Stakgold and Holst's (Stakgold & Holst (2011), §10.3, Thm 10.3.2, p. 701) Chapter 10, and 1 again in the 2π-in-the-exponent convention of Folland (Folland (1999), Thm 8.29, p. 252); see the conventions table of The Fourier transform.
Stakgold and Holst (Stakgold & Holst (2011), §2.4, (2.4.7)–(2.4.8), p. 147) give two formulas under the single name Parseval's formula, and they are not the same statement:
Formula (a) is Plancherel's theorem, quadratic in f and carrying a 2π; (b) is the multiplication formula, bilinear, carrying no constant, with no conjugation, and holding for any f,g∈L1 — it was proved as step 3 of the inversion proof in Fourier inversion, and it is Folland's (Folland (1999), Lemma 8.25, p. 251) Lemma 8.25. The two coincide only after inversion is used to replace one ^ by its inverse. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.9, (13.45), p. 451) reserve the name multiplication theorem for a third form again, conjugated like (a) but bilinear like (b). They write it constant-free, as ∫f∗gdx=∫[f~]∗g~dk, because their transform is the symmetric f~=(2π)−1/2f^; the two factors of (2π)−1/2 on the right supply the constant. In our hat the same statement reads
∫f∗(x)g(x)dx=2π1∫[f^(k)]∗g^(k)dk,
which is part 1 above. We keep Parseval's relation for part 1, Plancherel's theorem for parts 2 and 3, and multiplication formula for (b). Stakgold and Holst note that their derivation of (a) and (b) is a formal interchange of orders of integration and that they do not state precise conditions at that point.
Idea of the proofOriginal
A quadratic identity in f is a linear identity in the autocorrelation of f. By the Wiener-Kinchin half of The convolution theorem, the cross-correlation h=f⋆g has transform f^g^ — the product whose integral we want on the right-hand side. The integral we want on the left-hand side is h evaluated at displacement zero, since h(0)=∫fg. And the inversion formula of Fourier inversion, evaluated at the single point 0, says precisely that h(0) is 2π1 times the integral of h^. Putting the three sentences together is the whole proof: no new estimate is needed, only the two theorems already proved and the decision to look at one point of one auxiliary function.
The passage from the Schwartz class to L2 is then abstract. An isometry defined on a dense subspace extends to the whole space in exactly one continuous way, and the extension is still an isometry; it is onto because the transform already maps the dense subspace onto itself.
Step 1: h lies in S. By part 4 of the algebra statement of The convolution theorem, h=g~∗f with g~(x)=g(−x). Reflection and conjugation both preserve S, so g~∈S, and part 3 of the same statement gives g~∗f∈S. In particular h is continuous, so the value h(0) is unambiguous, and h^∈S⊂L1, so the inversion theorem applies to h at every point.
Step 2: the two sides, read off. The value at the origin is, directly from the definition of h,
h(0)=∫−∞∞g(x)f(x)dx=(f,g).
The transform is, by the Wiener-Kinchin part of The convolution theorem applied to the pair (g,f) in that order,
h^(k)=g^(k)f^(k).
Step 3: inversion at one point. Part 2 of Fourier inversion, applied to h∈S and evaluated at z=0, where eik⋅0=1, gives
h(0)=2π1∫−∞∞h^(k)dk.
Substituting the two expressions from step 2,
∫−∞∞f(x)g(x)dx=2π1∫−∞∞f^(k)g^(k)dk,
which is part 1. Setting g=f makes both integrands non-negative and gives part 2. The route is Stein and Shakarchi's (Stein & Shakarchi (2003), Ch. 5, Thm 1.12, pp. 143–144), who take g=f from the start and write the auxiliary function as f∗f♭ with f♭(x)=f(−x), the same object as our g⋆f with g=f; carrying two functions costs nothing and delivers the polarised identity without a separate argument. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.9, (13.44)–(13.45), p. 451) and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.4, (20.73)–(20.74), p. 987) organise the same calculation as setting the displacement variable to zero in the correlation and convolution theorems respectively. □
Proof(Part 3, the extension to L2)AdaptedAdapted from Folland (1999), §8.2, Prop. 8.17, p. 245 and §8.3, Thm 8.29, pp. 252–253; Reed & Simon (1980), §IX.2, Thm IX.6, p. 327
Write U=(2π)−1/2F, so that part 2 reads ∥Uf∥2=∥f∥2 for f∈S.
Step 1: U is an isometry on a dense subspace.S is dense in L2(R), by the density statement recorded in the definition above (Folland (Folland (1999), §8.2, Prop. 8.17, p. 245)), and U maps it into L2 isometrically by part 2.
Step 2: extend. A linear isometry defined on a dense subspace of a normed space, with values in a complete space, has exactly one continuous extension to the whole space, obtained by taking limits of Cauchy sequences; the extension is again linear and isometric, since both properties pass to limits. Apply this to U on S⊂L2, with L2 complete. Call the extension U again.
Step 3: U is onto. By part 3 of Fourier inversion, F(S)=S, so the range of U contains the dense set S. The range of an isometry defined on a complete space is closed: if Ufn→v then (Ufn) is Cauchy, hence so is (fn) by isometry, hence fn→f for some f∈L2 and v=Uf by continuity. A closed set containing a dense set is everything, so U is surjective. A surjective linear isometry of a Hilbert space is unitary, with U−1=U†; see Hilbert spaces.
Step 4: the extension agrees with the integral on L1∩L2. This is the one point that is not formal, since a function in L1∩L2 need not be in S and the extension was built by limits. Folland (Folland (1999), Thm 8.29, pp. 252–253) settles it by convolving with the Gaussian gt of the inversion proof: for f∈L1∩L2, f∗gt lies in the set X={f∈L1:f^∈L1}, on which the integral formula is itself an isometry, since for u∈X the inversion theorem and the multiplication formula give ∥u∥22=2π1∥u^∥22 directly. On X the integral formula therefore agrees with the extension built in step 2, and the reason is worth spelling out: X is a linear subspace of L2 containing S, and (2π)−1/2F is isometric on it, hence ∥⋅∥2-continuous there; so if u∈X and φn∈S with φn→u in L2, then (2π)−1/2φ^n→(2π)−1/2u^ in L2 by that isometry applied to the differences u−φn∈X, while (2π)−1/2φ^n=Uφn→Uu by construction of the extension. Finally f∗gt→f in L1 and in L2 as t→0, so f∗gt→f^ uniformly and in L2; the two limits must coincide. Reed and Simon (Reed & Simon (1980), §IX.2, p. 328) run the matching argument on the other side, truncating an L2 function to ∣x∣≤R to land in L1 and passing to the limit in the mean, which is where their l.i.m. representation comes from.
Step 5: from norms to inner products. Steps 2 and 3 deliver only that U preserves the norm and is onto. The inner-product form, part 1 for arbitrary f,g∈L2, follows by polarisation: in a complex inner-product space
(u,v)=41m=0∑3imu+imv22,
so a linear map preserving every norm preserves every inner product. Hence (Uf,Ug)=(f,g) for all f,g∈L2, which on multiplying by 2π is (f,g)=2π1(f^,g^). □
Remark(Where the hypotheses are used)Original
Parts 1 and 2 use rapid decrease only to guarantee that the three ingredients are legitimate: that h=g⋆f exists and is continuous at 0, that h^∈L1 so that inversion applies, and that both sides are finite. Any hypothesis delivering those three would do, and Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §1.7, p. 144) note that moderate decrease of f and of f^ suffices.
Continuity of h at the origin is what makes "evaluate at one point" meaningful. Both sides of part 1 are integrals, so neither changes if f is altered on a null set; the proof is consistent with that only because h, being a correlation of two L2 functions, is genuinely continuous and so has an unambiguous value at 0.
Square-integrability cannot be dropped and integrability cannot replace it. For f=∣x∣−1/2χ[−1,1]∈L1∖L2, part 2 would assert that 2π1∫∣f^∣2 equals a divergent integral; what is true instead is the Hausdorff-Young inequality, which Folland (Folland (1999), Thm 8.30, p. 253) states in his convention as ∥f^∥q≤∥f∥p for 1≤p≤2 and q conjugate to p, interpolating between the L1 bound and Plancherel.
Surjectivity is a separate fact from isometry, and it is what "unitary" adds. It comes from inversion, in step 3. The contrast with the L1 theory is sharp: Reed and Simon (Reed & Simon (1980), §IX.2, p. 327) observe that F maps L1 into but not onto the continuous functions vanishing at infinity, so there is no inverse theorem of this kind on L1.
For f∈L2∖L1 the symbol f^(k) denotes an equivalence class, not a number. The truncated integrals converge in ∥⋅∥2, which forces convergence at no individual k; this is why Reed and Simon (Reed & Simon (1980), §IX.2, p. 328) write l.i.m. rather than lim, and why statements about the pointwise behaviour of an L2 transform need separate work.
The constant 2π is a fact about the convention, not about the transform. Our F is not itself unitary; (2π)−1/2F is. In the symmetric convention the normalisation is built into the definition and the theorem reads ∥f∥2=∥f~∥2, which is the form in which Stakgold and Holst (Stakgold & Holst (2011), §10.3, Thm 10.3.2, p. 701) state it. A numerical answer obtained from a book in another convention must be rescaled before it is compared with one obtained here.
Method(Evaluating an integral by Parseval)Original
Match the integrand to a squared transform. The method applies when the integral has the shape ∫∣F(k)∣2dk or ∫F(k)G(k)dk for functions F,G that are recognisable transforms. Powers of sin over powers of the variable, and products of two Lorentzians, are the two common signatures.
Find the originals. Read the table of The Fourier transform backwards to get f with f^=F, using the operational rules and, where F is itself a product, The convolution theorem to write f as a convolution. A factor sin2 usually means f is a convolution of two pulses, hence piecewise linear.
Rescale the variable before anything else if the transform appears with an argument other than k; a substitution k=λu multiplies the answer by λ and is the commonest source of error.
Compute the left-hand side.∫∣f∣2 is almost always elementary, because f is a pulse, a triangle or an exponential even when F is not elementary at all. That asymmetry is the point of the method.
Equate, keeping the 1/2π on the transform side, and solve for the integral wanted.
Check. Both sides of part 2 are non-negative, so a negative answer is wrong; both are unchanged by translating f, so the answer may not depend on any phase; and a numerical quadrature of the original integral should reproduce the closed form to several digits. For a one-sided integral of an even integrand, halve at the end rather than at the start.
Example(Two integrals that are hard directly)Original
First: a fourth power of sinc. Let
T(x)=(1−∣x∣)+,
the triangle of height 1 and half-width 1. By The convolution theorem, T=χ[−1/2,1/2]∗χ[−1/2,1/2], since the convolution of a pulse with itself is the overlap of the pulse with its translate; and χ[−1/2,1/2](k)=2sin(k/2)/k, so
T^(k)=(k2sin(k/2))2=k24sin2(k/2),
in agreement with the triangle entry of the table of The Fourier transform at half-width 1. T is continuous with compact support, so T∈L1∩L2 — but it is not in S, being non-differentiable at x=0,±1. Part 2 as stated therefore does not cover it, and we use part 3, which extends the same identity to all of L1∩L2. The left-hand side is elementary:
∫−∞∞T(x)2dx=2∫01(1−x)2dx=32.
So
32=2π1∫−∞∞k416sin4(k/2)dk.
Substitute k=2u, so dk=2du and k4=16u4; the two factors of 16 cancel and the 2 from dk survives:
32=2π1∫−∞∞16u416sin4u2du=π1∫−∞∞u4sin4udu,
so
∫−∞∞u4sin4udu=32π=2.0943951024…
Oscillatory quadrature over (−∞,∞), taking the interval [0,π] separately and summing the tail over successive half-periods of sinu, returns 2.0943951024, confirming the closed form to ten digits. Nothing in the calculation required integrating sin4u/u4.
Second: a product of two Lorentzians, by the two-function form. Let f(x)=e−2∣x∣ and g(x)=e−5∣x∣. Neither is in S — both have a corner at the origin — but both are continuous and decay exponentially, so both lie in L1∩L2 and part 3 again supplies the identity of part 1 for them. They are real, so the conjugations may be dropped on the left, and from the table f^(k)=4/(4+k2) and g^(k)=10/(25+k2), both real. The left-hand side is
∫−∞∞e−2∣x∣e−5∣x∣dx=∫−∞∞e−7∣x∣dx=72,
and the right-hand side is 2π1∫(4+k2)(25+k2)40dk=π20∫(4+k2)(25+k2)dk. Equating,
∫−∞∞(4+k2)(25+k2)dk=72⋅20π=70π=0.0448798950…
This one can be checked in closed form independently: partial fractions give (4+k2)(25+k2)1=211[4+k21−25+k21], whose integral is 211(2π−5π)=70π. Numerical quadrature agrees to the full working precision, 0.0448798951.
What the two have in common. In both cases the integral asked for lives on the transform side, where it is awkward, and Parseval moves it to the function side, where it is a polynomial or an exponential. The second example also shows the converse reading: given the elementary partial-fraction evaluation, the same identity run backwards computes ∫e−7∣x∣dx, which is a consistency check on the pair of transforms rather than a new result.
ApplicationOriginal
Conservation of energy through an aperture. In Fraunhofer diffraction the far-field amplitude is the Fourier transform of the aperture function, and the observed intensity is its squared modulus. Plancherel then says that the total energy in the diffraction pattern equals the total energy transmitted by the aperture: Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.4, p. 988) point out that in this problem the Parseval relation is energy conservation. The same statement in a waveguide or an optical fibre says that redistributing power among modes cannot change the total.
Energy spectra of transient signals. For a decaying oscillation, ∣f(t)∣2 is the instantaneous energy and ∣f^(ω)∣2 the energy per unit frequency interval; Plancherel says the two accountings agree. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.9, p. 451) make this explicit for the damped harmonic oscillator f(t)=H(t)e−t/τsinω0t, whose transform has Lorentzian peaks of half-width 1/τ at ±ω0: the total energy may be computed by integrating over time or over frequency, and the line width observed in the spectrum is the inverse of the decay time measured in the signal. Combined with the Wiener-Kinchin theorem of The convolution theorem, this is the foundation of spectral density estimation in The periodogram and spectral estimation.
Quantum mechanics: position and momentum are the same state. If ψ∈L2(R) is normalised, ∥ψ∥2=1, then part 3 says the momentum-space wave function ψ^/2π is normalised too. So the total probability is the same whether the particle is described by where it is or by how fast it is going, and no normalisation constant has to be rechosen when changing representation. Unitarity is the statement that the change of representation is a rotation of the Hilbert space; it is why expectation values, overlaps and transition probabilities may all be computed in whichever of the two is more convenient, and it is the one-particle case of the framework of The mathematical framework of quantum mechanics.
The uncertainty principle and energy estimates. Plancherel is the step that converts a statement about ∫∣f′∣2 into one about ∫k2∣f^∣2, because f′=ikf^. Strauss (Strauss (2008), §12.3, (16)–(18), pp. 346–347) uses exactly this to prove the uncertainty inequality ΔxΔk≥21 from the Cauchy-Schwarz inequality, and the same move turns a Sobolev norm into a weighted L2 norm of the transform, which is how Stakgold and Holst (Stakgold & Holst (2011), §10.3, pp. 701–702) define the fractional spaces Hs(Rn) by the finiteness of ∥(1+∣⋅∣2)s/2u^∥2. See The uncertainty principle and Weak derivatives and Sobolev spaces.
Numerical spectral methods. The discrete analogue of part 2 is what makes a spectral method's error measurable: the squaredL2 error of a truncated expansion is exactly the tail 2π1∫∣k∣>K∣f^∣2dk, with no constant to estimate, so the decay of the computed spectrum is a direct read-out of the accuracy. See Spectral methods and, for the discrete version, The discrete Fourier transform and the FFT.
A normalised function and its Fourier transform cannot both be sharply concentrated: the product of their two variances is at least 1/16π2, with equality exactly for Gaussians.
Definition(Dispersion about a point, and the convention used here)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, §4, (11)–(13), pp. 159–160
Throughout this unit the Fourier transform of ψ∈S(R), the Schwartz space of smooth functions all of whose derivatives decay faster than every power of 1/∣x∣, is taken in the form
ψ^(ξ)=∫−∞∞ψ(x)e−2πixξdx,
the convention of Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §1, p. 134). In it the derivative rule reads ψ′(ξ)=2πiξψ^(ξ) and Plancherel's theorem reads ∫∣ψ∣2=∫∣ψ^∣2, with no stray constants in either. The other conventions in use change the numerical constant in the theorem but not its content; see the note below and the Fourier transform.
Call ψnormalised if ∫−∞∞∣ψ(x)∣2dx=1. Then ∣ψ(x)∣2dx and, by Plancherel, ∣ψ^(ξ)∣2dξ are both probability densities, and for x0,ξ0∈R we write
for the two dispersions, or variances. Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §4, (13), p. 160) call Vx the uncertainty attached to the expected position x=∫x∣ψ∣2dx, and record that x0=x is the choice that makes Vx(ψ;x0) smallest, since differentiating in x0 gives ∫(x−x0)∣ψ∣2=0. A small Vx means that nearly all of the mass of ∣ψ∣2 sits close to x0; a large one means it is spread out.
Theorem(Heisenberg inequality)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, Thm 4.1, p. 158
Let ψ∈S(R) satisfy ∫−∞∞∣ψ(x)∣2dx=1. Then
(∫−∞∞x2∣ψ(x)∣2dx)(∫−∞∞ξ2∣ψ^(ξ)∣2dξ)≥16π21,
with equality attained exactly when ψ(x)=Ae−Bx2, where B>0 and ∣A∣2=2B/π. More generally, for every x0,ξ0∈R,
Vx(ψ;x0)Vξ(ψ;ξ0)≥16π21.Note(The same inequality in other conventions)AdaptedAdapted from Strauss (2008), §12.3, (16), p. 346; Riley, Hobson & Bence (2006), §13.1.1, pp. 436–437, and §19.2, pp. 665–666; Kreyszig (1978), §11.2, Thm 11.2-2, p. 580
The constant 1/16π2 belongs to the convention above and to nothing else. Three other statements of the same fact appear in the sources, and they are reconciled here so that no two of them are mixed later.
Angular frequency.Strauss (Strauss (2008), §12.3, (16), p. 346) writes F(k)=∫f(x)e−ikxdx, so that Parseval reads ∫∣f∣2dx=∫∣F∣2dk/2π, sets Δx2=∫∣xf(x)∣2dx and Δk2=∫∣kF(k)∣2dk/2π for normalised f, and states the principle as Δx⋅Δk≥21. Since k=2πξ and ∣F(2πξ)∣2dk/2π=∣f^(ξ)∣2dξ, his Δk2 is 4π2 times our second integral, and ΔxΔk≥21 is the inequality above word for word.
Gaussian widths.Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.1, pp. 436–437) work with the symmetric convention and the normalised Gaussian f(t)=(τ2π)−1e−t2/2τ2, whose transform is (2π)−1/2e−τ2ω2/2; the root-mean-square widths of f and f~ satisfy ΔωΔt=1 whatever τ is, while the widths of the intensities∣f∣2 and ∣f~∣2 are τ/2 and 1/(2τ). With p=ℏk and E=ℏω they obtain ΔEΔt=ℏ/2 and ΔpΔx=ℏ/2, and note that the factor 21 is special to the Gaussian and is the least value attainable.
Operators.Kreyszig (Kreyszig (1978), §11.2, Thm 11.2-2 with (11), p. 580) takes Q to be multiplication by the coordinate q and Dψ=(h/2πi)dψ/dq, and writes the principle as sdψ(D)sdψ(Q)≥h/4π, that is ≥ℏ/2, which agrees with the previous item. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §19.2, (19.38), p. 665) give the general two-observable form (ΔA)2(ΔB)2≥41∣⟨ψ∣[A,B]∣ψ⟩∣2, of which ΔpxΔx≥21ℏ is the case [px,x]=−iℏ.
A difference between the sources.Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §4, p. 161), having rescaled units to remove Planck's constant, restore it by writing the physical conclusion as (uncertainty of position)×(uncertainty of momentum)≥ℏ/16π2, their "uncertainty" being the variance (13) rather than its square root. Kreyszig (Kreyszig (1978), §11.2, (11), p. 580) and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §19.2, p. 666) write ℏ/2 for the product of the two standard deviations. The two statements are not the same inequality, and the physical form used in the application below is the one with ℏ/2, taken from Kreyszig and from Riley, Hobson and Bence. The mathematical theorem proved here is Stein and Shakarchi's and is stated exactly as they give it.
Idea of the proofOriginal
A normalised ψ carries unit mass, and the whole of that mass can be produced by one integration by parts: 1=−∫xdxd∣ψ∣2dx. The derivative on the right brings in ψ′, so the number 1 has been written as a pairing of xψ against ψ′. Cauchy-Schwarz then says that this pairing cannot exceed the product of the two lengths ∥xψ∥ and ∥ψ′∥, and neither length can be small if their product must beat 1. The last move translates the second length into the frequency side: differentiation is multiplication by 2πiξ on the transform, so ∥ψ′∥=2π∥ξψ^∥ by Plancherel. That is the inequality. Equality demands equality in Cauchy-Schwarz, which forces ψ′ to be a constant multiple of xψ, and the only decaying solutions of that differential equation are the Gaussians.
Step 1: write the normalisation as a boundary-free integration by parts. Take u(x)=x and v(x)=∣ψ(x)∣2. Because ψ is a Schwartz function, x∣ψ(x)∣2→0 as ∣x∣→∞, so the boundary term vanishes and
1=∫−∞∞∣ψ(x)∣2dx=−∫−∞∞xdxd∣ψ(x)∣2dx.
This is the only place the normalisation is used, and the only place rapid decay is needed.
Step 2: differentiate the modulus squared. From ∣ψ∣2=ψψ and the product rule,
so step 1 reads 1=−2∫xRe(ψ′ψ)dx. A real number is at most its own modulus, and ∣Rez∣≤∣z∣, hence
1≤2∫−∞∞∣x∣∣ψ(x)∣∣ψ′(x)∣dx.
Step 3: Cauchy-Schwarz. Pair the factor ∣xψ∣ against ∣ψ′∣ in L2(R):
1≤2(∫−∞∞x2∣ψ(x)∣2dx)1/2(∫−∞∞∣ψ′(x)∣2dx)1/2.
Both integrals are finite because ψ∈S(R).
Step 4: move the second factor to the frequency side. Since ψ is Schwartz, so is ψ′, and the derivative rule gives ψ′(ξ)=2πiξψ^(ξ). Applying Plancherel's theorem to ψ′,
Step 5: combine. Substituting step 4 into step 3 and squaring both sides of the resulting inequality 1≤4π∥xψ∥∥ξψ^∥ gives
1≤16π2(∫x2∣ψ∣2dx)(∫ξ2∣ψ^∣2dξ),
which is the first inequality.
Step 6: arbitrary centres. Put φ(x)=e−2πixξ0ψ(x+x0). Translation and modulation are the two rules ψ(⋅+x0)(ξ)=e2πix0ξψ^(ξ) and e−2πixξ0g(ξ)=g^(ξ+ξ0), so φ^(ξ)=e2πix0(ξ+ξ0)ψ^(ξ+ξ0). Both operations leave moduli alone, so φ is again a normalised Schwartz function, ∫x2∣φ∣2dx=Vx(ψ;x0) and ∫ξ2∣φ^∣2dξ=Vξ(ψ;ξ0) after the substitutions x↦x−x0 and ξ↦ξ−ξ0. The first inequality applied to φ is the second.
Step 7: when equality holds. Suppose the product equals 1/16π2. Then step 3 must be an equality, so by the equality case of Cauchy-Schwarz the two functions xψ and ψ′ are proportional: there is a constant β with
ψ′(x)=βxψ(x)for all x.
Separating variables, (logψ)′=βx wherever ψ=0, so ψ(x)=Aeβx2/2 for a constant A, and this formula also covers the case A=0, excluded by normalisation. For ψ to be a Schwartz function β must be real and negative; writing β=−2B with B>0 gives ψ(x)=Ae−Bx2. Finally the normalisation fixes ∣A∣: by the Gaussian integral∫e−2Bx2dx=π/2B, so 1=∣A∣2π/2B, that is ∣A∣2=2B/π.
Conversely let ψ(x)=Ae−Bx2 with these constants. Then ψ′=−2Bxψ, so ∫∣ψ′∣2=4B2∫x2∣ψ∣2. The Gaussian moment ∫x2e−2Bx2dx=4B1π/2B gives ∫x2∣ψ∣2dx=1/4B, hence ∫∣ψ′∣2=B and, by step 4, ∫ξ2∣ψ^∣2dξ=B/4π2. The product is (1/4B)(B/4π2)=1/16π2 for every B, so every such Gaussian attains the bound. □
Remark(Where the hypotheses are used)Original
Rapid decay is used once, in step 1, to discard the boundary term [x∣ψ∣2]−∞∞. It is not a technicality: for a normalised ψ with ∣ψ(x)∣2∼c/∣x∣3 the boundary term still vanishes, but for slower decay the first integral ∫x2∣ψ∣2 is infinite and the inequality says nothing. The theorem is therefore vacuous exactly when it would be hardest to prove, which is why the Schwartz class costs nothing here. Strauss (Strauss (2008), §12.3, (19), p. 347) runs the same step for a test function and gets the same boundary-free value −21 for ∫xff′dx when f is real.
Normalisation sets the constant and nothing else. For a general ψ with ∥ψ∥2=N, running the proof on ψ/N gives (∫x2∣ψ∣2)(∫ξ2∣ψ^∣2)≥N2/16π2: the bound scales with the square of the mass, as it must, since both sides are quadratic in ∣ψ∣2.
Dilation invariance is the reason a bound on the product is the only bound available. If ψδ(x)=δ−1/2ψ(x/δ) then ψδ is still normalised, its position variance is δ2 times that of ψ and its frequency variance is δ−2 times that of ψ. So each factor separately can be made as small as we please, and only the product is constrained. Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, Cor. 1.5 and the following paragraph, p. 139) make the same observation for the Gaussian kernels Kδ(x)=δ−1/2e−πx2/δ: as δ→0 the kernel spikes and its transform flattens.
Equality needs two things, not one. The proof above follows Stein and Shakarchi in reading the equality case off Cauchy-Schwarz in step 3. Equality also has to hold in step 2, where Re(ψ′ψ) was replaced by ∣ψ′∣∣ψ∣ (up to sign). That costs nothing: for ψ′=βxψ with β real, ψ′ψ=βx∣ψ∣2 is already real and of constant sign against x, so step 2 is an equality automatically. Had β been complex, step 2 would be strict, which is a second way of seeing why the minimiser cannot carry a chirp eiγx2.
The two integrals must be computed about the right centres. The constant-centre version with x0=ξ0=0 is weaker, not stronger, than it looks: for a Gaussian shifted to sit at x=10 the integral ∫x2∣ψ∣2 is about 100, so the product beats 1/16π2 by a wide margin without saying anything about the width. Step 6 is what makes the inequality say something about spreads rather than about distances from the origin.
Nothing here is special to one dimension or to this pair of variables. The operator form of Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §19.2, (19.38), p. 665) shows that the obstruction is the commutator: two observables that commute carry no joint restriction, and if [A,B] is a non-zero constant then no state makes both spreads vanish.
Method(Testing a given profile against the bound)Original
Normalise: replace ψ by ψ/∥ψ∥L2, so that ∫∣ψ∣2=1.
Find the two centres x=∫x∣ψ∣2dx and ξ=∫ξ∣ψ^∣2dξ. For real even or real odd ψ both are 0 by symmetry, and the step can be skipped.
Compute the position variance Vx=∫(x−x)2∣ψ∣2dx directly.
Compute the frequency variance without computing ψ^, using step 4 of the proof: if ξ=0 then Vξ=4π21∫∣ψ′(x)∣2dx. This is usually much less work than transforming, and it is the standard route when ψ is given in closed form but its transform is not. If ξ=0, apply it to e−2πiξxψ(x) instead.
Form the ratio R=16π2VxVξ. Then R≥1, with R=1 precisely for Gaussians. The size of R is a dimensionless measure of how far the profile is from optimal; it is unchanged by translating, modulating or dilating ψ, so it is a property of the shape alone.
Check your arithmetic against the invariance in step 5: recompute R for ψ(x/δ) with some δ=1 and confirm that it does not move.
Example(A profile with R=9, and the Gaussian with R=1)Original
Take the odd profile
ψ(x)=cxe−πx2,c>0,
which is a Schwartz function and is not a Gaussian.
Normalisation.∫x2e−2πx2dx=4π1π/2π=4π21, using ∫x2e−αx2dx=21πα−3/2 with α=2π. Hence c2=4π2, and we may take c=(4π2)1/2=4.2156…
Centres.∣ψ∣2 is even, so x=0. Since ψ is real and odd, ψ^ is odd times i, so ∣ψ^∣2 is even and ξ=0. Step 2 of the method is therefore trivial.
Position variance. With α=2π and ∫x4e−αx2dx=43πα−5/2,
using ∫e−2πx2dx=1/2 together with the two moments already quoted. The first two terms cancel exactly, and the third collapses to 3π5/2/(2π)5/2=3/25/2=3/(42)=0.5303301. With c2=4π2 this gives ∫∣ψ′∣2=4π2⋅3/(42)=3π=9.424778. By step 4 of the method,
Vξ=4π23π=4π3=0.2387324…
The ratio.
VxVξ=(4π3)2=16π29=0.05699317,R=16π2VxVξ=9.
So this profile misses the bound by a factor of exactly nine: it is 3 times too wide in position and frequency together. (The pattern continues: the nth Hermite function has R=(2n+1)2, and the case n=0 below is the Gaussian.)
The Gaussian, for comparison. Take B=π in the equality case, so ∣A∣2=2π/π=2 and ψ0(x)=21/4e−πx2. From the closing paragraph of the proof, Vx=1/4B=1/4π and Vξ=B/4π2=1/4π, so
VxVξ=16π21=0.006332574,R=1,
and indeed ψ^0(ξ)=21/4e−πξ2=ψ0(ξ): this Gaussian is its own transform. The two computations differ by the factor 9, which is the whole content of the theorem in one number.
ApplicationOriginal
Position and momentum. In one-dimensional quantum mechanics a state is a normalised ψ, ∣ψ(x)∣2 is the probability density of position and, with p=hξ, ∣ψ^(ξ)∣2 is the density of momentum; this is the reading Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §4, pp. 159–161) set out. The theorem then says that the product of the two standard deviations is bounded below, and in the physical normalisation of Kreyszig (Kreyszig (1978), §11.2, Thm 11.2-2, p. 580) and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §19.2, p. 666) the bound is ΔxΔp≥ℏ/2. The equality case identifies the minimum-uncertainty states as the Gaussian wave packets, which is why coherent states of the harmonic oscillator are the ones used to model a classical trajectory as closely as quantum mechanics allows. Kreyszig stresses the correct reading: the limit is not an imperfection of instruments but a property of the state.
Pulse length against bandwidth. The same inequality with (t,ω) in place of (x,ξ) is the duration-bandwidth theorem of signal processing. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.2, (20.35), p. 973) make the point with a wave train of N cycles at frequency ω0: its spectrum has width Δω≈ω0/N, so clipping the pulse short necessarily broadens its frequency content. A radar or ultrasound designer who wants both a sharp range gate (short pulse) and a sharp Doppler estimate (narrow band) is asking for something the theorem forbids, and chirped pulses with matched filtering are the standard way to buy back resolution without violating it.
Spectral line widths. A transition with finite lifetime τ emits a wave train of finite duration, so its spectral line cannot be narrower than about 1/τ in frequency; this natural linewidth, and the analogous broadening of nuclear magnetic resonance lines by short relaxation times, is the uncertainty principle used as a measuring tool rather than as a limitation, since the observed width yields the lifetime.
Diffraction and beam divergence. For light passing an aperture, the transverse coordinate y and the transverse wave number ky=ksinθ are Fourier conjugates, so an aperture of width Δy produces an angular spread Δθ≳1/(kΔy)=λ/(2πΔy). Narrowing a slit to localise a beam therefore spreads it, and the Gaussian equality case is exactly the beam profile that diverges least, which is why laser cavities are designed around Gaussian modes.
Wrapping a rapidly decreasing function around the circle and sampling its Fourier transform at the integers produce the same periodic function, so a slowly convergent lattice sum can be traded for a rapidly convergent one.
the function obtained by wrapping f around a circle of circumference 1 and adding the overlapping layers. Whenever the series converges, Pf has period 1, because replacing x by x+1 only relabels the terms. Folland (Folland (1999), §8.5, p. 254) writes Pf=∑kTkf and treats it as an operator on L1.
Against this stands a second periodic function built from f, the sampled series
(Sf)(x)=n=−∞∑∞f^(n)e2πinx,
formed by reading off the Fourier transform at the integers and using those numbers as Fourier coefficients.
Convention, local to this unit. Throughout this unit f^ denotes the transform in the 2π-in-the-exponent convention,
f^(ξ)=∫Rf(x)e−2πixξdx,
which is the one Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §1, p. 134) and Folland (Folland (1999), §8.3, p. 249) use. It is not the handbook's default: the Fourier transform defines f^(k)=∫f(x)e−ikxdx, and the dictionary in its table of conventions converts between the two by ξ↦2πξ, that is, the f^(ξ) of this unit is the f^(2πξ) of that one. Every transform quoted below — the Gaussian and dilation rules of the first example, the Lorentzian pair of the second — is the corresponding entry of that unit's table read through this dictionary, and the check is shown where it is used. The reason for departing from the default here is that in this convention the formula carries no constants at all; the price, and the other conventions in which the sources state it, are in the note below.
Two decay classes are used. f belongs to the Schwartz space S(R) if it is smooth and xkf(m)(x) is bounded for all k,m≥0. The second is Stein and Shakarchi's (Stein & Shakarchi (2003), Ch. 5, §1.1, p. 131)moderate decrease, defined by them as: f is continuous and ∣f(x)∣≤A/(1+x2) for some constant A>0. Continuity is part of the definition and not an afterthought; a bounded function with a jump is not of moderate decrease however fast it decays. On f alone the condition is weaker than membership of S, but the form of the theorem that uses it — the one Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §3, p. 155) record as sufficient — asks for moderate decrease of f^ as well, which a Schwartz function supplies for nothing. It is a trade, not a straight loosening.
Theorem(Poisson summation formula)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, Thm 3.1 and the remark on p. 155, pp. 154–155
Let f∈S(R). Then for every x∈R
n=−∞∑∞f(x+n)=n=−∞∑∞f^(n)e2πinx,
both series converging absolutely. Setting x=0,
n=−∞∑∞f(n)=n=−∞∑∞f^(n).
Read the other way round: sampling f^ on the integer lattice returns precisely the Fourier coefficients of Pf, so Pf and Sf of the definition above are one and the same function. The conclusion survives, with the same proof, if f is only required to be of moderate decrease and f^ as well.
Theorem(The formula in Rn)AdaptedAdapted from Folland (1999), §8.5, Thm 8.32, p. 254
is Theorem 8.32 of Folland (Folland (1999), §8.5, Thm 8.32, p. 254), who assumes a continuous f on Rn obeying the two power bounds
∣f(x)∣≤C(1+∣x∣)−n−ε,∣f^(ξ)∣≤C(1+∣ξ∣)−n−ε,
with C and ε positive constants, and who records that under them each side converges absolutely, and uniformly on the torus Tn. This block quotes that theorem; the proof given below is written for n=1, and Folland's runs the same way on Tn, with the Fourier coefficients of the periodization computed over a fundamental cube instead of over [0,1]. The exponent is sharp in form rather than decorative: −n−ε is exactly what makes ∑k∈Zn(1+∣k∣)−n−ε finite, as a comparison with the corresponding integral over Rn shows.
Rescaling x turns the unit lattice into a lattice of spacing L>0. Applying the theorem to g(x)=f(Lx), whose transform is g^(ξ)=L−1f^(ξ/L), gives
n∈Z∑f(x+nL)=L1m∈Z∑f^(Lm)e2πimx/L.
The lattice of spacing L in x is paired with the reciprocal lattice of spacing 1/L in ξ: squeezing one spreads the other.
Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Exercise 13.9, p. 461) use f~(ω)=(2π)−1/2∫f(t)e−iωtdt and period X, and write the result as ∑nf(t+nX)=(2π/X)∑nf~(2nπ/X)e2πnit/X. Since f^(m/X)=2πf~(2πm/X), this is the displayed formula.
Strauss (Strauss (2008), Exercise 12.3.7, p. 349) uses F(k)=∫f(x)e−ikxdx and period 2π, and obtains ∑nf(2πn)=2π1∑nF(n), again the same statement because f^(m/2π)=F(m).
Stakgold and Holst (Stakgold & Holst (2011), §2.3, (2.3.26), p. 139) derive it distributionally from the completeness relation ∑kδ(x−2kπ)=2π1∑neinx applied to a test function, which is the same identity read as the Fourier series of a train of impulses; Riley, Hobson and Bence (Riley, Hobson & Bence (2006), Exercise 13.9, p. 461) start from the same impulse train. That route needs distributions to be legitimate, and is the quickest way to remember the formula; the proof below needs nothing beyond Fourier series.
Nothing in this unit mixes conventions: all displayed formulas use f^(ξ)=∫f(x)e−2πixξdx.
Idea of the proofOriginal
There are only two natural ways to manufacture a 1-periodic function out of a function on the line. One is geometric: cut the line into unit intervals and stack them, which is the periodization. The other is spectral: the inversion formula writes f as a continuous superposition ∫f^(ξ)e2πixξdξ of exponentials, and keeping only the integer frequencies — the ones that are already 1-periodic — leaves the sampled series. The theorem says these two constructions agree. The proof is correspondingly short: both candidates are continuous and 1-periodic, so it is enough to compare their Fourier coefficients, and the mth coefficient of the periodization is computed by unfolding a sum of integrals over consecutive unit intervals back into a single integral over R, which is f^(m) by definition.
ProofAdaptedAdapted from Stein & Shakarchi (2003), Ch. 5, Thm 3.1, p. 154, using Ch. 2, Thm 2.1, p. 39
Step 1: both sides are continuous and 1-periodic. Since f∈S(R), there is A with ∣f(y)∣≤A/(1+y2). For x in a compact set K the terms f(x+n) are therefore dominated by a summable sequence independent of x, so ∑nf(x+n) converges absolutely and uniformly on K and Pf is continuous. It is 1-periodic because n↦n+1 permutes the terms. The same argument applies to Sf: the Fourier transform of a Schwartz function is Schwartz, so ∑n∣f^(n)∣<∞ and Sf is a uniformly convergent series of continuous 1-periodic functions.
Step 2: the Fourier coefficients of the periodization. Fix m∈Z and compute Pf(m)=∫01(Pf)(x)e−2πimxdx. The series converges uniformly on [0,1] by step 1, so it may be integrated term by term:
Pf(m)=n=−∞∑∞∫01f(x+n)e−2πimxdx.
In the nth term substitute y=x+n. The interval [0,1] becomes [n,n+1], and e−2πimx=e−2πimye2πimn=e−2πimy because mn is an integer. This is the one place where it matters that the sampling points are the integers, that is, the lattice dual to the period.
Step 3: unfold the sum of integrals. The intervals [n,n+1], n∈Z, tile the line and overlap only at endpoints, so
the regrouping being justified by the absolute convergence of ∫∣f∣.
Step 4: the Fourier coefficients of the sampled series.Sf is a uniformly convergent trigonometric series, so again it may be integrated term by term, and ∫01e2πinxe−2πimxdx=δnm by orthonormality. Hence Sf(m)=f^(m) as well.
Step 5: equal coefficients force equal functions. The difference Pf−Sf is continuous, 1-periodic and has all Fourier coefficients zero, so by the uniqueness theorem for Fourier coefficients — an integrable function on the circle whose Fourier coefficients all vanish is zero at every point of continuity, which is Stein and Shakarchi's (Stein & Shakarchi (2003), Ch. 2, Thm 2.1, p. 39) Theorem 2.1 and is obtained in this handbook as a corollary of Fejér's theorem — it vanishes identically. That is the first formula; evaluating it at x=0 gives the second.
Step 6: moderate decrease. If instead f and f^ are only assumed to be of moderate decrease (continuity included, by the definition above), every step above goes through unchanged: step 1 used only the bound A/(1+y2) for f and the summability of ∣f^(n)∣, and steps 2–5 used only uniform convergence and integrability. This is the extension Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §3, p. 155) record. □
Remark(Where the hypotheses are used)Original
Decay of f is what makes the left side exist and converge uniformly, which is needed for the term-by-term integration in step 2. Decay of f^ does the same for the right side in step 1. Folland's hypotheses make this explicit: (1+∣x∣)−n−ε on both sides of the transform is exactly summability over the lattice in both variables.
The hypotheses cannot be dropped. Take f=1[0,1], the indicator of the unit interval. It obeys the size bound ∣f(x)∣≤A/(1+x2) but is discontinuous, so by the definition above it is not of moderate decrease; and f^(ξ)=(1−e−2πiξ)/(2πiξ) decays only like 1/∣ξ∣, so the hypothesis on f^ fails as well. Both sides can still be written down: f^(0)=1 and f^(n)=0 for every non-zero integer n, so the right side is 1, while the left side is f(0)+f(1)=2. The example does not separate the two hypotheses — for a compactly supported f a jump is precisely what forces f^ to decay no faster than 1/∣ξ∣ — but it does show where the proof stops: the periodization equals 1 off Z and 2 on Z, so it is discontinuous, and step 5 concludes nothing about a function that is not continuous. Assigning f the average of its one-sided limits at the two jumps replaces f(0)+f(1)=2 by 21+21=1 and restores the identity; this is the standard caution for a function with a jump. A corner is harmless: for the triangle f(x)=max(0,1−∣x∣) one has f^(ξ)=(sinπξ/πξ)2=O(ξ−2), both hypotheses hold, and both sides equal 1.
Periodization is meaningful more often than the formula is true.Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §3, p. 155) point out that the symmetric periodization of f(x)=1/x is πcotπx and that of 1/x2 is π2/sin2πx, although neither function is of moderate decrease and neither has a Fourier transform in the elementary sense. The lesson is that the left-hand side is the robust object; the equality with the right-hand side is what needs the hypotheses.
The lattice must be dual to the period. Step 2 works because e2πimn=1 for integers m,n. If one samples f^ at λZ with λ=1 while periodizing over Z, the phases no longer cancel and the computation of step 2 breaks down. It is not merely the argument that fails but the conclusion: for the Gaussian f(x)=e−πx2, which is its own transform, the left side at x=0 is ∑ne−πn2=1.0864348112, while the mismatched right side ∑nf^(λn)=∑ne−πλ2n2 equals 1.0098908590 at λ=1.3. (Both numbers are values of the theta function ϑ(s)=∑ne−πn2s studied in the first example below, at s=1 and s=λ2 respectively.) The correct pairing of spacings L and 1/L is in the note above.
Smoothness is not needed directly, only through decay. The formula is often quoted for Schwartz functions, but the proof uses no derivative of f. The role of smoothness is indirect: it is what guarantees decay of f^, by the smoothness-decay correspondence.
Method(Accelerating a lattice sum)Original
A sum ∑n∈Zf(n) whose terms die slowly can be replaced by one whose terms die fast, provided the transform is known.
Check the hypotheses: f continuous with ∣f(x)∣≲(1+∣x∣)−1−ε and the same for f^. A corner costs nothing, since it still leaves f^=O(ξ−2); a jump does, and at a jump f must be given the average of its one-sided limits, or that point treated separately. See the remark.
Compute f^ once, in closed form.
Decide which side converges faster. If f is concentrated on a scale a, then f^ is concentrated on the scale 1/a, so the f-side is slow exactly when a≫1 and fast when a≪1. Sum on the fast side.
If the parameter range makes neither side fast, split f=fshort+flong with fshort concentrated inside one lattice cell, sum fshort directly and apply Poisson summation to flong. This is the idea behind Ewald summation for Coulomb lattice sums.
Validate: evaluate both sides numerically at one convenient parameter value. They must agree to the accuracy of the slower side, and the discrepancy should shrink as more terms are kept.
For a lattice of spacing L, remember the prefactor 1/L and the sampling points m/L from the note above.
Example(The theta functional equation, and a numerical check)Original
For s>0 define the theta function
ϑ(s)=n=−∞∑∞e−πn2s.
Setting up. Apply the formula to f(x)=e−πsx2, a Schwartz function for each fixed s>0. In this unit's convention the Gaussian is its own transform, e−πx2(ξ)=e−πξ2: the table of pairs in the Fourier transform gives e−ax2(k)=π/ae−k2/4a in the handbook's default convention, and putting a=π and reading it at k=2πξ, as the dictionary of the definition block requires, turns it into e−πξ2. Dilating with f(x)=g(sx) where g(y)=e−πy2 — the dilation rule g(λ⋅)(ξ)=λ−1g^(ξ/λ) for λ>0 holds in both conventions, since the dictionary only relabels the variable — gives
f^(ξ)=s1g^(sξ)=s−1/2e−πξ2/s.
Applying the formula. The left side at x=0 is ∑ne−πn2s=ϑ(s) and the right side is ∑ns−1/2e−πn2/s=s−1/2ϑ(1/s). Hence
ϑ(s)=s−1/2ϑ(1/s)(s>0),
which is Stein and Shakarchi's (Stein & Shakarchi (2003), Ch. 5, Thm 3.2, p. 155) Theorem 3.2. Large s on the left corresponds to small s on the right: the identity exchanges a fast series for a slow one.
A number to check. Take s=4, so that 1/s=1/4 and the identity reads ϑ(1/4)=2ϑ(4).
with terms 0.4559381278, 0.0432139183, 0.0008514383, 0.0000034873, 0.0000000030, whose sum is 0.5000069747, giving
ϑ(41)=1+2(0.5000069747)=2.0000139494=2ϑ(4),
to all ten digits shown. The two series needed five terms and two terms respectively for the same accuracy, and the gap widens as s moves away from 1: near s=0 the left-hand series needs about s−1/2 terms and the right-hand one needs a handful.
Example(A slowly convergent sum done in closed form)Original
Let a>0 and f(x)=1/(x2+a2). Then f is continuous and of moderate decrease, and its transform is
f^(ξ)=aπe−2πa∣ξ∣,
again this unit's convention: the table in the Fourier transform gives (π/a)e−a∣k∣ for this f in the handbook's default convention, and the dictionary k=2πξ produces the exponent above. It is continuous and decays exponentially, so the extension of the theorem to moderate decrease applies to the pair (f,f^). At x=0,
The original series converges like ∑1/n2: its symmetric partial sums are 6.4890893 at N=5, 6.8111529 at N=50 and 6.8467583 at N=500, so after a thousand terms the relative error is still 5.8×10−4 and only the leading 6.8 is correct. The transformed series converges geometrically with ratio e−π=0.0432: keeping ∣m∣≤5 already gives 6.8507542143, correct to eight digits. This is step 3 of the method in action, and it is the whole practical point of the formula.
ApplicationOriginal
The heat kernel on a ring.Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §3.2, Thm 3.3, p. 156) deduce at once that periodizing the line's heat kernel produces the circle's: Ht(x)=∑nHt(x+n) with Ht(x)=(4πt)−1/2e−x2/4t. Physically this is the method of images — heat diffusing on a ring is heat diffusing on a line from infinitely many image sources — and the two representations are useful at opposite ends of the time range: the image sum converges in a term or two for short times, the eigenfunction sum ∑ne−4π2n2te2πinx for long times. See the heat kernel.
Diffraction from a crystal. A perfect lattice of scatterers has, as its Fourier transform, a lattice of delta functions on the reciprocal lattice; the general-period formula in the note is this statement for one dimension, and the Rn version of Folland is the three-dimensional statement used in crystallography. Bragg peaks appear only at reciprocal-lattice vectors for exactly this reason, and the inverse relation between the two spacings is why a finely spaced crystal gives widely spaced diffraction spots.
Ewald summation. The electrostatic energy of a periodic array of charges is a lattice sum of 1/r terms that converges only conditionally. Splitting each charge into a screened short-range part and a smooth long-range part, summing the first in real space and the second in reciprocal space by Poisson summation, turns an unusable sum into two exponentially convergent ones. This is the standard method in molecular dynamics and solid-state calculations, and step 4 of the method above is its skeleton.
Sampling. Poisson summation is the identity behind the reconstruction of a band-limited signal from its samples: periodizing the spectrum is exactly what sampling in time does, and aliasing is the overlap of the periodized copies. The statement and proof belong to the sampling theorem.
Number theory. The functional equation of the theta function proved above is the input to the functional equation of the Riemann zeta function, through the integral representation π−s/2Γ(s/2)ζ(s)=21∫0∞ts/2−1(ϑ(t)−1)dt recorded by Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 5, §3.1, p. 155); see the zeta function. It is also the mechanism behind the Euler-Maclaurin-style conversion of lattice sums into integrals plus exponentially small corrections, for which see the Euler-Maclaurin formula.
How to compute Fourier transforms on Rn: factorise when the function is a product, reduce to a single radial integral when it is spherically symmetric, and read off the transforms of the Yukawa and Coulomb potentials.
For f:Rn→C absolutely integrable, this unit uses throughout the symmetric convention of Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.10, (13.49)–(13.50), p. 452):
with k⋅r=k1r1+⋯+knrn and dnr Lebesgue measure. Writing k⋅r as a sum shows that the kernel factorises, e−ik⋅r=∏j=1ne−ikjrj: the n-dimensional transform is n one-dimensional transforms, one per coordinate. Composing the two formulas gives the delta representation δ(n)(r)=(2π)−n∫eik⋅rdnk of Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.10, (13.51), p. 452); see delta identities.
A function f on Rn is radial (spherically symmetric) if it depends only on ∣r∣, that is if there is f0 on [0,∞) with f(r)=f0(∣r∣). Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 6, §2, p. 182) note the equivalent formulation: f is radial exactly when f(Rr)=f(r) for every rotation R, because any two points of the same length are carried onto one another by some rotation. We write r=∣r∣ and k=∣k∣ and abuse notation by writing f(r) and f~(k).
Converting is mechanical: Strauss's F is (2π)n/2f~, Stein and Shakarchi's f^(ξ) is (2π)n/2f~(2πξ), and Arfken's g(k) is f~(−k). Arfken's opposite sign is invisible on radial functions, because their transforms are even. Two ways to see it: item 4 of the rules below applies to R=−I, which is orthogonal in every dimension (though a rotation only when n is even, since det(−I)=(−1)n), and gives f~(−k)=f(−⋅)(k)=f~(k) for even f; or directly, sin(k⋅r) integrates to zero against a radial density. Either way, every result quoted from Arfken below is radial, so the two conventions agree there. When a statement of Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 6, Prop. 2.1 and Thm 2.4, pp. 181–183) is used below it is restated in the convention of this unit, never mixed with it.
Let f lie in the Schwartz space S(Rn) of smooth functions all of whose derivatives decay faster than every power of 1/∣r∣. Then, in the convention of this unit:
Translation and modulation.f(r+h)⟼eik⋅hf~(k) and eih⋅rf(r)⟼f~(k−h).
Dilation. For δ>0, f(δr)⟼δ−nf~(k/δ).
Derivatives.∂f/∂rj⟼ikjf~(k), hence ∇f⟼ikf~ and ∇2f⟼−∣k∣2f~. Conversely −irjf(r)⟼∂f~/∂kj.
Rotations and reflections.f(Rr)⟼f~(Rk) for every orthogonal R. Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 6, Prop. 2.1 (vi), p. 181) state this for rotations; the proof below uses only ∣detR∣=1 and R−1=RT, so it covers the orthogonal matrices of determinant −1 as well, R=−I in odd dimensions among them.
Separability. If f(r)=∏j=1nfj(rj) then f~(k)=∏j=1nf~j(kj), each factor a one-dimensional transform.
Inversion and Plancherel. The inverse formula above recovers f, and ∫Rn∣f~∣2dnk=∫Rn∣f∣2dnr.
Convolution.f∗g(k)=(2π)n/2f~(k)g~(k), where (f∗g)(r)=∫f(y)g(r−y)dny.
Items 1–4 are Stein and Shakarchi's (Stein & Shakarchi (2003), Ch. 6, Prop. 2.1 (i)–(vi), p. 181) Proposition 2.1 (item 4 for rotations, as noted there); item 5 is the factorisation carried out in Step 1 of the proof of their Theorem 2.4 (Stein & Shakarchi (2003), Ch. 6, Thm 2.4, Step 1, pp. 182–183), where it is applied to the Gaussian, and in general it is immediate from the product form of the kernel recorded in the definition above; items 6–7 are that theorem and Step 4 of its proof. All are restated here in our convention; the prefactor in item 7 is the price of the symmetric normalisation and is absent in theirs.
and f~(0)=(2π)−3/2∫0∞4πr2f(r)dr, the value obtained by letting k→0.
Reduction in two dimensions. If f(ρ)=f(ρ) is radial on R2 and ∫0∞ρ∣f(ρ)∣dρ<∞ — the same absolute-convergence condition as above, with rn−1 in place of r2, needed to interchange the radial and angular integrations — then, by Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.1.10, (13.52), p. 453), who state the formula without a convergence hypothesis,
Two structural facts do all the work, and they pull in opposite directions. Factorisation says that if the function splits as a product over coordinates, so does its transform, which reduces an n-fold integral to n one-fold ones; this is how the Gaussian is handled. Rotation invariance says that the transform commutes with rotations, so a function with no preferred direction has a transform with no preferred direction either; this lets us align the axes with k after the fact, which collapses the angular part of the integral and again leaves a single radial integral. In three dimensions the angular integral is elementary — the exponential e−ikrcosθ integrates against sinθdθ to give 2sin(kr)/(kr) — and everything reduces to a sine transform of rf(r). The Coulomb potential then needs one more idea: its integral does not converge, so it is reached as the limit of the Yukawa potential, whose exponential factor supplies the missing decay.
Step 1: the transform commutes with orthogonal maps. Let R be orthogonal and g(r)=f(Rr). Substituting y=Rr, whose Jacobian has ∣detR∣=1 (the determinant itself is ±1, and only its modulus enters),
Because R is orthogonal, R−1=RT and hence k⋅R−1y=Rk⋅y, so the last integral is f~(Rk). This is item 4 of the rules, in the orthogonal generality stated there; rotations are the case detR=+1, which is the one Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 6, Prop. 2.1 (vi), p. 181) state.
Step 2: radial in, radial out. If f is radial then f(Rr)=f(r) for every rotation R, so step 1 gives f~(Rk)=f~(k) for every R; by the same equivalence, f~ is radial. The only property of f used is its rotation invariance, which is why this holds in every dimension.
Step 3: choose the axis after the fact. Now let n=3 and let f be radial. By step 2 we may compute f~ at any one k of length k and the answer will serve for all of them, so choose spherical polar coordinates with the polar axis along k. Then k⋅r=krcosθ and d3r=r2sinθdrdθdφ, giving
The hypothesis ∫0∞r2∣f∣dr<∞ makes the triple integral absolutely convergent, so by Fubini's theorem it may be evaluated in any order. This is the step that fails for the Coulomb potential.
Step 4: do the angles. Nothing depends on φ, so that integral contributes 2π. For the polar angle substitute u=cosθ, du=−sinθdθ:
which is the claim. As k→0 the factor sin(kr)/(kr)→1 pointwise and is bounded by 1 everywhere, so ∣4πr2f(r)sin(kr)/(kr)∣≤4πr2∣f(r)∣, which is integrable by hypothesis; dominated convergence then gives f~(0)=(2π)−3/2∫4πr2f(r)dr, the total mass of f up to the prefactor. □
Remark(Where the hypotheses are used)Original
Absolute convergence is used in step 3, to interchange the three integrations. For f(r)=1/r the radial integral ∫0∞r2⋅r−1dr=∫0∞rdr diverges, so the hypothesis fails and the reduction is not available. Both Hassani (Hassani (2013), §9.2, Example 9.2.3, p. 282) and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.2, Example 20.2.5, p. 974) say so explicitly and obtain the Coulomb transform only as a limit. Note that the trouble is at infinity, not at the origin: r2⋅r−1 is perfectly integrable near 0, which is why the singular-looking 1/r is in fact the easy end.
Rotation invariance of f is used in step 3, and nowhere else. Without it, aligning the axis with k would change the integrand, and the angular integral would not collapse. For a function with an axis of symmetry but not full spherical symmetry the same trick gives a double, not a single, integral.
Orthogonality of R, not just invertibility, is used in step 1, twice: once for ∣detR∣=1, so the volume element is unchanged, and once for k⋅R−1y=Rk⋅y. A general linear change of variables A gives instead f(Ar)⟼∣detA∣−1f~(A−Tk), with the inverse transpose rather than A itself; the two coincide exactly for orthogonal A.
The limit α→0 in the Yukawa formula is not automatic. The pointwise limit of 2/π(k2+α2)−1 is 2/πk−2, and the sources take this as the transform of 1/r; making it a theorem rather than a convention requires reading 1/r as a tempered distribution and the convergence as convergence of distributions, for which see tempered distributions. The limit is not uniform near k=0, which is the analytic trace of the fact that the Coulomb potential has infinite range.
Separability needs a genuine product, not merely a radial function. The Gaussian is the one function that is both radial and a product, which is exactly why it can be transformed either way and why it is the test case for every convention.
In two dimensions the angular integral is not elementary, because ∫02πe−ikρcosφdφ=2πJ0(kρ) rather than a ratio of elementary functions. Three dimensions is the lucky case; the general pattern is the Hankel transform of order n/2−1.
Method(Computing a transform on Rn)Original
Fix the convention first and write it down. Nearly every discrepancy between two published transforms of the same function is a factor of (2π)n/2 or a sign; see the table above. Carry the convention through to the final answer.
Look for a product. If f(r)=∏jfj(rj), transform each factor separately in one dimension and multiply. This covers Gaussians, products of exponentials, and anything built by separation of variables.
Otherwise look for spherical symmetry. If f is radial, check ∫0∞rn−1∣f∣dr<∞, then use the radial reduction: in R3 compute the single integral k4π∫0∞rf(r)sin(kr)dr and divide by (2π)3/2; in R2 use the J0 formula.
If the radial integral diverges, insert a convergence factor. Replace f(r) by f(r)e−αr with α>0, transform, then let α→0 in the answer. This is how 1/r is handled, and the parameter α usually has a physical reading as an inverse screening length.
Differentiate in a parameter to get new pairs for free. Differentiation with respect to a parameter commutes with the transform, since they integrate over different variables; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.2, (20.44), p. 974) use exactly this to get the transform of e−Zr from that of e−αr/r.
Check three things. (i) f~(0) must equal (2π)−n/2∫fdnr. (ii) If f is real and even, f~ is real. (iii) Plancherel: ∫∣f~∣2dnk must equal ∫∣f∣2dnr when both are finite. A numerical quadrature of the radial integral at one value of k is the fastest check of all.
Example(Yukawa and Coulomb potentials, and the hydrogenic 1s orbital)Original
The Yukawa potential. Take f(r)=e−αr/r with α>0 on R3. It is radial, and ∫0∞r2∣f∣dr=∫0∞re−αrdr=α−2<∞, so step 3 of the method applies. The radial integral is elementary because the factor r cancels:
∫0∞r⋅re−αrsin(kr)dr=∫0∞e−αrsin(kr)dr=α2+k2k,
the standard integral obtained by writing sinkr=Imeikr and integrating e−(α−ik)r. Hence
using 4π/(2π)3/2=2/2π=2/π. This is Hassani's (Hassani (2013), §9.2, Example 9.2.3, p. 283) result V~α=(2q/2π)(k2+α2)−1 for charge q, and Arfken, Weber and Harris' (Arfken, Weber & Harris (2013), §20.2, (20.41), p. 974) equation (20.41), who note that without the prefactor the answer is the familiar 4π/(k2+α2).
A number. At α=1 and k=2,
π2⋅4+11=50.7978845608=0.1595769122,
equivalently 4π/5=2.5132741229 in the unnormalised convention. A direct numerical evaluation of k4π∫0∞e−rsin2rdr returns 2.5132741229, so the two agree to ten digits.
The Coulomb limit. Letting α→0,
(r1)(k)=π2k21,
which is Arfken, Weber and Harris' (Arfken, Weber & Harris (2013), §20.2, (20.42), p. 974) equation (20.42), (2π)−3/24π/k2, and Hassani's (Hassani (2013), §9.2, p. 283)V~Coul=(2q/2π)k−2. At k=2 its value is 0.7978845608/4=0.1994711402. The warning of the remark applies: the direct integral does not converge, and this equality is a statement about distributions.
One more pair for free. Differentiating the Yukawa result with respect to α under the integral sign, and noting −∂α(e−αr/r)=e−αr,
(e−αr)(k)=−∂α∂[π2k2+α21]=π2(k2+α2)22α,
the transform of the unnormalised hydrogenic 1s orbital, matching Arfken, Weber and Harris' (Arfken, Weber & Harris (2013), §20.2, (20.44), p. 974)8πZ/[(2π)3/2(k2+Z2)2] since 8π/(2π)3/2=22/π. At α=1, k=2 this is 0.7978845608×2/25=0.0638307649, again confirmed to ten digits by numerical quadrature of k4π∫0∞re−rsin2rdr=1.0053096491 divided by (2π)3/2=15.7496099457.
Check by separability. For the Gaussian f(r)=e−∣r∣2 both routes are open. As a product, f=∏j=13e−rj2 and each one-dimensional transform is (2π)−1/2∫e−r2e−ikrdr=2−1/2e−k2/4, so f~(k)=2−3/2e−∣k∣2/4. As a radial function, the reduction gives (2π)−3/2(4π/k)∫0∞re−r2sinkrdr, and at k=1.5 the integral evaluates to 0.3787177, so the formula returns (15.7496099)−1×(4π/1.5)×0.3787177=0.2014486; the product route gives 2−3/2e−0.5625=0.3535534×0.5697828=0.2014486. The two agree, as they must.
ApplicationOriginal
Screened interactions. The pair e−αr/r⟷const/(k2+α2) is the most used transform in many-body physics. Read forwards it says that a screened Coulomb interaction has a momentum-space form regular at k=0; read backwards it says that any interaction whose momentum dependence is (k2+α2)−1 has exponential range 1/α in real space. Hassani (Hassani (2013), §9.2, Example 9.2.3, p. 283) records that this is why Yukawa introduced α: the short range of the nuclear force demanded it, and for electromagnetism, where the range is unbounded, α=0.
Scattering and the Born approximation. In the Born approximation the scattering amplitude is proportional to the Fourier transform of the potential evaluated at the momentum transfer q. For the Coulomb potential this gives an amplitude ∝1/q2 and hence a cross-section ∝1/q4, which is Rutherford's law; the 1/k2 computed above is the whole of that derivation. See the Born approximation.
Form factors.Hassani (Hassani (2013), §9.2, Example 9.2.4, pp. 283–284) extends the calculation to a charge spread with density ρ: the transform of the potential factorises as V~(k)=4πqρ~(k)/∣k∣2, so a plot of k2V~(k) against k displays ρ~, the form factor, directly. A constant means a point charge; any departure means structure. Hassani notes that the relativistic version of this argument was the evidence for point-like constituents inside the proton in 1968.
Solving equations with constant coefficients. Because ∇2⟼−∣k∣2, Poisson's equation ∇2ϕ=−4πρ becomes ∣k∣2ϕ~=4πρ~ and the screened equation (∇2−α2)ϕ=−4πρ becomes (∣k∣2+α2)ϕ~=4πρ~. Inverting the first recovers the Coulomb Green function 1/(4π∣r∣) and the second the Yukawa one; this is the standard route to the fundamental solutions in Laplace's equation and the Helmholtz equation.
Crystallography and imaging. A diffraction experiment measures ∣f~(k)∣2 for the electron density f of a sample, and the rotation rule is what makes powder diffraction work: averaging over orientations of the sample is averaging f~ over rotations of k, so for a radial density nothing is lost. The same rule underlies the reconstruction of radially symmetric objects from a single one-dimensional profile in small-angle scattering.
The one-sided transform fˉ(s)=∫0∞f(t)e−stdt, the half-plane on which it converges, its standard pairs and its derivative, shift and scaling rules, together with the Bromwich contour integral that inverts it.
whenever the integral exists. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.2, (13.53), p. 453) write it this way and take s real; Evans (Evans (2010), §4.3.3, (43), p. 192) defines it for u∈L1(R+) and s≥0; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.7, p. 1008 and footnote 6) allow complex s. Only the values of f on (0,∞) enter, which is why the transform suits an initial-value problem, where data are given at t=0 and nothing is known for t<0. It is linear: L[af1+bf2]=afˉ1+bfˉ2 wherever both transforms exist.
f is of exponential order if there are constants s0, M and t0≥0 such that
e−s0tf(t)≤Mfor all t>t0,
the condition of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.7, (20.127), p. 1008). Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.2, p. 453) record the resulting picture: for a given f there is in practice a real number s0, the abscissa of convergence, with the integral existing for s>s0 and diverging for s≤s0. Each line of the table below carries its own s0.
Existence. If f is of exponential order with constant s0, and is integrable on (0,t0), then fˉ(s) exists for s>s0. Allowing complex s, Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.10, p. 1039) require Res≥β for convergence when f grows no faster than eβt, so the natural domain is a right half-plane.
Decay of the transform.Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.7, p. 1010) state that whenever fˉ is a Laplace transform, fˉ(s)→0 as s→∞, and conclude that a function of s growing like a positive power of s is the transform of nothing. Step 1 of the proof below establishes this for every f of exponential order, which is the case the statement is about. It does not cover the last row of the table: the same book's convention L[δ(t)]=1 for an impulse at the endpoint produces a transform that does not decay, and δ is not a function of exponential order. The two statements sit side by side in Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.7, pp. 1010–1011) without being reconciled; we report both and read the decay claim as a statement about transforms of functions.
Standard pairs. With s0 as given by Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.2, Table 13.1, p. 455):
f(t)
fˉ(s)
s0
1
1/s
0
tn(n>−1)
Γ(n+1)/sn+1
0
eat
1/(s−a)
a
tneat
n!/(s−a)n+1
a
sinbt
b/(s2+b2)
0
cosbt
s/(s2+b2)
0
sinhat
a/(s2−a2)
∣a∣
coshat
s/(s2−a2)
∣a∣
eatsinbt
b/[(s−a)2+b2]
a
t−1/2
(π/s)1/2
0
H(t−t0)
e−st0/s
0
δ(t−t0)
e−st0
0
The entry for tn is Arfken, Weber and Harris' (Arfken, Weber & Harris (2013), §20.7, (20.136), p. 1009) (20.136), valid for every real n>−1 through the gamma function; for integer n it reduces to n!/sn+1. The entries for H and δ use the Heaviside step, and L[δ(t)]=1 is the convention Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.7, (20.140), p. 1011) adopt for an impulse sitting at the endpoint. Strauss (Strauss (2008), §12.5, (10)–(12), p. 353) adds the pairs a(4πt3)−1/2e−a2/4t↔e−as and (πt)−1/2e−a2/4t↔s−1/2e−as, which are what make the transform useful for diffusion on a half-line.
Definition(Operational rules)AdaptedAdapted from Riley, Hobson & Bence (2006), §13.2.1–§13.2.2, (13.57)–(13.64), pp. 455–459; Strauss (2008), §12.5, properties (i)–(ix), p. 354
Let fˉ=L[f] and gˉ=L[g], let a be a real constant, b>0 and n∈N. Only rule 5 needs a>0; rule 3 and the table entries built on eat hold for either sign, and the half-plane moves with a, which is why the table's s0 column reads a and not ∣a∣ for those rows.
Derivatives.L[f′]=sfˉ(s)−f(0), and in general
L[f(n)]=snfˉ(s)−sn−1f(0)−sn−2f′(0)−⋯−f(n−1)(0).
Integral.L[∫0tf(u)du]=s1fˉ(s).
Shift in s.L[eatf(t)]=fˉ(s−a): multiplying by eat moves the origin of s.
Delay in t.e−bsfˉ(s)=L[g] where g(t)=0 for 0<t≤b and g(t)=f(t−b) for t>b; equivalently L[H(t−b)f(t−b)]=e−bsfˉ(s).
Scaling. For a>0, L[f(at)]=a−1fˉ(s/a).
Multiplication by tn.L[tnf(t)]=(−1)ndnfˉ/dsn.
Division by t.L[tf(t)]=∫s∞fˉ(u)du, provided limt→0f(t)/t exists.
Convolution.L[∫0tf(u)g(t−u)du]=fˉ(s)gˉ(s).
Rules 1–7 are Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §13.2.1–§13.2.2, pp. 455–457) (13.57)–(13.63) and rule 8 their convolution theorem. Strauss' (Strauss (2008), §12.5, properties (i)–(ix), p. 354) table of properties carries seven of the eight — all but the integral rule 2 — as (ii)–(ix), with rule 1 written out only for n=1,2 and rule 6 only for n=1; his remaining entry (i) is linearity, which we have put in the definition block above.
Let f be of exponential order, so that fˉ(s)=∫0∞e−stf(t)dt converges for Res>s0, and let f in addition be piecewise smooth on every bounded subinterval of (0,∞). Then for every t>0
2f(t+)+f(t−)=2πi1∫λ−i∞λ+i∞estfˉ(s)ds,
the Bromwich integral, taken as a symmetric limit up the vertical line Res=λ, where λ is chosen so that every singularity of fˉ lies strictly to its left. At every t where f is continuous the left-hand side is f(t), so the formula recovers f at every point of continuity. At a jump it returns the midpoint, which is in general not f(t); the remark below explains why no formula built from fˉ could do better.
The regularity hypothesis is not decoration. Both sources hedge at exactly this point: Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.10, (20.199)–(20.201), pp. 1038–1039) ask of the damped function exactly the restriction its Fourier-integral representation needs, which at the opening of the section they name as the Dirichlet conditions, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.5, p. 884) call a proof valid for arbitrary f and fˉ hard to come by, and offer a verification in its place. Piecewise smoothness is the hypothesis under which Fourier inversion is available in this handbook, and it is what step 4 of the proof below uses.
On the contour: Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.5, (25.25), p. 884) state the formula with λ>0 and all singularities of fˉ to the left of the line; Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.10, (20.203), p. 1039) write β for λ and note that fˉ is analytic in Res≥β, so the integrand can be singular only to the left. It is also called the Fourier-Mellin integral.
Evaluation by residues. If t>0 and fˉ is analytic apart from isolated singularities, has no branch points, and is small for large ∣s∣, the contour may be closed by a large semicircle in the left half-plane which contributes nothing, and then
2f(t+)+f(t−)=∑(residues of estfˉ(s) at the singularities with Res<λ),
which is Arfken, Weber and Harris' (Arfken, Weber & Harris (2013), §20.10, (20.204), p. 1039) equation (20.204), written by them with F(t) on the left under the same reading as above. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.5, p. 885) describe the completions to use when instead fˉ has infinitely many poles (a sequence of arcs threading between them, giving f as a series) or branch points (a cut plane with a keyhole contour).
Idea of the proofOriginal
The Laplace transform is the Fourier transform in disguise, and the disguise is a single exponential factor. Given f on (0,∞) that may grow like eβt, the product e−βtf(t), extended by zero to t<0, decays and so has a Fourier transform. Writing that transform out and renaming the variable s=β+iu turns it into fˉ(s): the Laplace transform is the Fourier transform evaluated along a vertical line in the complex plane rather than along the real axis. Everything then follows. Existence needs only that the damping beat the growth, which is the exponential-order condition. Inversion is the Fourier inversion formula with the same change of variable, and the real line of frequencies becomes the vertical Bromwich line, which must be far enough to the right that the damping still works. The derivative rule is integration by parts, and the boundary term at t=0 — absent for the Fourier transform, which has no boundary — is exactly what lets initial data enter an algebraic equation.
Proof(Existence, the derivative rule, and the Bromwich formula)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.7, p. 1008 and §20.10, (20.199)–(20.203), pp. 1038–1039; Riley, Hobson & Bence (2006), §13.2.1, (13.57), p. 455
Step 1: existence. Suppose ∣f(t)∣≤Mes0t for t>t0 and f is integrable on (0,t0). For real s>s0,
and the integral over (0,t0) is finite by assumption. The same bound holds for complex s with Res>s0, since ∣e−st∣=e−tRes, so the transform converges absolutely on a right half-plane. Differentiating under the integral sign, which the same bound justifies, shows fˉ is analytic there. The estimate also gives fˉ(s)→0 as Res→∞.
Step 2: the derivative rule. Integrate by parts, assuming f is differentiable with f′ again of exponential order:
the boundary term at infinity vanishing for Res large enough. Applying this to f′ in place of f gives L[f′′]=s2fˉ−sf(0)−f′(0), and induction gives rule 1 in general. The initial values appear because the integral starts at t=0: the transform of a derivative remembers where the function started.
Step 3: damp the function. Let f be of exponential order and piecewise smooth, and choose β>s0. Put
g(t)={e−βtf(t),0,t>0,t≤0.
Then g is absolutely integrable on R, since ∫0∞e−βt∣f∣dt converges by step 1. It is also piecewise smooth on every bounded interval: multiplying by the smooth factor e−βt and cutting off at t=0 creates at worst one further jump, at the origin.
Step 4: apply Fourier inversion to g. Integrability of g by itself does not license a pointwise inversion formula — the inversion integral need not even converge — so we use the piecewise-smooth version of Fourier inversion, whose hypotheses are exactly what step 3 supplies: g∈L1(R) satisfying Dirichlet conditions. It returns the midpoint value and requires the integral to be read as a symmetric limit. In the angular-frequency convention of the Fourier transform,
Step 5: undo the damping and change variable. Fix t>0. The factor eβt is continuous, so multiplying the last display by it turns the left-hand side into the midpoint value of f at t:
Now set s=β+iu, so that ds=idu and u running from −∞ to ∞ corresponds to s running up the vertical line Res=β. Then du=ds/i and
2f(t+)+f(t−)=2πi1∫β−i∞β+i∞estfˉ(s)ds,
the Bromwich integral, the contour integral again being a symmetric limit. This is the route of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.10, (20.200)–(20.203), pp. 1038–1039). For t<0 the same computation returns 0, because g vanishes there; the formula reconstructs the one-sided function, not a two-sided one.
Step 6: residues. Fix t>0 and suppose fˉ is analytic except at finitely many isolated singularities, all with real part less than β, and that ∣fˉ(s)∣→0 uniformly as ∣s∣→∞. Close the vertical segment from β−iR to β+iR with the semicircular arc ΓR of radius R to the left, parametrised by s=β+Reiθ, 2π≤θ≤23π. The crude bound ∣est∣≤eβt is not enough, because the arc has length πR; what saves it is that ∣est∣=eβtetRcosθ is exponentially small except near the two ends. Writing θ=2π+φ and using sinφ≥2φ/π on [0,2π],
a bound independent of R. With M(R)=maxΓR∣fˉ∣→0 the arc contributes at most M(R)πeβt/t→0. This is Jordan's lemma; see infinite real integrals and Jordan's lemma. The closed contour encircles every singularity once, counter-clockwise, so the residue theorem gives
2f(t+)+f(t−)=j∑Ress=sj(estfˉ(s)),
which is f(t) at every t>0 where f is continuous. When fˉ is rational the right-hand side is a finite sum of terms pj(t)esjt and hence smooth: it is then the continuous representative of f, which is what the method below reads off. □
Remark(Where the hypotheses are used)Original
Exponential order is what makes step 1 work, and it can fail.Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.7, p. 1008) give f(t)=et2: no s0 dominates it, and L[et2] exists for no s at all. The transform can also fail at the other end, from a singularity at the origin rather than growth at infinity: ∫0∞tne−stdt diverges at 0 for n≤−1, so L[t−1] does not exist, which is why the table entry for tn carries the restriction n>−1.
The abscissa is a property of f, not of the formula.L[eat]=1/(s−a) was derived for s>a, but the right-hand side is defined and analytic for all s=a. The transform is the analytic continuation of the integral, and the integral's half-plane of convergence is recorded only so that the Bromwich line can be placed correctly. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.10, p. 1040) address exactly this apparent paradox when closing the contour through Res<0: the restriction Res≥β defined fˉ, and once defined it may be used wherever it is analytic.
The contour must be to the right of everything. If λ is placed to the left of a pole of fˉ, step 3 fails — the damping e−λt no longer beats the growth of f — and the integral returns a different function. For fˉ(s)=1/(s−2) the correct choice λ>2 gives e2t, while closing to the right of the pole with λ<2 gives −e2tH(−t).
Closing to the left needs t>0 and decay of fˉ. Both are used in step 6: for t<0 the exponential est grows in the left half-plane and the arc does not vanish, which is the mechanism by which the formula returns 0 there. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.10, Example 20.10.2, p. 1041) show a case, fˉ(s)=(1−e−as)/s, where the two terms of the integrand must be closed on opposite sides depending on whether t<a or t>a, producing a square pulse.
No branch points is a real restriction. Transforms such as e−as from Strauss' (Strauss (2008), §12.5, (10), p. 353) diffusion table have a branch point at s=0 and are not inverted by a finite residue sum; a keyhole contour around the cut is needed, as Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §25.5, p. 885) describe.
Inversion is unique only up to negligible sets, which is why the Bromwich statement is a midpoint statement.Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.2, footnote, p. 454) note that two functions differing at finitely many isolated points have the same transform. No formula built from fˉ can therefore distinguish them, and the one the Bromwich integral picks is the midpoint: at a jump it returns 21[f(t+)+f(t−)], not f(t). A table lookup likewise recovers f only where f is continuous. This is also why the statement block assumes f piecewise smooth: that is what makes the one-sided limits exist and the symmetric limit converge, by the piecewise-smooth form of Fourier inversion used in step 4.
Finally, the derivative rule needs f to be differentiable with f′ of exponential order, and it needs f(0) to exist. A forcing term with a jump at t=0 must be written with a Heaviside factor and handled by rule 4, not by rule 1.
Method(Solving a linear constant-coefficient initial-value problem)AdaptedAdapted from Riley, Hobson & Bence (2006), §15.1.5, (15.31)–(15.32) and the surrounding scheme, p. 501
For any(n)+⋯+a1y′+a0y=g(t) with y(0),…,y(n−1)(0) given:
Transform both sides, using rule 1 on each derivative. The initial data enter as the polynomial ∑ksn−1−ky(k)(0) terms; no separate step is needed to impose them. This is the advantage Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §15.1.5, p. 501) stress: the solution arrives in one piece rather than as a complementary function plus a particular integral.
Solve the resulting algebraic equation for yˉ(s). It is always yˉ=[P(s)+gˉ(s)]/Q(s), with Q the characteristic polynomial and P carrying the initial data.
Put yˉ over a common denominator and record the poles. The poles of Q are the natural frequencies of the system; those coming from gˉ are the driven response.
Invert. If yˉ is rational, split it by partial fractions and read off each piece from the table; equivalently, and often faster when the poles are simple, evaluate the Bromwich integral as a sum of residues. If yˉ contains a factor e−bs, strip it off, invert the rest, and reinstate it with the delay rule 4.
Check three things: y(0) and the other initial values against the data; the residual any(n)+⋯−g at one value of t; and lims→∞syˉ(s), which must equal y(0).
Example(A driven damped system, inverted by residues)Original
Solve
y′′+3y′+2y=e−3t,y(0)=1,y′(0)=0.
Step 1: transform. By rule 1, L[y′′]=s2yˉ−s⋅1−0 and L[y′]=syˉ−1, while L[e−3t]=1/(s+3) for s>−3. Hence
obtained by putting the two terms over the common denominator. The poles are at s=−1,−2,−3, all simple, so the abscissa of convergence is s0=−1 and any λ>−1 serves as the Bromwich line.
Step 3: invert by residues. The numerator is quadratic and the denominator cubic, so yˉ(s)=O(1/∣s∣) at infinity and the contour may be closed to the left for t>0. Each pole is simple, so the residue of estyˉ(s) at sj is esjt times the numerator over the product of the other two factors:
Step 4: the checks of the method.y(0)=2.5−2+0.5=1, as required. y′(t)=−25e−t+4e−2t−23e−3t, so y′(0)=−2.5+4−1.5=0. Substituting the e−3t term alone into the left-hand side gives 21(9−9+2)e−3t=e−3t, and the e−t and e−2t terms are annihilated by (d/dt+1)(d/dt+2), so the equation holds. Finally syˉ(s)→1=y(0) as s→∞.
A second, shorter check of the table. The entry L[tn]=Γ(n+1)/sn+1 at n=21 gives L[t1/2]=Γ(3/2)/s3/2=21π/s3, which is Riley, Hobson and Bence's (Riley, Hobson & Bence (2006), §13.2, Table 13.1, p. 455) table entry. At s=4 this is 21π/64=0.1107784, and the direct integral ∫0∞t1/2e−4tdt=(π/2)/8=π/16=0.1107784 agrees.
ApplicationOriginal
Circuit and mechanical transients. A series RLC circuit driven by a source obeys a constant-coefficient equation in the charge, and transforming it turns the circuit into an algebraic divider: the transform of the response is the transform of the drive divided by Ls2+Rs+1/C. The roots of that quadratic are the natural frequencies, the residues at them are the amplitudes of the transient, and the poles contributed by the drive are the steady state. Switching an input on at time b costs only a factor e−bs, by rule 4, and that is why a system started from rest is handled with this one-sided transform instead of its two-sided relative.
Transfer functions and stability. In control theory the ratio yˉ(s)/gˉ(s)=1/Q(s) is the transfer function, and the residue calculation of the worked example shows why its poles decide everything: each pole sj contributes a term esjt, so the response decays exactly when every pole has Resj<0. The Bromwich line, which must lie to the right of all poles, is then in the right half-plane only for an unstable system.
Diffusion and waves on a half-line. Transforming in time rather than in space turns a partial differential equation into an ordinary one in x with s as a parameter. Strauss (Strauss (2008), §12.5, Examples 4–5, pp. 355–356) solves utt=c2uxx on x>0 with boundary data f(t) this way: the transformed equation gives U(x,s)=fˉ(s)e−sx/c, and the delay rule reads off u(x,t)=H(t−x/c)f(t−x/c), a signal travelling at speed c. The same method on the diffusion equation produces e−xs/k and, through the s pairs in the table, the error-function solution. This is the technique that handles inhomogeneous boundary conditions most cleanly.
Resolvents and semigroups.Evans (Evans (2010), §4.3.3, Example 8, pp. 192–193) observes that transforming the heat equation vt=Δv, v∣t=0=f, in time yields −Δu+su=f for u=vˉ(⋅,s): the Laplace transform of the heat flow is the resolvent of the Laplacian. This is the bridge between the transform and operator semigroup theory, where the Bromwich integral reappears as the contour formula recovering a semigroup from its resolvent.
Probability and renewal problems. For a non-negative random variable with density p, pˉ(s)=E[e−sT] is the moment generating function at −s, and rule 6 turns differentiation in s at s=0 into the moments. Because the convolution rule 8 matches the sum of independent waiting times, renewal equations become algebraic, which is the standard route through queueing and reliability calculations.
Shows that every rational function splits uniquely into a polynomial plus one singular part per pole, gives the cover-up rule and its derivative extension for the coefficients, the real form that keeps complex-conjugate pairs together as quadratics, and turns the decomposition into a recipe for inverting Laplace and z-transforms term by term.
An inverse Laplace transform is hard in general and easy for a handful of standard shapes: 1/(s−p), 1/(s−p)m, and the two real forms (s−α)/[(s−α)2+β2] and β/[(s−α)2+β2]. The point of a partial-fraction decomposition is that every rational transform is a finite sum of exactly those shapes, so the inversion is a table lookup once the algebra is done. The Bromwich integral of The Laplace transform is not needed; the decomposition is what the Bromwich contour would have produced anyway, pole by pole.
Everything here is algebra over C and applies verbatim to any transform whose inverse is read off a table: Laplace, z- and discrete-time transforms, and generating functions (Generating functions).
Definition(Proper rational functions, poles and singular parts)AdaptedAdapted from Ahlfors (1979), Ch. 2, §1.4, pp. 30–32; Riley, Hobson & Bence (2006), §1.4, p. 18 and footnote
A rational function is a quotient R=g/h of polynomials with complex coefficients, h≡0. We always assume the quotient is in lowest terms: g and h have no common root. Following Riley, Hobson and Bence we also normalise h to be monic, dividing both polynomials by the leading coefficient of h if necessary.
which the fundamental theorem of algebra always permits over C, the numbers pk are the poles of R and mk is the order of the pole pk. A pole of order 1 is simple.
The singular part of R at pk is a sum
Sk(s)=j=1∑mk(s−pk)jak,j,
that is, a polynomial in 1/(s−pk) with no constant term. Ahlfors calls this the singular part of R at pk and builds the decomposition out of these pieces, one for each pole and one more at ∞.
Let R=g/h be rational and in lowest terms, with the factorisation of h above.
1. Existence and uniqueness. There is exactly one polynomial G and exactly one family of complex numbers ak,j (1≤k≤q, 1≤j≤mk) with
R(s)=G(s)+k=1∑qj=1∑mk(s−pk)jak,j.
The polynomial G is the quotient of g by h: it is zero exactly when R is proper, and otherwise has degree degg−degh. The number of unknown coefficients ak,j is degh, one per root counted with multiplicity.
Two cases are worth naming. The cover-up rule is the case j=mk, in which no derivative is taken:
ak,mk=s→pklim(s−pk)mkR(s),
so for a simple pole ak,1=lims→pk(s−pk)R(s), which is Arfken's Eq. (1.84) and, in Riley, Hobson and Bence's wording, the instruction to delete the factor s−pk from h and evaluate what is left at s=pk. When every pole is simple this can be written with the derivative of the denominator,
R(s)=G(s)+k=1∑qh′(pk)g(pk)s−pk1,
which is Ahlfors's Exercise 2 on p. 32. Each ak,1 is the residue of R at pk in the sense of The residue theorem.
3. The real form. If g and h have real coefficients then ak,j and al,j are complex conjugates whenever pl=pk, and the two singular parts of a conjugate pair of poles of order m combine into
j=1∑m[(s−α)2+β2]jBjs+Cj,pk=α+iβ,Bj,Cj∈R,
so a real R has a decomposition with only real constants, in which each irreducible real quadratic factor of h carries a numerator of degree at most one. Riley, Hobson and Bence give this real form for factors x2+ai2 (their p. 23), and Arfken for a general quadratic x2+px+q (his p. 42).
Method(Heaviside expansion theorem)AdaptedAdapted from Arfken, Weber & Harris (2013), §20.10, Exercise 20.10.11, p. 1045
If a Laplace transform can be written as F(s)=g(s)/h(s) with g and h analytic and h having simple, isolated zeros at the points si, then
f(t)=L−1{F}(t)=i∑h′(si)g(si)esit.
Arfken states this as an exercise and in this generality: g and h need only be analytic, not polynomials, so the sum may be infinite — the transform 1/(scoshs) of a heat-conduction problem is of this type. For polynomial g and h with degg<degh it is exactly part 2 of the decomposition theorem above, inverted term by term with 1/(s−p)↦ept.
Idea of the proofOriginal
A rational function is determined, up to a polynomial, by how it blows up. Near a pole pk of order mk the function behaves like a finite string of powers (s−pk)−mk,…,(s−pk)−1, and that string is the singular part Sk. Subtract all of the singular parts: what remains is rational with no poles anywhere, hence a polynomial. That is the whole of existence, and it is Ahlfors's argument.
Uniqueness is a counting statement: two decompositions would differ by a rational function that is identically zero, and multiplying the difference by (s−pk)mk and letting s→pk kills every term except one, forcing the coefficients to agree. Running that same limit on R itself, rather than on the difference, is exactly the formula for the coefficients — which is why the cover-up rule is not a trick but the uniqueness proof read forwards.
Proof(Existence, uniqueness and the coefficient formula)AdaptedAdapted from Ahlfors (1979), Ch. 2, §1.4, Eqs. (12)–(14), pp. 31–32
Step 1: isolate the behaviour at one pole. Fix k and put Φk(s)=(s−pk)mkR(s). Because the quotient is in lowest terms, g(pk)=0, and the factor (s−pk)mk cancels the whole of that factor in h; so Φk is rational with no pole at pk, and Φk(pk)=0. Hence Φk is differentiable at pk to all orders and has a Taylor expansion there (Taylor series of analytic functions),
Φk(s)=i≥0∑ck,i(s−pk)i,ck,i=i!Φk(i)(pk).
Dividing by (s−pk)mk,
R(s)==:Sk(s)i=0∑mk−1(s−pk)mk−ick,i+no pole at pki≥mk∑ck,i(s−pk)i−mk.
Setting j=mk−i identifies the coefficients of Sk as ak,j=ck,mk−j=Φk(mk−j)(pk)/(mk−j)!, which is the formula in part 2 of the statement. For j=mk the derivative order is zero and the formula reads ak,mk=Φk(pk): the cover-up rule.
Step 2: subtract every singular part. Consider
D(s)=R(s)−k=1∑qSk(s).
D is rational, so its only possible poles are the pk. At pk, Step 1 shows R−Sk has no pole, and the other Sl (l=k) are finite at pk because pl=pk; so D has no pole at pk either. A rational function with no poles in C has a denominator with no roots, hence a constant denominator: D is a polynomial, which we name G. This gives the decomposition.
Step 3: G is the quotient of the division. Each Sk(s)→0 as ∣s∣→∞, so R−G→0. Multiplying the decomposition by h gives g=Gh+r with r=h⋅(R−G) a polynomial of degree less than degh — because r/h→0 at infinity. By uniqueness of polynomial division, G is the quotient and r the remainder of g by h. In particular G=0 precisely when R is proper, and degG=degg−degh otherwise. This is the preliminary long division of Riley, Hobson and Bence's Eq. (1.45).
Step 4: uniqueness. Suppose two decompositions of R exist, with data (G,ak,j) and (G~,a~k,j). Subtracting, the rational function
0=(G−G~)(s)+k∑j∑(s−pk)jak,j−a~k,j
holds for all s outside the poles. Fix k, multiply by (s−pk)mk and let s→pk: every term with l=k, and the polynomial term, carries the factor (s−pk)mk→0; the terms with j<mk carry (s−pk)mk−j→0; only ak,mk−a~k,mk survives, so it is zero. Remove that term and repeat with (s−pk)mk−1 to get ak,mk−1=a~k,mk−1, and so on down to j=1. Doing this at every pole leaves G−G~=0. The count ∑kmk=degh of coefficients is immediate from the factorisation.
Step 5: the real form. Let g,h have real coefficients. Taking complex conjugates in R(sˉ)=R(s) and using uniqueness, the decomposition of R obtained by conjugating all coefficients and swapping each pk for pk is again a decomposition of R; so akˉ,j=ak,j, where kˉ indexes the pole pk. For a conjugate pair p=α+iβ, pˉ=α−iβ with β=0, the two jth terms add to
using (s−p)(s−pˉ)=(s−α)2+β2. The numerator is 2Re[aj(s−pˉ)j]=:Nj(s), a real polynomial of degree ≤j, and we write Q(s)=(s−α)2+β2 for the real quadratic in the denominator. It remains to reduce Nj/Qj to the advertised shape. Divide Nj by Q with remainder: Nj=QUj+(Bjs+Cj) with Bj,Cj∈R and Uj real of degree ≤j−2, so that
QjNj=QjBjs+Cj+Qj−1Uj.
The second term has the same shape as the first with the exponent lowered by one and the numerator degree lowered by two, so the division can be repeated on it. After r repetitions the exponent is j−r and the numerator degree is at most j−2r; by the time the exponent has fallen to 1, at r=j−1, the numerator degree is at most 2−j≤1 and nothing is left to divide. So Nj/Qj is a real combination of terms (Bs+C)/Qi with i≤j. Summing over j=1,…,m and collecting the terms with equal i gives the form in part 3. □
Remark(Where the hypotheses are used)Original
Lowest terms. Step 1 uses g(pk)=0 only to say that pk really is a pole of order mk. If g and h share the root pk, the cover-up rule applied to the uncancelled h returns lims→pk(s−pk)mkg(s)/h(s)=0, because the numerator still carries the factor s−pk; and the ansatz of "one term per factor of h" asks for more terms than the function has poles. Cancel common factors first.
Properness. Nothing fails without it — G simply is not zero. What fails is the naive ansatz: (s3+2)/(s2−1) is in lowest terms (the numerator is 3 at s=1 and 1 at s=−1), and if one writes it as A/(s−1)+B/(s+1) the right side tends to 0 as ∣s∣→∞ while the left side grows like s, so no choice of A,B can repair it. Riley, Hobson and Bence make the division the first step for exactly this reason.
Distinctness of the pk. Step 2 needs Sl to be finite at pk for l=k, which is the statement pl=pk. This is why the factorisation must be organised by distinct roots with multiplicities, not by repeating a factor.
The full string of powers at a repeated pole. Dropping the lower powers is the standard error. Take R(s)=1/[(s+1)(s−3)2] and suppose it could be written as A/(s+1)+B/(s−3)2, with no 1/(s−3) term. Clearing denominators, 1=A(s−3)2+B(s+1) would have to hold for every s. Putting s=3 gives B=41 and putting s=−1 gives A=161; but then the coefficient of s2 is A=161 on the right and 0 on the left. No pair (A,B) works, and the missing 1/(s−3) term is what repairs it. Riley, Hobson and Bence issue the same warning on their p. 23.
The cover-up rule at a repeated pole. With mk>1 the rule gives only ak,mk, the coefficient of the highest inverse power; the remaining ak,j need the derivatives. Riley, Hobson and Bence state this restriction explicitly on their p. 24. The reason is visible in Step 1: a derivative of Φk is needed to reach the next Taylor coefficient.
Real coefficients for the real form. Step 5 is the only place realness is used. For a genuinely complex R there is no real form, and one works with the complex poles directly.
β=0 in the real form. A real quadratic with β=0 is (s−α)2, not irreducible; it belongs to the repeated-real-pole case.
Cancel and normalise. Remove common factors of g and h; divide through so that h is monic.
Divide if improper. If degg≥degh, divide to get F=G+r/h with degr<degh. For a Laplace transform a non-zero G signals a delta function and its derivatives in F(t), since L{δ}=1; a transform of an ordinary function is always proper.
Factor h and list the distinct poles pk with their orders mk. Decide now whether to work over C (every factor linear) or over R (complex pairs kept as quadratics). For hand inversion of a real transform the real route is usually shorter, because it lands directly on eαtcosβt and eαtsinβt.
Write the ansatz with one term per factor: ak,j/(s−pk)j for j=1,…,mk at each real or complex pole, and (Bjs+Cj)/[(s−α)2+β2]j for each conjugate pair kept real.
Fix the coefficients. Cover-up first: it gives every ak,1 at a simple pole and the top coefficient ak,mk at a repeated one, with no simultaneous equations. For the rest, either differentiate as in part 2 of the statement, or subtract the terms already known from F and simplify, or substitute a few convenient values of s and solve. Riley, Hobson and Bence list these as methods (i)–(iii) and mix them freely.
Check. Evaluate both sides at one value of s that is not a pole — s=0 is usually easiest — or compare the coefficient of 1/s as s→∞: for a proper F, ∑kak,1 must equal lims→∞sF(s).
Invert term by term, using linearity and the standard pairs of The Laplace transform:
To use the last two, first write the numerator of a quadratic term as a multiple of s−α plus a constant; completing the square in the denominator supplies α and β. This is the step Arfken's Example 20.7.2 performs, reading each piece off a table.
8. Check the answer at t=0. The initial-value theorem f(0+)=lims→∞sF(s) catches most sign errors in one line.
For a z-transformX(z)=∑n≥0xnz−n, apply the same steps to X(z)/z rather than to X(z), then multiply each resulting term by z. The reason is that the useful tabulated pair is z/(z−p)↦pn, not 1/(z−p); decomposing X(z)/z puts a factor z on every term automatically and keeps the n=0 value right. A repeated pole gives z/(z−p)2↦(n/p)pn, and a conjugate pair gives a damped oscillation rncos(nθ+φ) with p=reiθ.
Example(A fourth-order transform, and a z-transform)Original
(a) Laplace, with a double pole and a complex pair. Invert
F(s)=(s+2)2(s2+2s+5)2s2+11s+19.
Step 1–3. The numerator has no root at s=−2 (8−22+19=5=0) nor at the roots of s2+2s+5, so the quotient is in lowest terms; it is proper (degree 2 over degree 4), so G=0. The poles are s=−2, double, and s=−1±2i, simple. We keep the pair real: s2+2s+5=(s+1)2+4, so α=−1, β=2.
The numerator is 2s2+11s+19−(s2+2s+5)−(s3+4s2+9s+10)=−s3−3s2+4=−(s+2)2(s−1), so the quotient is (−s+1)/(s2+2s+5) and C=−1, D=1.
Step 6: check. At s=0 the left side is 19/(4⋅5)=0.95 and the right side is 41+21+51=0.95. Also lims→∞sF(s)=0 and the residues sum to B+C=1−1=0, as they must.
Step 7. Split the quadratic term as (−s+1)=−(s+1)+2:
which gives x1=1, x2=1.5, x3=1.75, x4=1.875 — the same as 2−21−n. The sequence climbs to its limit x∞=2, and the rate at which it does so is set by the pole at z=21: the transient is −2⋅2−n, halving each step.
Application(Where the decomposition does the work)Original
Modes of a linear system. Laplace-transforming a constant-coefficient ODE (Solving initial-value problems with the Laplace transform) turns it into P(s)X(s)=Q(s), with P the characteristic polynomial and Q carrying the forcing and the initial data. The response is the rational function X=Q/P, and its partial-fraction decomposition is the modal decomposition of the system: one term per root of P, each inverting to a mode ept, tjept or eαtcosβt. Stability is read off the poles — all in Res<0 — and the slowest mode is the pole nearest the imaginary axis, as in part (a), where the pair at −1±2i decays more slowly than the double pole at −2 and therefore dominates for large t. A repeated pole is the critically damped case, and the factor t in te−2t is where it shows itself.
Circuits and transfer functions. A network of resistors, capacitors and inductors has an impedance that is a rational function of s; so is the transfer function from any input to any output. Partial fractions convert it into a sum of first-order lags and resonant second-order sections, which is both how such a filter is analysed and how it is physically built, section by section. The coefficient ak,1 at a pole is the residue, so the same numbers appear if one instead closes the Bromwich contour of The Laplace transform and sums residues.
Digital filters and difference equations. A linear recurrence ∑jαjyn−j=∑jβjun−j becomes Y(z)=H(z)U(z) with H rational, and part (b) is the pattern for getting the impulse response back: poles inside ∣z∣<1 give decaying terms, a pole on the unit circle gives a persistent oscillation or a constant, and a pole outside means the filter diverges. The same algebra on a probability generating function (Generating functions) turns a rational generating function into a closed form for the nth term — the Fibonacci numbers being the smallest example, with two simple poles.
Rational integrands and series. The decomposition is also what makes ∫R(x)dx elementary for every rational R: each term integrates to a logarithm, a power, or an arctangent. Arfken uses it in the opposite direction, to convert an awkward series into telescoping pieces, and it is the finite counterpart of the pole expansions of a meromorphic function in Pole expansions and infinite products.
Defines the Hankel transform of order ν and the Mellin transform in the normalisation of the cited books, shows that each is the Fourier transform seen through a change of variable — Hankel is the plane Fourier transform restricted to radial functions, Mellin is the line Fourier transform after t=ex — and reads the inversion formula for each off that identification.
The integral converges absolutely for every k. The reason is that Jν is bounded on [0,∞): it is continuous there, with Jν(0)=0 for ν>0 and J0(0)=1, and it is O(x−1/2) as x→∞ by Asymptotic forms of Bessel functions, so it tends to 0 at infinity. A continuous function on [0,∞) with limit 0 at infinity is bounded, say by Mν. Hence ∣f(r)Jν(kr)r∣≤Mν∣f(r)∣r, which is integrable by hypothesis.
The inversion formula has the same shape with the roles of r and k exchanged:
f(r)=∫0∞Fν(k)Jν(kr)kdk.
So Hν is, where it is valid, its own inverse. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.3, example (i), p. 459) list exactly this pair, and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.1, (20.8), p. 965; §21.2, (21.21), p. 1054) give the same forward transform and the same self-reciprocity. The statement and proof below establish the inversion formula for ν=0 on a class of functions for which every step can be justified.
Two remarks on the weight rdr. It is the Sturm–Liouville weight of Bessel's equation, the one that makes the Bessel functions orthogonal (Orthogonality and Fourier-Bessel series); Stakgold and Holst (Stakgold & Holst (2011), §7.3, (7.3.31)–(7.3.39), pp. 450–451) obtain the pair precisely as the spectral expansion of that singular problem on 0<x<∞, with the endpoint condition limx→0xu′=0 selecting the solution regular at the origin, and Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.1, p. 965) describe the transform as the continuum limit of a Bessel series. It is also, up to the factor 2π from the angular integration, the area element of the plane in polar coordinates, which is the content of the first theorem below.
Let f be measurable on (0,∞). The Mellin transform of f is
M[f](z)=F(z)=∫0∞f(t)tz−1dt,
defined at those complex z for which the integral converges absolutely. Here tz−1=e(z−1)lnt with the real logarithm of the positive number t, so no branch cut is involved. Since ∣tz−1∣=tσ−1 with σ=Rez, convergence depends only on σ. The set of good σ is an interval, and the reason is one inequality: if the integral converges absolutely at σ1 and at σ2>σ1, then for σ1≤σ≤σ2 we have tσ−1≤tσ1−1+tσ2−1 at every t>0 — the first term dominating on (0,1] and the second on [1,∞) — so it converges absolutely at σ as well. The domain of F is therefore a vertical strip. What the next statement adds is a sufficient condition, in terms of power bounds on f, for a prescribed open strip a<σ<b to lie inside that domain.
This unit uses throughout the inversion formula of Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §21.2, (21.20), p. 1054) and Hassani (Hassani (2013), §16.3, Problem 16.1, p. 505),
f(t)=2πi1∫c−i∞c+i∞F(z)t−zdz,
the integral running up a vertical line Rez=c lying inside the strip. Arfken writes the abscissa γ and Hassani writes it α; the value of c is free within the strip, and the two proofs below show both why it is free and why some choice must be made.
Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §20.1, p. 966) point out that the gamma function is a Mellin transform: M[e−t](z)=Γ(z) is exactly Euler's integral of The Gamma function, convergent for Rez>0.
The Hankel pair. All three books that write it out — Arfken, Weber and Harris (Arfken, Weber & Harris (2013), (20.8), p. 965), Riley, Hobson and Bence (Riley, Hobson & Bence (2006), example (i), p. 459) and, at order zero, Stakgold and Holst (Stakgold & Holst (2011), (7.3.38)–(7.3.39), p. 451) — agree on the normalisation used above, weight rdr forward and kdk back, with no numerical prefactor anywhere. Hassani (Hassani (2013), Example 16.0.1(e), p. 494) lists the kernel as yJn(xy), which is the same thing with the weight absorbed into the kernel. There is nothing to reconcile.
The Mellin inversion contour. Here the sources genuinely differ. Arfken, Weber and Harris (Arfken, Weber & Harris (2013), (21.20), p. 1054) and Hassani (Hassani (2013), Problem 16.1, p. 505) put the line at a general abscissa, as above. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), example (ii), p. 459) write the limits as −i∞ to i∞, that is, with the line through the origin. The two are not interchangeable, because the strip of convergence need not contain the imaginary axis: for f(t)=1/(1+t) the strip is 0<Rez<1, computed in the second example below, and the integral defining F diverges at z=0. This unit follows Arfken and Hassani and keeps the free abscissa c. Where the strip does straddle the imaginary axis and the inversion hypotheses below hold at c=0 as well, the two prescriptions give the same number, since both give f(t).
A variable on the imaginary axis.Stakgold and Holst (Stakgold & Holst (2011), Exercise 7.3.4, (7.3.52)–(7.3.53), p. 456) reach the pair from the singular spectral problem −(rR′)′−(λ/r)R=0 on 0<r<∞, with p(r)=r and s(r)=1/r, and write it with a real spectral variable ν in a purely imaginary exponent, FM(ν)=∫0∞f(r)r−iν−1dr. Setting z=−iν turns their forward formula into ours restricted to the line Rez=0, and their inversion integral over ν∈R into ours over that line: dν=idz, the limits reverse because z=−iν runs downwards as ν runs upwards, and the resulting constant −i/2π is 1/2πi. So theirs is the special case c=0 in a different variable, not a different transform.
The sign in the plane.Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 20.2.15, p. 979) write the two-dimensional Fourier pair with e+i(ux+vy) forward. As recorded in Fourier transform in several dimensions, the sign is invisible on radial functions, because their transforms are even; so their statement that the order-zero Hankel pair is a special case of the Fourier pair is the same statement as ours, proved below in the convention of that unit.
Let f:(0,∞)→C be such that g(ρ)=f(∣ρ∣) belongs to the Schwartz class S(R2), and write g~ for its two-dimensional Fourier transform in the symmetric convention of Fourier transform in several dimensions,
g~(k)=2π1∫R2g(ρ)e−ik⋅ρd2ρ.
Then:
Identification.g~ is radial and its profile is the order-zero Hankel transform of the profile of g:
g~(k)=F0(k)=∫0∞f(r)J0(kr)rdr,k=∣k∣.
No prefactor survives: in this one convention the two-dimensional Fourier transform of a radial function is the order-zero Hankel transform of its profile.
Inversion.F0 is itself the profile of a Schwartz function on R2 — because the Fourier transform maps S(Rd) into itself, which is Stein and Shakarchi's (Stein & Shakarchi (2003), Ch. 6, §2, Cor. 2.2, p. 181) Corollary 2.2 — and for every r>0f(r)=∫0∞F0(k)J0(kr)kdk.
Thus H02=I on this class: the order-zero Hankel transform is an involution.
Parseval.∫0∞∣f(r)∣2rdr=∫0∞∣F0(k)∣2kdk.
Three dimensions, order 21. If instead g(r)=f(∣r∣) lies in S(R3), then with the same convention
g~(k)=k−1/2H1/2[r1/2f(r)](k).Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Exercise 20.2.15 and its Note, p. 979) record that the construction continues through ν=0,21,1,23,… — the two cases proved here, n=2 giving ν=0 and n=3 giving ν=21, fit the pattern ν=2n−1, which is our reading of it and not something they assert — and that at the two half-integral orders ν=±21 what results is a Fourier sine or cosine transform. Their ν=−21 lies outside the range ν≥0 fixed in the definition above, so of the two orders they name only ν=+21 is a Hankel transform in the sense defined here — it is the one in part 4, and the sine transform is the one that appears in its proof. (The cosine transform is the companion statement, since J−1/2(x)=2/πxcosx; but H−1/2 is not defined in this unit, and nothing below uses it.) They refer to Sneddon for the general case, and this unit proves only the two cases n=2 and n=3 displayed above.
Definition(The strip, and Mellin inversion)Original
Let f be measurable on (0,∞) and suppose there are real numbers a<b and a constant C with
∣f(t)∣≤Ct−a(0<t≤1),∣f(t)∣≤Ct−b(t≥1).
The strip. These bounds make the integral defining F(z)=M[f](z) converge absolutely on the open strip a<Rez<b, and there F is holomorphic, with
F′(z)=∫0∞f(t)(lnt)tz−1dt.
It is a Fourier transform. Fix c with a<c<b and put g(x)=ecxf(ex) for x∈R. Then g∈L1(R) and, in the convention g^(k)=∫−∞∞g(x)e−ikxdx of The Fourier transform,
F(c+iξ)=g^(−ξ)(ξ∈R).
Inversion. If in addition f is continuous on (0,∞) and ξ↦F(c+iξ) is integrable on R, then for every t>0f(t)=2πi1∫c−i∞c+i∞F(z)t−zdz=2π1∫−∞∞F(c+iξ)t−c−iξdξ.
The value does not depend on which admissible c is used: if the hypotheses hold at two abscissas, both integrals equal f(t) and so equal each other.
This is the statement Hassani's (Hassani (2013), §16.3, Problem 16.1, p. 505) substitution k=lnt asks for, with hypotheses supplied; the decay bounds and the holomorphy are ours.
Definition(Operational rules for the Mellin transform)Original
Let F=M[f] with strip a<Rez<b. For λ>0, p>0 and μ∈R:
Scaling.M[f(λt)](z)=λ−zF(z), on the same strip.
Weighting.M[tμf(t)](z)=F(z+μ), on the strip shifted by −μ.
Powers.M[f(tp)](z)=p−1F(z/p), on the strip pa<Rez<pb.
Reciprocal.M[f(1/t)](z)=F(−z), on the strip −b<Rez<−a.
The Euler operator. Write D=td/dt. If f is continuously differentiable and tzf(t)→0 as t→0+ and as t→∞, then
M[Df](z)=−zF(z),henceM[Dnf](z)=(−z)nF(z)
whenever the corresponding limits vanish at each stage.
Rule 5 is the reason the transform exists: it converts the operator D, which generates dilations t↦λt, into multiplication by −z, exactly as the Fourier transform converts d/dx into multiplication by ik. Rules 1 and 5 together say that the Mellin transform diagonalises dilations of the half-line, as the Fourier transform diagonalises translations of the line.
Idea of the proofOriginal
Neither transform is new. Each is the Fourier transform read in coordinates adapted to a symmetry, and each inversion formula is Fourier inversion carried back through the change of variable.
For Mellin the symmetry is dilation. The map t=ex carries the multiplicative half-line (0,∞) onto the additive line R, and under it the Mellin kernel together with its measure, tz−1dt, becomes ezxdx. On the vertical line z=c+iξ that is a Fourier kernel eiξx times a fixed real weight ecx, so the Mellin transform on that line is the Fourier transform of ecxf(ex). The free abscissa c is not an artefact: it is the exponential weight one is allowed to insert to make f(ex) integrable, and different admissible weights give the same answer because each of them returns f.
For Hankel the symmetry is rotation about an axis. The Fourier transform on R2 commutes with rotations, so it maps radial functions to radial functions; and when the angular integral in the plane is done, the plane-wave kernel collapses to J0. The two-dimensional pair is symmetric between ρ and k, so the resulting one-dimensional pair is symmetric between r and k — which is precisely the self-reciprocity of H0. The weight rdr is the polar area element with the angle integrated out, which is why it appears on both sides.
Step 1: the angular integral. Let g(ρ)=f(∣ρ∣) lie in S(R2). Two facts from Fourier transform in several dimensions apply at once: the transform of a radial Schwartz function is radial, and in the plane the radial reduction reads
g~(k)=2π1∫0∞2πρf(ρ)J0(kρ)dρ.
Its hypothesis, ∫0∞ρ∣f(ρ)∣dρ<∞, holds because a Schwartz function decays faster than any power. The two factors of 2π — one from the symmetric normalisation, one from integrating e−ikρcosφ over φ — cancel, leaving
g~(k)=∫0∞f(ρ)J0(kρ)ρdρ=F0(k),
which is part 1. The cancellation is special to n=2 and to this normalisation; it is what makes the order-zero pair as clean as it is.
Step 2: inversion. Because g∈S(R2), so is g~. This is the one ingredient not among the rules of Fourier transform in several dimensions, and it is Stein and Shakarchi's (Stein & Shakarchi (2003), Ch. 6, §2, Cor. 2.2, p. 181) Corollary 2.2, that the Fourier transform maps S(Rd) to itself; they deduce it from their Proposition 2.1 (iv)–(v), which are items 3 of the rules in our unit, since those turn kα∂kβg~ into a constant times the transform of ∂α(ρβg) and a transform is bounded by the L1 norm of what it transforms. Their convention differs from the one used here only by the fixed factor (2π)−n/2 and the dilation ξ=k/2π of the transform variable, and S(R2) is closed under both. Two consequences are used below: ∫0∞k∣g~(k)∣dk<∞, which is the hypothesis step 1 needs before it can be applied to g~ in place of g; and the first clause of part 2.
The right-hand side is the forward transform of g~ evaluated at −ρ. Now g~ is radial, so step 1 applies to it word for word, and it says that the forward transform of g~ at any point of length r equals ∫0∞g~(k)J0(kr)kdk. Since ∣−ρ∣=∣ρ∣=r, that is the value of the right-hand side. Since the left-hand side is f(r) and g~(k)=F0(k), this is the displayed formula of part 2. Its first clause, that F0 is the profile of a Schwartz function, is the statement g~∈S(R2) settled at the start of this step by Stein and Shakarchi's corollary; it is an input to the argument, not a consequence of the inversion formula just proved.
Step 3: Parseval. Plancherel's identity in the plane, item 6 of the rules in Fourier transform in several dimensions, says ∫R2∣g~∣2d2k=∫R2∣g∣2d2ρ. Both sides are integrals of radial functions, so each reduces by the polar area element: ∫R2∣g∣2d2ρ=2π∫0∞∣f(r)∣2rdr, and likewise on the other side with F0. Dividing by 2π gives part 3.
Step 4: three dimensions. For g∈S(R3) radial, the same unit gives the reduction
g~(k)=(2π)3/21k4π∫0∞rf(r)sin(kr)dr.
Collect the constants: 4π/(2π)3/2=2/(2π)1/2=2/π. Now use the closed form J1/2(x)=2/πxsinx of Spherical Bessel functions, which gives sin(kr)=πkr/2J1/2(kr). Substituting,
since 2/π⋅πk/2=k and one factor r1/2 is left over after writing the measure as rdr. The bracket is the function whose order-21 Hankel transform is taken, which is part 4. □
Proof(The Mellin theorems)AdaptedAdapted from Hassani (2013), §16.3, Problem 16.1, p. 505
Step 1: the strip. Let σ=Rez. Split the integral at t=1 and use the two bounds:
The first right-hand integral is finite exactly when σ−a>0, the second exactly when σ−b<0. So both bounds are finite precisely on a<σ<b, and there the defining integral converges absolutely. Outside that range the bounds say nothing; a particular f may of course decay better than its bounds, which is why the strip is only guaranteed, not determined, by a and b. What the hypothesis a<b supplies is that the range is not empty, and emptiness really does happen: for f(t)=1/t, where the two bounds are attained with a=b=1, the integral is ∫0∞tσ−2dt, which diverges for every σ.
Step 2: holomorphy. The Leibniz rule for differentiating under the integral sign (Differentiation under the integral sign) is not available here, and exhibiting a dominating function does not make it so: that rule is stated for a bounded interval, a real parameter and a real integrand continuous together with its partial derivatives, whereas our interval is (0,∞), our parameter is complex, and the integrand is unbounded at t=0 whenever a>0. The Gamma function records exactly this obstruction for Euler's integral — which is the case f(t)=e−t of the present one, as noted above — and the route taken there is the route taken here: Morera's theorem.
A dominating function on each closed substrip. Fix a<σ1≤σ2<b and let S={σ1≤Rez≤σ2}. For z∈S and t>0,
the first power dominating on (0,1] and the second on [1,∞). Step 1, applied at the two abscissas σ1 and σ2, says h is integrable. Note that h does not depend on z.
Continuity. If zm→z inside S then f(t)tzm−1→f(t)tz−1 for each fixed t, because z↦e(z−1)lnt is continuous, and every term is bounded by h. So Dominated convergence gives F(zm)→F(z). Every point of the open strip is interior to some such S, so F is continuous on the strip.
Closed-contour integrals vanish. Let γ be a closed contour in the strip. Its image is compact, hence contained in some S as above, and
each inner integral vanishing by Cauchy's theorem, since for fixed t>0 the function z↦tz−1=e(z−1)lnt is entire. A continuous function all of whose closed-contour integrals vanish is holomorphic, by Morera's theorem (part 4 of the consequences in Cauchy's integral formula). Hence F is holomorphic on the strip. Morera delivers holomorphy outright, so no separate appeal to the Cauchy–Riemann equations is needed, and none is made.
The derivative. Fix z in the strip and a circle C={∣w−z∣=r} whose closed disc lies in the strip. Cauchy's formula for the first derivative (Cauchy's integral formula) gives F′(z)=(2πi)−1∮CF(w)(w−z)−2dw, and the exchange above is available again — the extra factor ∣w−z∣−2 is the constant r−2 on C — so
the inner integral being the same Cauchy formula applied to the entire function w↦tw−1, whose derivative at z is (lnt)tz−1. That the final integral converges absolutely is the one place the logarithm has to be absorbed into the powers. For u≥1 and ε>0 one has lnu≤uε/ε: the difference uε/ε−lnu equals 1/ε>0 at u=1 and has derivative uε−1−u−1=(uε−1)/u≥0 there and beyond; applying it to u=1/t and to u=t gives ∣lnt∣≤ε−1t−ε on (0,1] and ∣lnt∣≤ε−1tε on [1,∞). Choose ε small enough that a<σ1−ε and σ2+ε<b; then on S
∣f(t)(lnt)tz−1∣≤ε−1∣f(t)∣(tσ1−ε−1+tσ2+ε−1),
which step 1, applied now at the abscissas σ1−ε and σ2+ε, makes integrable. This is part 1 of the statement.
Step 3: the substitution. Put t=ex, so dt=exdx and x runs over all of R as t runs over (0,∞). This is a smooth increasing bijection, so for z=c+iξ,
By step 1 at σ=c this converges absolutely, which says exactly that g∈L1(R). Comparing with g^(k)=∫g(x)e−ikxdx gives F(c+iξ)=g^(−ξ), which is part 2 of the statement. Note which way the sign goes: the Mellin variable runs up the line while the Fourier variable runs down it.
Step 4: inversion. Suppose further that f is continuous and ξ↦F(c+iξ) is integrable. Then g∈L1 and g^∈L1, so the integrable version of Fourier inversion applies: g agrees almost everywhere with a continuous function given by the inversion integral. But g(x)=ecxf(ex) is already continuous, and two continuous functions equal almost everywhere are equal everywhere, so for every x
the second form obtained by the substitution k=−ξ, which reverses the limits and the sign of the measure, leaving the integral over R unchanged. Multiply by e−cx and move it inside:
f(ex)=2π1∫−∞∞F(c+iξ)e−(c+iξ)xdξ.
Finally write t=ex, so e−(c+iξ)x=t−z with z=c+iξ, and parametrise the vertical line by z: along it dz=idξ, so dξ=dz/i and the limits ξ=∓∞ become z=c∓i∞. This turns the display into f(t)=(2πi)−1∫c−i∞c+i∞F(z)t−zdz, which is part 3. Independence of c needs no further argument: if the hypotheses hold at two abscissas in the strip, each of the two integrals equals f(t). (Moving the contour directly, by applying Cauchy's theorem to a long rectangle, is a different and stronger statement: it needs F(σ+iT)→0 as ∣T∣→∞ uniformly for σ between the two abscissas, so that the two horizontal sides contribute nothing. That is the kind of hypothesis checked by hand in the second example below.)
Step 5: the operational rules. Each is one substitution.
Scaling.∫0∞f(λt)tz−1dt=λ−z∫0∞f(u)uz−1du after u=λt, since tz−1dt=λ−zuz−1du.
Weighting.∫0∞tμf(t)tz−1dt=∫0∞f(t)t(z+μ)−1dt, which converges whenever z+μ lies in the original strip.
Powers.u=tp gives t=u1/p and tz−1dt=p−1uz/p−1du, whence M[f(tp)](z)=p−1F(z/p).
Reciprocal.u=1/t gives dt=−u−2du and reverses the limits, so ∫0∞f(1/t)tz−1dt=∫0∞f(u)u−z−1du=F(−z).
The Euler operator.Df=tf′, so M[Df](z)=∫0∞f′(t)tzdt. Integrate by parts: the boundary term is [f(t)tz]0∞, which vanishes by hypothesis, and what is left is −z∫0∞f(t)tz−1dt=−zF(z). Applying this to Df in place of f gives M[D2f]=−zM[Df]=z2F, and induction gives (−z)nF. □
Remark(Where the hypotheses are used)Original
a<b is used in step 1 of the Mellin proof, and is the one condition that can fail outright. The function f(t)=1/t has a=b=1: the strip is empty and f has no Mellin transform at any z. The transform exists only for functions whose behaviour at the two ends leaves a gap between the two exponents, and the width of that gap is the whole freedom in the inversion contour.
Integrability of F along the line is used only in step 4, to invoke the L1 inversion theorem. Without it the inversion integral must be read as a symmetric principal value, exactly as in the corresponding Fourier statement of Fourier inversion; the formula is not thereby false, but it is a different theorem.
Continuity of f is used in step 4 to upgrade "almost everywhere" to "everywhere". It cannot be dropped: changing f at one point changes neither F nor the right-hand side, so the pointwise formula must fail at that point.
The vanishing boundary term in rule 5 is a real restriction. For f(t)=e−t we have tzf(t)→0 at both ends when Rez>0, so rule 5 applies. Here Df=tf′=−te−t, whose Mellin transform is −Γ(z+1) by rule 2, while rule 5 says it is −zΓ(z). Equating the two gives Γ(z+1)=zΓ(z): the functional equation of The Gamma function is rule 5 in disguise. For f≡1 the boundary term never vanishes, and indeed f has no strip at all.
The Schwartz hypothesis in the Hankel theorem is used twice, and more strongly than necessary. Step 1 needs only ∫0∞ρ∣f∣dρ<∞, which is what makes F0 defined. Step 2 needs the transform to be integrable too, so that the inverse Fourier integral converges absolutely; Schwartz decay is the cheapest hypothesis that delivers both, and it is preserved by the transform — Corollary 2.2 of Stein and Shakarchi (Stein & Shakarchi (2003), Ch. 6, §2, Cor. 2.2, p. 181), used in step 2 — which is what makes the involution statement clean. Part 3 says that H0 preserves the norm of L2((0,∞),rdr), but it says this only on the profiles of radial Schwartz functions. Promoting it to an isometry of the whole space is a further step, in two parts: those profiles must be dense in L2((0,∞),rdr), and a norm-preserving map defined on a dense subspace must be extended to the closure. Neither is carried out here. Granted them, the inversion formula holds on L2 in the mean-square sense rather than pointwise.
Radial symmetry is used in step 1 and again in step 2. In step 2 it is applied to g~, which is legitimate only because step 2 of the proof in Fourier transform in several dimensions — the step after the orthogonal-commutation lemma, and deduced from it — shows the transform of a radial function is radial. Without that, the inverse transform would not reduce to a one-dimensional integral and there would be no inversion formula of Hankel type.
The cancellation of 2π in step 1 depends on the convention, not on the mathematics. In the convention f^(k)=∫fe−ik⋅xd2x of The Fourier transform the same calculation gives g^(k)=2πF0(k), and the inversion formula then carries a compensating 1/(2π). The Hankel pair itself is unaffected: it has no free constant, because the forward and inverse formulas are the same formula.
Method(Choosing and using one of these transforms)AdaptedAdapted from Stakgold & Holst (2011), §7.3, (7.3.40)–(7.3.43), pp. 451–452, and Exercise 7.3.4(d), p. 456
Match the transform to the symmetry of the domain. A problem on 0<r<∞ that is invariant under rotations about an axis — the radial part of a cylindrical or plane-polar problem — calls for a Hankel transform. A problem on 0<r<∞ that is invariant under dilations r↦λr — a wedge, a cone, an equation whose radial part is of Euler type — calls for a Mellin transform. The test is which operator the transform diagonalises: r1∂r(r∂r)−ν2/r2 for Hankel, r∂r for Mellin.
Fix the order, or the strip. For Hankel, the order ν is set by the angular dependence: separating the Helmholtz equation in plane polars with an angular factor eiνφ leaves Bessel's equation of order ν for the radial factor, so a term eiνφ in the data forces Hν. For Mellin, read the strip off the behaviour of the data at r→0 and r→∞ using step 1 of the proof, and keep c inside it.
Transform the equation. The two key identities are these.
Hankel. If u and u′ decay fast enough at both ends that the boundary terms below vanish, then
Hν[r1drd(rdrdu)−r2ν2u](k)=−k2Uν(k).
This is two integrations by parts, with boundary term [ru′Jν(kr)−rudrdJν(kr)]0∞, followed by the observation that Jν(kr) satisfies Bessel's equation in the form (ry′)′−(ν2/r)y=−k2ry, so the surviving integral is −k2∫0∞uJν(kr)rdr. Stakgold and Holst (Stakgold & Holst (2011), §7.3, (7.3.41), p. 452) perform exactly this step on the ring-source problem.
Mellin. Multiply the equation by r2 first. In plane polars r2∇2u=D2u+∂φ2u with D=r∂r, so applying ∫0∞(⋅)rs−1dr and rule 5 turns ∇2u=0 into the ordinary differential equation
∂φ2∂2U(s,φ)+s2U(s,φ)=0,U(s,φ)=∫0∞u(r,φ)rs−1dr.
Solve the transformed problem. It is an ordinary differential equation in the remaining variable, or an algebraic equation if only one variable is left. Impose the boundary conditions, transformed the same way.
Invert. For Hankel, apply Hν again. For Mellin, evaluate (2πi)−1∫c−i∞c+i∞F(z)t−zdz — usually by closing the contour to one side and summing residues, which is legitimate when the integrand's contribution over the closing sides tends to zero; see The residue theorem and the second example below for how that is checked.
Check the strip at the end. The poles of the transformed solution tell you where the contour may be placed and hence which behaviour at r=0 and r=∞ the solution has. A pole crossed when the contour is shifted contributes a power of r to the answer; this is how Mellin inversion produces asymptotic expansions, discussed in Asymptotic series.
Example(The Gaussian is its own order-zero Hankel transform)Original
Take f(r)=e−αr2 with α>0. Then g(ρ)=e−α∣ρ∣2 is radial and Schwartz on R2, so the theorem applies and F0=g~. Compute g~ by factorising rather than by the radial integral.
Step 1: factorise. With ρ=(x,y) and k=(k1,k2),
e−α(x2+y2)=e−αx2e−αy2,
a product of one-dimensional factors, so by the separability rule — item 5 of Fourier transform in several dimensions — the two-dimensional transform is the product of the two one-dimensional ones.
Keep the normalisation in view, because it is the one place this calculation can go wrong. In item 5 each factor is itself a transform in the symmetric convention, so each carries its own (2π)−1/2, and the two together supply (2π)−1/2⋅(2π)−1/2=1/2π — the same 1/2π that the two-dimensional symmetric convention carries, as it must be. Steps 2 and 3 below therefore compute the two un-normalised integrals and apply that single factor 1/2π once, at the end. Dividing by 2π inside each factor and by 2π at the end would divide by 2π twice.
Step 2: one dimension, un-normalised. By Gaussian integrals, completing the square in the exponent gives
∫−∞∞e−αx2e−ik1xdx=απe−k12/4α.
Step 3: multiply and normalise. Multiplying the two un-normalised answers, one in k1 and one in k2, gives (π/α)e−(k12+k22)/4α; the two prefactors set aside in step 1 now divide this by 2π in total, and by nothing further, so with k2=k12+k22
F0(k)=g~(k)=2π1⋅απe−k2/4α=2α1e−k2/4α.
Step 4: read off the fixed point. The transform is again a Gaussian, with α replaced by 1/4α and multiplied by 1/2α. It equals the original exactly when 1/2α=1 and 1/4α=α, both of which say α=21. So
∫0∞e−r2/2J0(kr)rdr=e−k2/2,
and e−r2/2 is a fixed point of H0, as the Gaussian is of the Fourier transform. Consistency check on the involution: applying the formula of step 3 twice sends α↦1/4α with prefactor 1/2α, then 1/4α↦α with prefactor 1/(2⋅4α1)=2α, and the two prefactors multiply to 1 — which is part 2 of the theorem, verified by hand.
A number. At α=21 and k=2 the formula gives
∫0∞e−r2/2J0(2r)rdr=e−2=0.1353352832…,
and numerical quadrature of the left-hand side returns 0.135335283236612…, agreeing to every digit shown.
Example(A Mellin pair, and recovering the function by residues)Original
Take f(t)=1/(1+t) on (0,∞).
Step 1: the strip. As t→0+, f(t)→1, so ∣f(t)∣≤1=1⋅t0 on (0,1] and the exponent at the origin is a=0. As t→∞, f(t)≤1/t, so the exponent at infinity is b=1. By step 1 of the proof the strip is
0<Rez<1,
and it does not contain the imaginary axis — which is the point made in the note on conventions above.
Step 2: evaluate the transform. For 0<Rez<1 the integral ∫0∞tz−1(1+t)−1dt is the half-line form of Euler's Beta integral, part 4 of The Beta function at p=z and q=1−z: there p+q=1, so the denominator is (1+t)1, and the two hypotheses Rep>0 and Req>0 are precisely the two edges of our strip. So the integral is B(z,1−z)=Γ(z)Γ(1−z), and part 6 of the same unit evaluates that by the reflection formula of The Gamma function as π/sinπz — for complex z in the strip, which is what is needed here, and not merely for real exponents.
The keyhole-contour computation of Integrals around branch cuts returns the same value, but it is carried out there for a real exponent 0<p<1, so on its own it settles only the real segment of the strip. It reaches the rest by continuation: both sides are holomorphic on the strip, the left by part 1 of the strip theorem and the right because sinπz has no zero there, and two holomorphic functions agreeing on a segment agree on the connected open set containing it, by the identity theorem of Analytic continuation. Either way
F(z)=M[1+t1](z)=sinπzπ,0<Rez<1.
Observe the consistency: π/sinπz is holomorphic exactly on C∖Z, and its poles nearest the strip sit at z=0 and z=1, the two endpoints predicted by step 1. This is the usual picture: the ends of the strip are set by how f behaves at 0 and at ∞, and they are where the analytic continuation of F first becomes singular.
A number. At z=21 the formula gives F(21)=π/sin(π/2)=π, that is
∫0∞t(1+t)dt=π=3.14159265…,
which is checked in three lines without any of the above: substituting t=u2 turns the left-hand side into 2∫0∞du/(1+u2)=2⋅2π=π. Numerical quadrature returns 3.1415926535897932…
Step 3: invert, and get a series for free. Fix c with 0<c<1 and take t>1. In the inversion integral (2πi)−1∫c−i∞c+i∞πt−z(sinπz)−1dz the factor t−z=e−zlnt decays as Rez→+∞, because lnt>0; and 1/sinπz decays like e−π∣Imz∣ away from the real axis. So close the contour to the right. Take the rectangle with vertical sides Rez=c and Rez=N+21 and horizontal sides Imz=±T, and let T→∞ first and then N→∞ through the integers. On the horizontal sides ∣sinπz∣≥sinhπT, and the sides have length N+21−c, so their contribution is at most (N+21−c)πt−c/sinhπT, which tends to 0 as T→∞ with N fixed. On the right-hand side, ∣sinπ(N+21+iy)∣=coshπy, so that edge contributes at most
πt−(N+1/2)∫−∞∞coshπydy=πt−(N+1/2),
using ∫−∞∞sech(πy)dy=1; and this tends to 0 as N→∞ because t>1. For each fixed N the rectangle encloses the poles z=1,…,N and no others, so the residue theorem applies to it; traversed anticlockwise it runs down the side Rez=c, so letting T→∞ and then N→∞ leaves
the geometric series converging for t>1. The inversion has returned the function we started from, and in the form of its expansion at large t. Closing to the left instead picks up the poles at z=0,−1,−2,… with residues (−1)ntn and reproduces the expansion 1−t+t2−⋯ valid for t<1. This is the mechanism of step 6 of the method: which way the contour is closed decides which end of the range the resulting expansion describes.
Application(Axisymmetric potentials in cylindrical coordinates)AdaptedAdapted from Stakgold & Holst (2011), §7.3, (7.3.40)–(7.3.45), pp. 451–452
The Hankel transform is the natural tool for a field produced by a source with an axis of symmetry. Stakgold and Holst (Stakgold & Holst (2011), §7.3, (7.3.40)–(7.3.43), pp. 451–452) take steady heat conduction with absorption, −∇2u+k2u=source, driven by a unit ring of sources at radius a in the plane z=0. Because the source and the domain are invariant under rotation about the z axis, an order-zero Hankel transform in the cylindrical radius r removes that variable, exactly as in step 3 of the method; what is left is a one-dimensional problem in z, solved by requiring decay as z→±∞, and the answer is read back as an inverse Hankel transform.
The payoff of that calculation is a representation theorem. Letting the ring shrink to a point source at the origin, the solution must agree with the familiar screened potential e−kr2+z2/4πr2+z2, and comparing the two expressions gives what they attribute in essence to Sommerfeld: a spherically symmetric function is written as a superposition of cylindrical waves. At k=0 it reduces to
r2+z21=∫0∞J0(γr)e−γ∣z∣dγ,
which is their (Stakgold & Holst (2011), §7.3, (7.3.45), p. 452) formula (7.3.45). Read as a Hankel transform it says that the order-zero transform of γ↦e−γ∣z∣/γ is the Coulomb potential at height z above the plane; read physically it resolves a point charge into plane-parallel modes, which is what one wants near a conducting plate or a layered medium, where each mode can be reflected separately. As a check at r=1, z=21, both sides equal 2/5=0.8944271909999…
The same identity is the reason the Hankel transform appears in Green's function constructions for layered problems; compare the direct constructions of Green's functions for Laplace and Poisson problems.
Application(Diffraction by a circular aperture, and a wedge)Original
The Airy pattern. In the Fraunhofer regime the amplitude diffracted into a direction is proportional to the two-dimensional Fourier transform of the aperture function, as the one-dimensional version in The Fourier transform explains. For a circular hole of radius a the aperture function is radial — it is 1 for r<a and 0 beyond — so by part 1 of the Hankel theorem the amplitude is an order-zero Hankel transform:
F0(k)=∫0aJ0(kr)rdr=kaJ1(ka).
The evaluation is one line from the recurrence dxd[xJ1(x)]=xJ0(x) of Bessel functions of the first kind: substituting x=kr turns the integral into k−2∫0kaxJ0(x)dx=k−2[xJ1(x)]0ka, and xJ1(x) vanishes at 0.
The intensity is ∣F0∣2, that is the Airy pattern [2J1(ka)/(ka)]2 up to normalisation, and its first dark ring sits at the first positive zero of J1, namely ka=3.831705970… With k=(2π/λ)sinθ and aperture diameter D=2a this is
sinθ=π3.831705970⋅Dλ=1.2197Dλ,
the Rayleigh resolution limit of a telescope or a microscope objective, usually quoted as 1.22λ/D. The factor 1.22 is nothing but the first zero of J1 divided by π. Had the aperture been a slit rather than a hole, the same argument in one dimension would have given the first zero of sin, and the factor would have been 1: the extra 22% is the price of the circular geometry, and it enters through J1 rather than through any physics.
A wedge. The Mellin counterpart is the Laplace equation in a wedge 0<φ<α of infinite extent, with u=0 on one face and prescribed data h(r) on the other — the problem set as Exercise 7.3.4(d) by Stakgold and Holst (Stakgold & Holst (2011), §7.3, Exercise 7.3.4(d), p. 456). The domain is invariant under r↦λr and has no length scale, so by step 1 of the method the Mellin transform is the right one. Step 3 turns the problem into Uφφ+s2U=0 with U(s,0)=0 and U(s,α)=H(s), whose solution is
U(s,φ)=H(s)sinsαsinsφ,
and u is recovered by one inverse Mellin integral. The poles of 1/sinsα, at the points s=mπ/α with m∈Z, are the physics. Near the corner r is small, so ∣r−s∣=r−σ is small for very negative σ and the contour is pushed to the left; each pole it crosses, at s=−mπ/α with m=1,2,…, contributes a term in rmπ/α. The pole at s=0 contributes nothing, because the zero of sinsφ there cancels it. So the leading behaviour at the tip is rπ/α, set by the opening angle alone. For a re-entrant wedge, α>π, the exponent is less than 1 and the gradient blows up at the tip — the standard corner singularity that forces mesh refinement in numerical work and concentrates stress or electric field at a sharp edge.
Defines the discrete Fourier transform of N samples and its competing normalisations, proves the discrete orthogonality relation and from it the inversion formula, the discrete Parseval identity and the circular convolution theorem, derives the aliasing formula that says exactly which Fourier coefficients a finite sample confuses with which, and sets out the even-odd splitting that turns the O(N2) sum into the O(NlogN) radix-2 FFT.
Fix N∈N and write ζ=e2πi/N, a primitive Nth root of unity. A sample vector is a function F:Z→C with F(k+N)=F(k); equivalently an N-tuple (F(0),…,F(N−1)), with all indices read modulo N. The space V of such functions is N-dimensional, with the Hermitian inner product and norm
(F,G)=k=0∑N−1F(k)G(k),∥F∥2=(F,F).
The discrete exponentials are eℓ(k)=ζℓk=e2πiℓk/N for ℓ∈Z; only N of them are different, since eℓ+N=eℓ.
The discrete Fourier transform (DFT) of F is the sample vector F^ given by
F^(n)=N1k=0∑N−1F(k)e−2πink/N,n∈Z,
which is again N-periodic in n. This is the convention of Stein and Shakarchi, chosen here because it makes F^(0) the mean of the samples and because it is the one that matches the Fourier coefficients f^(n) of Fourier series in the limit of fine sampling; the statement below makes both of those precise.
Conventions. Three choices are free: the sign in the exponent, where the factor N goes, and whether the index runs over 0,…,N−1 or symmetrically.
Arfken's sign is the opposite of everyone else's, which conjugates the result; his and Hassani's symmetric N−1/2 makes the transform a unitary matrix, which is convenient in theory and a nuisance in bookkeeping. Stakgold and Holst use an odd number of symmetric indices, which avoids the awkward lone index at n=N/2. Their N is not ours: writing it NSH to keep them apart, their data sit at the 2NSH+1 points xk=2πk/(2NSH+1) with k=−NSH,…,NSH, so NSH is the highest harmonic index while their number of samples is 2NSH+1. Our N is their 2NSH+1. Converting between conventions is a matter of a factor N±1/2 and possibly a complex conjugation; nothing below depends on the choice except the constants we state.
Method(Orthogonality, inversion, Parseval and circular convolution)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 7, Lemma 1.1 and Theorem 1.2, pp. 222–223; Arfken, Weber & Harris (2013), §20.6, Eqs. (20.119)–(20.120) and (20.123)–(20.125), pp. 1003–1005; Stakgold & Holst (2011), §2.3, Eqs. (2.3.30)–(2.3.31), pp. 140–141
Arfken writes this as N∑r∈Zδm−ℓ,rN, which keeps the congruence visible; restricted to 0≤m,ℓ≤N−1 it is Nδmℓ. Consequently {N−1/2e0,…,N−1/2eN−1} is an orthonormal basis of V (Inner products and norms).
2. Inversion. For every k,
F(k)=n=0∑N−1F^(n)e2πink/N.
Equivalently: there is exactly one trigonometric polynomial with the N prescribed frequencies that takes the values F(k) at the N points xk=2πk/N, and its coefficients are the F^(n). This is how Stakgold and Holst introduce the transform, as the solution of an interpolation problem.
3. Parseval.N1k=0∑N−1∣F(k)∣2=n=0∑N−1∣F^(n)∣2, and more generally N1(F,G)=∑nF^(n)G^(n).
4. Shift. For j∈Z, the transform of k↦F(k−j) is n↦e−2πijn/NF^(n).
5. Circular convolution. Define (F∗G)(k)=∑j=0N−1F(j)G(k−j), the index k−j read modulo N. Then
F∗G(n)=NF^(n)G^(n).
The factor N is an artefact of our normalisation; with Arfken's symmetric convention, and his N−1/2 in front of the convolution sum as well, the rule is simply the product of the transforms, his Eq. (20.125).
6. Aliasing. Let f:R→C be continuous and 1-periodic, with an absolutely convergent Fourier series, ∑m∈Z∣f^(m)∣<∞. Sample it at N equally spaced points, F(k)=f(k/N). Then
So the DFT returns each Fourier coefficient contaminated by every coefficient whose frequency differs from it by a multiple of N: the sampling cannot tell those frequencies apart, because the corresponding exponentials agree at all N sample points. In particular, if f^(m)=0 for all ∣m∣≥N/2, then F^(n)=f^(n) for −N/2<n<N/2 and the DFT is exact.
Idea of the proofOriginal
Everything comes from one finite geometric sum. Put q=ζm−ℓ, which satisfies qN=1. Whenever q=1 the sum ∑k<Nqk collapses to (1−qN)/(1−q), and its numerator is zero. That single line gives orthogonality, and orthogonality gives everything else, exactly as on the circle: inversion is expansion in an orthogonal basis, Parseval is Pythagoras, and convolution factorises because the exponentials satisfy en(k+j)=en(k)en(j).
Aliasing is the same identity read the other way. The sum is 0 when m≡ℓ, but it is N whenever m≡ℓ — not only when m=ℓ. A finite sample has no way to distinguish the frequency n from the frequency n+N, so the discrete transform reports their sum.
1. Orthogonality. Put q=ζm−ℓ=e2πi(m−ℓ)/N, so (em,eℓ)=∑k=0N−1qk. If m≡ℓ modulo N then q=1 and the sum is N. Otherwise q=1, so the geometric sum equals (1−qN)/(1−q); and qN=e2πi(m−ℓ)=1 because m−ℓ is an integer, so the sum is 0. The N vectors e0,…,eN−1 are therefore orthogonal and non-zero, hence linearly independent, hence a basis of the N-dimensional space V; dividing by ∥eℓ∥=N makes it orthonormal.
2. Inversion. Expand F in that orthonormal basis: F=∑n(F,N−1/2en)N−1/2en=∑nN1(F,en)en. Now N1(F,en)=N1∑kF(k)ζnk=F^(n), which is the definition. So F=∑nF^(n)en, which is the claim.
3. Parseval. Writing cn=(F,N−1/2en)=NF^(n) and likewise dn for G, Parseval's identity in an orthonormal basis gives (F,G)=∑ncndn=N∑nF^(n)G^(n). Taking G=F gives the norm statement. The asymmetry — a factor 1/N on the sample side and none on the transform side — is forced by our normalisation; it is the discrete counterpart of the 1/L in the Fourier coefficient of Fourier series.
4. Shift. Substituting k′=k−j and using that both F and e−2πink/N are N-periodic in k, so that the sum over any N consecutive indices is the same,
where the exponential was split by en(k)=en(j)en(k−j). Take the inner sum first. Substituting k′=k−j and using, as in part 4, that both G and e−2πink′/N are N-periodic in k′, so that any N consecutive indices give the same total,
again because ∑jF(j)e−2πinj/N=NF^(n) by definition. So F∗G(n)=NF^(n)G^(n).
6. Aliasing. First, f really is the sum of its Fourier series at every point. Absolute convergence makes ∑mf^(m)e2πimx converge uniformly, by the Weierstrass M-test, to a continuous 1-periodic function u; integrating the uniformly convergent series term by term against e−2πimx gives u^(m)=f^(m) for every m. So f−u is continuous with all Fourier coefficients zero, hence ∥f−u∥2=0 by Parseval (Mean-square convergence and Parseval's identity), and a continuous function of zero mean square vanishes identically. This is where the continuity of f is used, and it is used nowhere else; without it the two sides could differ at the N sample points. Now insert the series at x=k/N:
The interchange of the two sums is legitimate because ∑m∣f^(m)∣<∞ — this is the one place absolute convergence is used — so the double family is absolutely summable and Absolute and conditional convergence permits the rearrangement. By part 1 the inner average is 1 when m≡n modulo N and 0 otherwise, leaving F^(n)=∑m≡nf^(m). If f^ vanishes outside ∣m∣<N/2 then for each n in that range the congruence class of n meets the support only at m=n, since two distinct members of a class differ by at least N while the support has width less than N. □
Let N=2M be even and let F∈V. Split the samples by parity,
F0(r)=F(2r),F1(r)=F(2r+1),r=0,…,M−1,
two vectors of length M, and write F^0,F^1 for their M-point transforms. Then, with ω=e−2πi/N,
F^(n)=21[F^0(n)+ωnF^1(n)],n=0,1,…,N−1,
in which F^0(n) and F^1(n) are read M-periodically, so that the two half-length transforms are each used twice.
Cost. Applying the splitting recursively when N=2n gives the fast Fourier transform. Stein and Shakarchi count one addition or one multiplication of complex numbers as one operation; their Lemma 1.4 is the recursion #(2M)≤2#(M)+8M for the minimum number #(⋅) of operations, and their Theorem 1.3 concludes
#(N)≤4Nlog2N=O(NlogN)
against the bound 2N2+N they give for evaluating the defining sums directly. Arfken, Weber and Harris count only multiplications and report (N/2)log2N for the Cooley–Tukey procedure against about N2 for the direct calculation; the two statements are consistent, counting different things, and we follow Stein and Shakarchi's accounting below. The ratio Arfken quotes belongs to his count and not to theirs: on multiplications alone, N=1024 gives (N/2)log2N=5120 against about N2=1.05×106, a saving of a factor 204.8 — "better than 200".
Proof(The splitting formula and the cost)AdaptedAdapted from Stein & Shakarchi (2003), Ch. 7, Lemma 1.4, p. 225
The formula. Separate the defining sum into even and odd indices, k=2r and k=2r+1:
Now ω2=e−2πi/M is the corresponding root of unity for length M, so ω2nr=e−2πinr/M. Pulling ωn out of the second sum and inserting a factor M/M in each,
The M-point transforms are M-periodic in n, which is why the M values of each suffice for all N=2M values of n: this is the entire saving, and it is where the evenness of N is used.
The cost. Suppose the two half-length transforms are available at a cost of #(M) operations each. Producing the powers ω,ω2,…,ω2M costs at most 2M operations, and gives ω2, the root of unity needed at length M, for free. Each of the 2M outputs then costs at most three further operations — two multiplications and one addition — so
#(2M)≤2M+2#(M)+3⋅2M=2#(M)+8M,
which is Stein and Shakarchi's Lemma 1.4. Induct on n where N=2n. At n=1 the two outputs are 21(F(0)±F(1)), costing at most five operations, and 4⋅2⋅1=8≥5. If #(2n−1)≤4⋅2n−1(n−1) then
#(2n)≤2⋅4⋅2n−1(n−1)+8⋅2n−1=4⋅2n(n−1)+4⋅2n=4⋅2nn,
which is 4Nlog2N. □
Remark(Where the hypotheses are used, and what goes wrong without them)Original
Periodicity of the index. Parts 4 and 5 need F and G extended N-periodically. This is not a convenience: the convolution the DFT diagonalises is the circular one, in which data runs off the right end and reappears on the left. If one wants the ordinary (linear) convolution of two finite sequences of lengths L1 and L2, both must first be padded with zeros to a common length N≥L1+L2−1; otherwise the tail wraps around and corrupts the head. This single point is the commonest bug in fast convolution code.
Absolute convergence in the aliasing formula. It is used only to exchange the two sums. Without it the formula can fail: the right-hand side need not converge. For a merely continuous f whose Fourier series is not absolutely convergent, the left-hand side is still perfectly well defined — it is a finite sum — but it need not be the sum of the coefficients in a congruence class.
The band-limit ∣m∣<N/2. The last sentence of part 6 needs the support of f^ to have width less than N, so that no two frequencies in it are congruent modulo N. At the boundary the conclusion fails: the frequencies m=N/2 and m=−N/2are congruent modulo N, so their coefficients are added and cannot be separated. The worked example below shows this happening, with a component that vanishes entirely. This boundary case is the Nyquist frequency, and it is why The sampling theorem asks for a strict band limit.
Even N in the FFT. The splitting needs N=2M. For a prime N no such splitting exists and the radix-2 algorithm has nothing to say; general N is handled by splitting on any factor, and a prime length needs a different idea.
No smoothness anywhere. Parts 1–3 are statements about N complex numbers and assume nothing at all. Parts 4 and 5 assume one thing only, and it is bookkeeping rather than analysis: that the index is read modulo N, as the first bullet above says. Only part 6, which relates the DFT to a function on the circle, puts a hypothesis on a function — continuity and an absolutely convergent Fourier series, the first used to identify f with the sum of its series and the second to exchange the sums.
Why the DFT is not the Fourier transform. The DFT of a sampled signal estimates the Fourier coefficients of its periodic extension. If the data does not in fact come from a periodic function, the implied jump between the last sample and the first is real as far as the transform is concerned, and its coefficients decay only like 1/n (Smoothness and decay of Fourier coefficients). This is spectral leakage, and it is a separate effect from aliasing: leakage comes from truncating in time, aliasing from sampling in time.
Suppose a signal f(t) is sampled N times at spacing τ, giving fk=f(kτ), k=0,…,N−1, over a total record length T=Nτ.
Fix the frequency grid. The N available frequencies are ωn=2πn/T, equivalently νn=n/(Nτ) cycles per unit time. Hassani derives this spacing by requiring the discrete orthogonality to hold, which is what forces ΔωΔt=2π/N. The resolution 1/T is set by the record length and the highest frequency by the sample spacing — the two are independent.
Fold the index to a frequency. The output index n runs 0,…,N−1, but the frequencies it represents run 0,1,… up to N/2 and then negative: for n>N/2, bin n carries the frequency (n−N)/(Nτ). For real data this makes F^(N−n)=F^(n), so only the first half of the output is independent. Plotting bins N/2 to N−1 as high frequencies is a standard and badly misleading error.
Identify the Nyquist frequencyνNyq=1/(2τ), half the sampling rate. Any content above it is reported at a false, lower frequency by part 6 above. If the signal may contain such content, filter it out before sampling; no processing afterwards can undo it.
Scale to the quantity you want. The F^(n) are estimates of Fourier coefficients (dimensionless fractions of the signal). To estimate the continuous Fourier transform instead, read the defining sum as a Riemann sum: ∫f(t)e−iωntdt≈NτF^(n)=TF^(n). This is Riley, Hobson and Bence's Exercise 13.10(a) up to their (2π)−1/2; they also note there that the estimate is good for n<N/2 and poor above it, which is step 2 restated.
Choose N a power of two if the transform is to be computed many times, so that the radix-2 recursion applies all the way down. Zero-padding to the next power of two is harmless for spectral estimation — it interpolates the frequency grid without adding information — but it is not harmless if a circular convolution of a specific length is wanted.
Check the answer.F^(0) must be the mean of the samples; for real data F^(N−n)=F^(n); and Parseval must hold, N1∑k∣fk∣2=∑n∣F^(n)∣2. These three cost nothing and catch almost every normalisation error.
Example(Eight samples, the radix-2 step, and an alias)Original
Take the 1-periodic signal
f(x)=1+2cos2πx−sin4πx,
whose Fourier coefficients we know exactly: writing the trigonometric functions as exponentials, f^(0)=1, f^(±1)=1, f^(2)=−2i1=2i, f^(−2)=−2i, and all other coefficients zero.
(a) Sampling at N=8. With xk=k/8 the samples fk=1+2cos(πk/4)−sin(πk/2) are
F=(3,2,1,2−2,−1,−2,1,2+2)
(to six decimals, 3,1.414214,1,0.585786,−1,−1.414214,1,3.414214). Their sum is 8, so F^(0)=1, already matching f^(0).
(b) The radix-2 step, done by hand. Split by parity: F0=(3,1,−1,1) from k=0,2,4,6 and F1=(2,2−2,−2,2+2) from k=1,3,5,7. Their four-point transforms, computed from the definition with e−2πi/4=−i, are
the second entry of F^1 because 41[2−i(2−2)+2+i(2+2)]=41[22+22i]. Now recombine with ω=e−iπ/4=22(1−i), remembering that F^0,F^1 are 4-periodic:
n
F^0(n)
F^1(n)
ωn
F^(n)=21[F^0+ωnF^1]
0
1
1
1
1
1
1
22(1+i)
22(1−i)
21[1+1]=1
2
0
−1
−i
21[i]=i/2
3
1
22(1−i)
22(−1−i)
21[1−1]=0
4
1
1
−1
0
5
1
22(1+i)
22(−1+i)
21[1−1]=0
6
0
−1
i
−i/2
7
1
22(1−i)
22(1+i)
21[1+1]=1
The products are all of the form 21(1±i)(1∓i)=21⋅2=1 or 21(1∓i)(−1∓i)=21⋅(−2)=−1, which is why the arithmetic closes so neatly.
Reading the table with the folding rule of step 2 — bins 5,6,7 are the frequencies −3,−2,−1 — gives F^(0)=1, F^(±1)=1, F^(2)=i/2, F^(−2)=−i/2 and nothing else: exactly the Fourier coefficients of f. This is part 6 in the good case, since f^ is supported in ∣m∣≤2<4=N/2. Parseval checks the whole computation in one line. On the sample side (2−2)2+(2+2)2=2(4+2)=12 and the remaining squares are 9+2+1+1+2+1=16, so
81k=0∑7fk2=812+16=828=3.5,
while on the transform side ∑n=07∣F^(n)∣2=1+1+1+41+41=3.5. The two agree.
(c) Sampling at N=4: an alias and a vanishing act. Now take xk=k/4. The samples are F=(3,1,−1,1), the vector F0 above, and its transform is F^=(1,1,0,1). Compare with the aliasing formula:
Every entry matches. The last line is the interesting one: the component −sin4πx, of amplitude 1, has disappeared completely. Its frequency 2 is exactly the Nyquist frequency N/2 for four samples, so the bins m=2 and m=−2 have collapsed onto one another and the two coefficients, being negatives, have cancelled. The signal 1+2cos2πx and the signal 1+2cos2πx−sin4πx have the same four samples, and no amount of analysis of those four numbers can separate them. A sine at exactly the Nyquist frequency is sampled at its zeros.
Application(What the fast transform is used for)Original
Spectral estimation. Every digital spectrum analyser, vibration monitor and audio equaliser computes a DFT of a windowed block of samples and plots ∣F^(n)∣ against νn=n/(Nτ). The two facts a user must carry are those of the method above: the resolution is 1/T and the ceiling is 1/(2τ). Doubling the record length sharpens the peaks; doubling the sampling rate raises the ceiling; neither does the other's job.
Fast convolution and filtering. Part 5 turns an O(N2) convolution into three transforms and N multiplications, so O(NlogN). This is how long digital filters are applied, how large integers and polynomials are multiplied in computer algebra systems, and how cross-correlations are computed in radio astronomy and in seismology. The zero-padding caveat of the remark is what separates a working implementation from a broken one.
Spectral methods for differential equations. On a periodic domain, differentiation is diagonal in the Fourier basis: f′(n)=2πinf^(n) for the series of Smoothness and decay of Fourier coefficients, and the same multiplication acts on the DFT of the samples. A pseudo-spectral solver therefore steps a PDE by transforming, multiplying by in, transforming back, and evaluating the nonlinear terms in physical space — each step O(NlogN) instead of the O(N2) of a dense differentiation matrix. The same idea solves Poisson's equation on a periodic box in one pass, by dividing by −4π2∣n∣2.
Crystallography and imaging. A diffraction pattern is the squared modulus of a Fourier transform, and the reconstruction of a density from measured intensities is done on a grid with the FFT. In magnetic resonance imaging the measured data are samples of the transform, and the image is a single inverse DFT; aliasing there appears as the familiar wrap-around artefact, in which anatomy outside the chosen field of view is drawn folded back inside it — part 6, in two dimensions, made visible.
Proves that a function whose Fourier transform vanishes outside [−a,a] is recovered exactly from its values at the points k/2a by a sum of shifted sinc functions, identifies 2a as the Nyquist rate below which distinct frequencies become indistinguishable, and shows what the same series returns when the function is not band-limited.
Definition(Band-limited functions, bandwidth and the sinc)AdaptedAdapted from Folland (1999), §8.3, Exercise 15, p. 255; Stakgold & Holst (2011), §2.4, p. 147
We use the Fourier transform in the convention that puts the 2π in the exponent,
f^(ξ)=∫−∞∞f(x)e−2πixξdx,
so that the inversion formula carries no constant at all. This is Folland's convention and the one in which the sampling theorem is cleanest, but it is local to this unit and is not the handbook's default: the Fourier transform defines f^(k)=∫f(x)e−ikxdx, whose inversion formula carries 2π1. The dictionary is ξ↦2πξ: the f^(ξ) of this unit is the f^(2πξ) of that one. Read through that dictionary, the inversion formula of Fourier inversion and Plancherel's identity of Plancherel's theorem lose their constants, which is what every step below uses. The translation to the angular-frequency convention ω=2πξ is given at the end of the statement.
The normalised sinc is
sincx=πxsinπx(x=0),sinc0=1,
Folland's definition. It is continuous, even, equal to 1 at x=0 and to 0 at every other integer, and it is the transform pair of an interval: χ[−a,a](x)=2asinc2ax. Warning on conventions: Stakgold and Holst, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §13.4, Exercise 13.4, p. 460), write sincx=sinx/x without the π, so their sampling formula reads sincΩ(x−xk) where ours reads sinc(2ax−k). Every statement below uses the normalised one.
For a>0 let
Ha={f∈L2(R):f^(ξ)=0 for almost every ∣ξ∣>a}.
A member of Ha is band-limited, and in the language of signal analysis, which Folland adopts, its bandwidth is 2a — the full width of the band, counting negative frequencies. The spacing h=1/(2a) is the Nyquist spacing: it is the reciprocal of the bandwidth, which is how Folland puts it on his p. 255. The corresponding rate 1/h=2a samples per unit of x is the Nyquist rate; numerically it equals the bandwidth 2a, and it is twice the highest frequency a present. Rate and spacing are reciprocal to one another, so only one of the two can be the reciprocal of the bandwidth, and it is the spacing.
Stakgold and Holst state the same condition in angular frequency: ϕ is Ω-band-limited when ϕ^(ω)=0 for ∣ω∣≥Ω, and Ω is called the bandwidth. Their Ω is our 2πa, and their sample spacing π/Ω is our 1/(2a).
Warning on the word bandwidth. The two books use it for two different numbers, and we keep both because we cite both. Folland's bandwidth is the two-sided width 2a in cycles per unit of x; Stakgold and Holst's bandwidth is the one-sided cut-off Ω=2πa in radians per unit of x. Their ratio is Ω/(2a)=π: one factor 2π from the change of frequency unit, one factor 21 from counting one side of the band instead of two. In this unit "bandwidth" unqualified always means Folland's 2a; wherever Stakgold and Holst's Ω is meant it is written Ω.
Theorem(Shannon sampling theorem)AdaptedAdapted from Folland (1999), §8.3, Exercise 15b–c, p. 255; Stakgold & Holst (2011), §2.4, Eqs. (2.4.10)–(2.4.10a), p. 148; Strauss (2008), §12.3, Exercise 6, p. 348
Let a>0.
1.Ha is a Hilbert space, and the functions
sk(x)=2asinc(2ax−k),k∈Z,
form an orthonormal basis of it.
2. (Sampling theorem.) Every f∈Ha agrees almost everywhere with a function in C0(R) — continuous and tending to 0 at ±∞ — and, taking that continuous version,
f(x)=k=−∞∑∞f(2ak)sinc(2ax−k)for every x∈R,
the series converging both uniformly on R and in L2(R). In words: a signal of bandwidth 2a is determined by its values on a lattice of spacing 1/2a, the reciprocal of the bandwidth.
3. (Energy.)k∈Z∑f(2ak)2=2a∥f∥22, so the samples are square-summable and carry exactly the energy of the signal.
4. (Aliasing.) Suppose f∈L1∩L2 does not lie in Ha — either because it is not band-limited at all, or because its band is wider than [−a,a]. Suppose moreover that f^∈L1(R) and that the periodised series ∑mf^(ξ+2am) converges absolutely and uniformly on [−a,a], and take for f the continuous representative
f(x)=∫−∞∞f^(ξ)e2πixξdξ,
which exists because f^∈L1 and which makes the individual samples f(k/2a) meaningful numbers rather than values of an equivalence class. Then the series of part 2, built from those samples, does not reproduce f. It reproduces the function fal with
f^al(ξ)=χ[−a,a](ξ)m∈Z∑f^(ξ+2am),
the spectrum of f folded into the band [−a,a] in strips of width 2a. Content at frequency ξ0 with ∣ξ0∣>a is reported at the frequency ξ0−2am that lies in the band: this is aliasing, and the discrete version of the same statement is part 6 of The discrete Fourier transform and the FFT.
In Stakgold and Holst's angular-frequency notation, part 2 is their Eq. (2.4.10a): an Ω-band-limited ϕ satisfies ϕ(x)=∑kϕ(xk)sinc[Ω(x−xk)] with xk=kπ/Ω and their unnormalised sinc. Strauss states the case Ω=π, in which the samples fall at the integers.
Idea of the proofOriginal
The hypothesis is that f^ lives on a bounded interval. A function on a bounded interval has a Fourier series, so f^ has one — and the Fourier coefficients of f^ on [−a,a] are, by the inversion formula read backwards, the values of f at the lattice points k/2a — up to the factor 1/(2a) and a reflection of the index, as Step 2 works out. So the samples are not an arbitrary list of numbers: they are a complete set of coefficients for f^, and therefore determine f^, and therefore determine f.
To get the interpolation formula, transform the series back. Each exponential e−πikξ/a, cut off to the band, inverts to a shifted sinc, because the transform of an interval is a sinc. That is the whole proof: a Fourier series in the frequency variable becomes a sinc series in the space variable.
The shape of the result is also forced by the first statement. The sincs sk are an orthonormal basis of Ha, so f=∑k(f,sk)sk for some coefficients; the content of the theorem is that the coefficient (f,sk) is nothing but the sample f(k/2a) divided by 2a.
Throughout, g=f^, which by hypothesis is supported in [−a,a] and lies in L2[−a,a].
Step 1: f is continuous and bounded. By the Cauchy–Schwarz inequality on an interval of finite length,
∫−aa∣g∣≤2a∥g∥2<∞,
so g∈L1(R) as well as L2. The inversion theorem of Fourier inversion is not what justifies the next line: every version of it assumes f∈L1, and a band-limited f need not be integrable — g=χ[−1/2,1/2] gives f=sinc, which is not. The L2 theory is in Plancherel's theorem instead, and it is enough. Its part 3 extends the transform to a unitary map of L2(R) onto itself, so the inverse map F−1 is defined on all of L2 and f=F−1g, by the definition g=f^. On the overlap L1∩L2 — where g lies, by the display above — that extension is still given by the defining integral rather than by a limit in the mean, as the definition block of the same unit records; and in the convention local to this unit the inverse transform is the transform followed by a reflection, (F−1g)(x)=(Fg)(−x), so the same is true of it. Hence, for almost every x,
f(x)=∫−aag(ξ)e2πixξdξ.
The right-hand side is continuous in x by dominated convergence and tends to 0 at ±∞ by the Riemann–Lebesgue lemma, which for an L1 function is the statement F(L1)⊂C0 of The Fourier transform; so it is the C0 representative promised in part 2. From now on fis that integral. Note for later that this formula is a bounded linear map: writing (Tg)(x) for the integral,
∣(Tg)(x)∣≤∫−aa∣g∣≤2a∥g∥L2[−a,a]for every x,
again by Cauchy–Schwarz. So convergence of functions in L2[−a,a] forces uniform convergence of their images under T. This is the only analytic point in the proof, and it is where the boundedness of the band is used.
Step 2: the Fourier coefficients of g are the samples of f. Regard g as a function on the interval [−a,a], of length 2a, and expand it in the Fourier series of Fourier series with period L=2a:
the series converging in L2[−a,a] because the normalised exponentials {(2a)−1/2eπikξ/a} are not merely orthonormal there but complete: that is the mean-square convergence of Mean-square convergence and Parseval's identity, and it is the only fact about Fourier series this proof needs. Compare the integral for ck with the inversion formula of Step 1 evaluated at x=−k/(2a): the exponents match, since 2πixξ=−πikξ/a there. Hence
ck=2a1f(−2ak).
The samples of f have appeared without being asked for; this is the heart of the matter.
valid at u=0 too, where both sides equal 2a. This is Folland's pair χ[−a,a]=2asinc2a⋅ read as an inverse transform. Applying T to the kth term of Step 2, with u=x+k/(2a),
Step 4: pass to the limit. The partial sums gK=∑∣k∣≤Kckeπikξ/a converge to g in L2[−a,a] by Step 2. By the bound in Step 1, TgK→Tg=funiformly on R. By Step 3,
relabelling k↦−k, which is legitimate because the sum is symmetric and finite. Letting K→∞ gives part 2 with uniform convergence. For the L2 convergence, Plancherel's theorem (Plancherel's theorem), in the convention local to this unit — where the transform is unitary, the (2π)−1/2 of that unit having been absorbed by the dictionary — gives ∥f−TgK∥L2(R)=∥g−gK∥L2[−a,a]→0, since the transform of TgK is gK extended by zero.
Step 5: the orthonormal basis. Applying T to the function 2a1e−πikξ/a and using Step 3's integral with u=x−k/(2a) shows that the transform of x↦sinc(2ax−k) is 2a1e−πikξ/aχ[−a,a](ξ). These lie in Ha, and by Plancherel
(sk,sl)=2a∫−aa4a21e−πi(k−l)ξ/adξ=δkl,
using Step 3's integral with u=(l−k)/(2a), which is 2asinc(l−k)=2aδkl. They are complete in Ha because their transforms are, up to the constant (2a)−1/2, the exponential basis of L2[−a,a], and the transform maps Ha onto L2[−a,a] unitarily. Since Ha is the preimage under the unitary transform of the closed subspace of L2 of functions vanishing off [−a,a], it is closed in L2(R), hence a Hilbert space.
Step 6: the energy identity. Expanding f in the orthonormal basis of Step 5 (Orthonormal bases in Hilbert space), f=∑k(f,sk)sk; comparing with part 2, (f,sk)=f(k/2a)/2a. Parseval's identity along that basis then reads ∥f∥22=∑k∣f(k/2a)∣2/(2a), which is part 3.
Step 7: aliasing. Drop the band-limitation and assume instead the hypotheses of part 4. Since f∈L1, f^ is continuous and bounded, so each term of P(ξ)=∑mf^(ξ+2am), the periodisation of f^ with period 2a, is continuous; by hypothesis the series converges uniformly on [−a,a], so P is continuous and bounded there and in particular lies in L2[−a,a]. Uniform convergence also lets the series be integrated term by term against the bounded factor e−πikξ/a, and the resulting strips reassemble into one integral over R, absolutely convergent because f^∈L1(R):
the exponential being unchanged by the shift ξ↦ξ+2am because e−πik(2am)/a=e−2πikm=1. The last integral is the inversion integral of part 4 at the point −k/(2a) — again f^∈L1 is what makes it an honest absolutely convergent integral and its value the sample of the continuous representative — so the coefficient equals 2a1f(−k/2a). (This chain is exactly the Poisson summation formula of Poisson summation formula applied to f^.) So P has on [−a,a] exactly the coefficients ck=2a1f(−k/2a) that Step 2 produced, and running Steps 3 and 4 with P in place of g shows that the sinc series converges to the function whose transform is Pχ[−a,a]. That is part 4; and when fis band-limited, P=f^ on the band and nothing changes, which is the second, Poisson-summation route to part 2 that Stakgold and Holst mention. □
Remark(Where the hypotheses are used)Original
Boundedness of the band. Used twice, both times in Step 1: to get f^∈L1, which makes f continuous so that pointwise samples mean anything at all, and to get the bound ∣Tg(x)∣≤2a∥g∥2 that upgrades L2 convergence of the Fourier series to uniform convergence of the sinc series. For f∈L2 alone, "the value of f at k/2a" is not even defined, since f is an equivalence class.
Sampling at or above the Nyquist rate. Step 2 expands f^ with period exactly 2a. Sampling more slowly, at spacing h>1/2a, amounts to using a period 1/h<2a, which does not cover the support of f^; the strips then overlap and part 4 takes over. Sampling faster is harmless: if f∈Ha then f∈Ha′ for every a′>a, and the theorem applies with the finer lattice.
The strictness at the band edge. Part 2 says nothing about a function whose spectrum has an atom exactly at ±a — in L2 a single point is a null set, so such a function is simply not distinguished. The reason engineers insist on sampling strictly faster than the Nyquist rate is visible in the classical counterexample f(x)=sin(2πax): every sample f(k/2a)=sinπk is zero, yet f=0. This is not a counterexample to the theorem as stated, because sin is not in L2(R); but it is why Stakgold and Holst define Ω-band-limited with the closed condition ϕ^(ω)=0 for ∣ω∣≥Ω, which excludes it, and why in practice one oversamples.
Convergence is absolute but slow. By part 3 the samples lie in ℓ2, and for each fixed x one also has ∑k∣sinc(2ax−k)∣2=1. This is part 3 applied to one particular signal, not Step 5, which proves only orthonormality: fix x and put u(y)=sinc(2a(y−x)), which by Step 3 with u-variable y−x equals 2a1∫−aae−2πixξe2πiyξdξ and so lies in Ha with transform 2a1e−2πixξχ[−a,a](ξ); Plancherel gives ∥u∥22=2a/(2a)2=1/(2a), so part 3 applied to u reads ∑k∣u(k/2a)∣2=2a∥u∥22=1; and u(k/2a)=sinc(k−2ax)=sinc(2ax−k), the sinc being even. Then Cauchy–Schwarz makes the series absolutely convergent at every point. But sinc decays only like 1/∣k∣, so truncating at ∣k∣≤K leaves an error controlled only by the ℓ2 tail of the samples. Reconstruction from a short record is correspondingly poor near the ends of the record, which is why practical interpolators use a windowed or tapered kernel rather than the bare sinc.
No band-limited signal has compact support. If f∈Ha then Step 1's formula makes sense for complex x: the integrand is continuous on the compact set [−a,a] and entire in x, so differentiation under the integral sign (Differentiation under the integral sign) shows f extends to an entire function. An entire function vanishing on an interval vanishes identically, by the identity theorem (Analytic continuation). So a non-zero signal cannot be both band-limited and time-limited, and every real measurement — which is time-limited — is strictly speaking undersampled. The theorem is an idealisation, and the practical question is always how much aliasing one is willing to tolerate.
What part 4 needs. Two things, and the same hypothesis buys both. First, f^∈L1(R): this makes the inversion integral f(x)=∫f^(ξ)e2πixξdξ absolutely convergent, hence gives a continuous representative of f whose individual samples f(k/2a) are defined — without it the statement has no subject, by the first bullet above — and it is what makes the unfolded integral in Step 7 legitimate. Second, absolute and uniform convergence of the periodised series on [−a,a]: this is what licenses the term-by-term integration there, and it puts P in L2[−a,a] so that Steps 3 and 4 may be rerun. Decay of f^ faster than ∣ξ∣−1 — say ∣f^(ξ)∣≤C(1+∣ξ∣)−1−ε — gives both at once, by comparison with ∑m∣m∣−1−ε; a compactly supported f^, as in the example, gives them trivially, the periodised series then having only finitely many non-zero terms on the band. Without any decay the folded series can diverge and nothing can be said.
Method(Sampling and reconstructing a signal)AdaptedAdapted from Stakgold & Holst (2011), §2.4, pp. 147–148
Find the highest frequencyνmax present, in cycles per unit of x. If the signal is not band-limited — and no real signal is — choose the frequency above which its content is negligible for the purpose at hand.
Band-limit it physically, before sampling. Pass the signal through a low-pass (anti-alias) filter with cut-off at νmax. This step cannot be postponed: once the samples are taken, the out-of-band content has already been folded into the band and no processing can separate it from the genuine content there, by part 4.
Sample at a rate at or above 2νmax. The Nyquist rate 2νmax is the theoretical minimum, and it is attained rather than merely approached — part 2 holds at the rate 2a itself; because no filter cuts off abruptly, practice uses 10–20% more. The sample spacing is h=1/(2νmax) at the Nyquist rate.
Reconstruct with f(x)=∑kf(kh)sinc((x−kh)/h), which is part 2 written with h=1/2a. In software this is done by zero-padding the discrete transform and transforming back, which is the same sum computed with the FFT of The discrete Fourier transform and the FFT; in hardware it is done by holding each sample and low-pass filtering, which approximates the sinc kernel.
Check the energy. The identity of part 3, h∑k∣f(kh)∣2=∥f∥22, is a cheap test that the sampling rate and scaling are consistent.
If the rate is fixed and too low, predict the alias rather than be surprised by it: content at ν is reported at the ν′ with ν−ν′ an integer multiple of 1/h and ∣ν′∣<1/(2h), and that ν′ is unique. The strict inequality matters: if ν is an odd multiple of the Nyquist frequency 1/(2h) then both ν′=+1/(2h) and ν′=−1/(2h) satisfy the congruence, the two cannot be told apart by the samples, and the pair of frequencies is reported as one — the degenerate band-edge case of part 4, which part 6 of The discrete Fourier transform and the FFT shows in its discrete form. Away from that edge the fold is unambiguous: a sinusoid just above the Nyquist frequency appears just below it, running backwards — the wagon-wheel effect.
For quadrature, use the trapezoidal rule. On a band-limited function the trapezoidal rule is not an approximation but an identity: ∫f=h∑kf(kh). Poisson summation (Poisson summation formula) gives h∑kf(kh)=∑mf^(m/h), and every term with m=0 vanishes as soon as 1/h>a, so the identity holds for every h<1/a — twice as coarse as the Nyquist spacing 1/(2a), because quadrature needs only the value f^(0) and not the whole band. Stakgold and Holst state it in their notation as 0<λ≤2π/Ω, that is h≤1/a with the non-strict inequality their closed band condition ϕ^(ω)=0 for ∣ω∣≥Ω allows. The example below verifies the identity at the Nyquist spacing 1/(2a), twice as fine as the quadrature bound requires.
Example(A root-raised-cosine pulse, and a carrier that lies about its frequency)Original
(a) Exact reconstruction. Let f^(ξ)=cos(πξ/2) for ∣ξ∣≤1 and 0 outside, so f∈H1: the bandwidth is 2, the Nyquist spacing is h=21. (This spectrum is the square root of the raised cosine 21[1+cosπξ], whence the name used in the application; the root and the full raised cosine behave quite differently at the sampling instants, as noted there.) The pulse itself is
since sin(2π±θ)=cosθ and A+B1+A−B1=A2−B22A with A=π/2, B=2πx. (At x=±41 both numerator and denominator vanish; l'Hopital's rule gives the removable value f(±41)=1, which also follows directly from ∫−11cos2(πξ/2)dξ=1.)
The samples. At xk=k/2,
f(2k)=π(1−4k2)4cosπk=π(1−4k2)4(−1)k,
that is 1.273240 at k=0, then 0.424413, −0.084883, 0.036378, −0.020210 for k=1,2,3,4 and the mirror images at negative k. They decay like k−2, so the series converges respectably.
Splitting 4k2−11=21[2k−11−2k+11] and using the two Leibniz sums 1−31+51−⋯=4π and 31−51+71−⋯=1−4π, the bracketed series is 21[4π−(1−4π)]=4π−21, so the whole expression is π2[1+2π−1]=1. The trapezoidal rule with spacing 21 gives the integral exactly, as it must.
(b) An alias. A pure cosine is not in L2(R) and so is not in any Ha; to put a carrier inside the theorem we give it a band-limited envelope. Rescale the pulse of part (a): q(x)=21f(x/2) has q^(ξ)=f^(2ξ)=cosπξ on ∣ξ∣≤21 and 0 outside, that is
q(x)=π(1−4x2)2cosπx,q(0)=π2,q(±21)=21,
the value at ±21 being the removable one, equal to ∫−1/21/2cos2πξdξ=21. Now modulate it onto a carrier of frequency 7 and onto a carrier of frequency 1:
g(x)=q(x)cos(2π⋅7x),g~(x)=q(x)cos(2π⋅1x).
Multiplying by cos2πνx splits the spectrum into two copies at ±ν: g^(ξ)=21[q^(ξ−7)+q^(ξ+7)] is supported in [6.5,7.5]∪[−7.5,−6.5] and g~^ in [0.5,1.5]∪[−1.5,−0.5]. Both q and the two modulated signals decay like x−2, so all three are in L1∩L2, and
Sample both at 6 per unit, h=61 — far too slow for g, comfortably fast for g~. Step 6 predicts the alias: 1/h=6 and 1/(2h)=3, and 7−6=1 with ∣1∣<3, so the carrier at 7 is reported at 1. The samples confirm it exactly, since cos(2π⋅7⋅6k)=cos(3πk+2πk)=cos(2π⋅1⋅6k):
(the value −0.5 at k=3 is q(21)cosπ), with the mirror images g(−k/6)=g(k/6) at negative k, both functions being even.
What the reconstruction returns. Apply part 2 to g~ with a=3, which is legitimate because g~∈H1.5⊂H3 and 1/(2a)=61=h. It gives g~(x)=∑kg~(k/6)sinc(6x−k). Since the samples of g are the same numbers, the sinc series built from g's samples converges to g~, not to g: the reconstruction has silently moved the carrier from 7 to 1. (Part 4 says the same thing in the frequency domain, and its hypotheses hold here: g∈L1∩L2 and g∈/H3; g^ is continuous — q^ vanishes at ξ=±21, so the two copies join the axis without a jump — and, being continuous with support in [6.5,7.5]∪[−7.5,−6.5], it lies in L1(R); and on [−3,3] the periodised series ∑mg^(ξ+6m) has only the two terms m=1 and m=−1 non-zero, so it converges absolutely and uniformly there. The periodisation of g^ with period 2a=6, cut to [−3,3], keeps only the copy at 7 shifted down by 6 and the copy at −7 shifted up by 6, giving 21[q^(ξ−1)+q^(ξ+1)]=g~^.)
The damage, at one point. Take x=121, halfway between the samples at k=0 and k=1. There q(121)=35π72cos(π/12)=35π18(6+2)=0.6324968, so
the same magnitude with the opposite sign — the reconstruction is not merely inaccurate, it is in antiphase with the truth. Summing the sinc series symmetrically over ∣k∣≤K confirms which of the two it finds:
K
4
10
50
200
partial sum at x=121
0.5471683
0.5477579
0.5477583
0.5477583
against g~(121)=0.5477583 and g(121)=−0.5477583. To recover gby the sinc series of part 2 one must sample at 15 per unit or faster: g∈H7.5, so part 2 applies with a=7.5 and the Nyquist rate 15 itself suffices. Fifteen is not, however, the slowest rate at which these samples contain g at all. The calculation just made shows that at 6 per unit the two copies of q^ fold into [−3,3] without overlapping, so a listener who knows that the signal lives in [6.5,7.5]∪[−7.5,−6.5] can shift them back and recover g exactly; that is bandpass sampling, and it needs a reconstruction built for the known band rather than the lowpass series of part 2, which has no way to tell g from g~. The alternative, when the band is not known in advance, is to filter the content near frequency 7 away before sampling. In practice one samples faster still, by the margin of step 3.
Application(Where the theorem is used)Original
Digital audio. Human hearing stops near 20 kHz, so the compact-disc standard samples at 44.1 kHz — above the Nyquist rate 40 kHz, with the margin of step 3 left for the anti-alias filter to roll off. The sincs of part 2 are the reconstruction filter; the theorem is the reason a sequence of numbers can be an exact representation of a continuous waveform rather than an approximation to one, and the reason the question "but what happens between the samples?" has the answer "nothing that is not already determined".
Instrumentation and the anti-alias filter. Any digitiser — an oscilloscope, a seismometer, a satellite radiometer, a climate data logger — carries an analogue filter ahead of its converter, and part 4 is its justification. Without it, a 60 Hz mains hum sampled at 100 Hz would appear as a 40 Hz oscillation indistinguishable from signal. The same arithmetic explains the wagon-wheel effect in film, where 24 frames per second alias the rotation of a wheel down to a slow forward or backward drift.
Magnetic resonance imaging and radio astronomy. Here the measured data are samples of the Fourier transform rather than of the signal, so the theorem applies in the transposed form: sampling f^ on a lattice of spacing Δ determines f exactly provided f vanishes outside an interval of length 1/Δ. The field of view is the reciprocal of the sample spacing in k-space, and anatomy outside it is folded back in — the wrap-around artefact, which is part 4 seen in an image.
Numerical quadrature and spectral accuracy. Step 7 and the example show that the trapezoidal rule is exact on band-limited functions. For a function that is merely very smooth and rapidly decaying, the spectrum is not compactly supported but is extremely small outside a band, and the trapezoidal error is correspondingly tiny — which is why the trapezoidal rule, usually a second-order method, converges faster than any power of h for smooth decaying integrands on the whole line, and is the method of choice there.
Information theory. Part 3 says that a signal occupying a band of full width W=2a, observed over a time T, is carried by about WT numbers: the sampling rate is 2a=W samples per unit time, so T units of time hold WT of them. The familiar "2WT" of the engineering literature is the same count with W read as the one-sided bandwidth — the highest frequency a — rather than as Folland's full width 2a. Under the convention of this unit the count is WT, equivalently 2aT; the warning in the definition about the two meanings of "bandwidth" is exactly what the factor 2 turns on. This count is the starting point of Shannon's channel-capacity theorem, and it is why bandwidth, not merely power, is the scarce resource in a communication link. The pulse of the example is not an accident either, though it must be named carefully. Its spectrum cos(πξ/2) on [−1,1] is the root-raised cosine with roll-off 1, and that is the standard transmit shape in digital communication: the receiver applies a matched filter with the same response, the two multiply in the frequency domain, and the cascade realises the full raised cosine cos2(πξ/2)=21[1+cosπξ]. It is the full raised cosine, not the transmit pulse, that has the zero-intersymbol-interference property: its inverse transform equals 1 at x=0 and vanishes at every other integer, so symbols spaced one unit of time apart are read off cleanly. The transmit pulse of the example does not vanish at those instants — its samples f(k/2)=4(−1)k/[π(1−4k2)] are all non-zero, which is exactly what makes the reconstruction in the example work — and its band of full width 2 is twice the width 1 that a symbol rate of 1 needs at minimum. That factor is the 100% excess bandwidth bought by the roll-off, and it is paid for the pulse's rapid decay; the trade between excess bandwidth and decay rate is the design question the raised-cosine family answers.