Poisson summation formula (${\sum }_{n\in \mathbb{Z}}f(n)={\sum }_{k\in \mathbb{Z}}\hat{f}(k)$)

The Poisson summation formula equates sample sums of a function and its Fourier transform. Exchanging broad and narrow sums can save work. Use the conventionFor every Schwartz function $f\in \mathcal{S}(\mathbb{R})$,Both series have absolute convergence. A sufficient hypothesis is that $f$ is smooth and ${\sup }_x(1+|x|{)}^M|f^{(r)}(x)|<\infty$ for all nonnegative integers $M,r$; its Fourier transform is then also a Schwartz function.^[ref] This assumption controls both decay and regularity; mere membership of an L1 space does not control values at individual sampling points, as periodization of an integrable function illustrates.

Here $\xi$ counts cycles per unit length. In the angular convention used elsewhere in Fourier analysis, $\tilde{f}(\omega )=\int f(x)e^{-i\omega x}dx$, so $\hat{f}(\xi )=\tilde{f}(2\pi \xi )$ and the right-hand side becomes ${\sum }_k\tilde{f}(2\pi k)$. Keeping the Fourier transform convention explicit prevents stray $2\pi$ factors. Changing the sample spacing leads to the Poisson summation formula for a Euclidean lattice.

Take the Gaussian function $f_a(x)=e^{-\pi a{x}^2}$ with $a>0$. Every derivative is a polynomial times the same Gaussian function, so $f_a$ is a Schwartz function. The Fourier transform of a Gaussian gives, in the cycles convention,Thus the Poisson summation formula gives the exact identityThe Gaussian Fourier transform is given in scaled form in Exercise 41 of Terence Tao's 245C, Notes 2: The Fourier transform.^[ref] The identity is also the imaginary-axis transformation of the Jacobi theta function. At $a=1$ the Gaussian function is its own Fourier transform; at reciprocal values of $a$ its spatial and spectral widths exchange.

For a short derivation of the Gaussian Fourier transform, write $H(\xi )={\hat{f}}_a(\xi )$. Differentiation under the integral and integration by parts, using $f_a^{\prime }(x)=-2\pi ax{f}_a(x)$, giveThe Gaussian integral gives $H(0)=a^{-1/2}$, so this ordinary differential equation yields the claimed formula. The decay of the Gaussian function justifies these steps and removes the boundary term.^[ref]

For $a=0.04$, the physical Gaussian function has width parameter $1/\sqrt{a}=5$, while its Fourier transform has width parameter $\sqrt{a}=0.2$. Many physical samples matter, but the dual zero sample already givesFor $a=25$, use the physical side instead:These are bounds on the full positive tails, obtained below, rather than estimates based only on the first omitted term. Ordinary floating-point arithmetic rounds these sums to $5$ and $1$; the exact real numbers are slightly larger. A Gaussian sum truncation bound lets us choose the cheaper side and certify its truncation error before evaluating it.

For $b>0$ and a nonnegative integer $N$, define the symmetric partial sum and its positive remainder bySuccessive omitted terms have ratio $e^{-\pi b(2n+1)}\le e^{-\pi b(2N+3)}$ for $n\ge N+1$. Bounding them by a geometric series givesConsequently the physical approximation to $S(a)$ has error at most $B_N(a)$, while the dual approximation $a^{-1/2}{S}_N(1/a)$ has error at most $a^{-1/2}{B}_N(1/a)$. The prefactor must be included when setting an absolute error tolerance.

For a truncation error below $10^{-12}$, at $a=0.04$ the physical side first passes this bound with $N=15$: $31$ lattice points give $B_{15}(0.04)<2.172\times 10^{-14}$, whereas $N=14$ gives $B_{14}(0.04)>10^{-12}$. The dual side needs only $k=0$. At $a=25$ the physical side needs only $n=0$, whereas the dual side first passes the bound with $N=14$: $29$ lattice points give $25^{-1/2}{B}_{14}(0.04)<2.146\times 10^{-13}$. These counts measure symmetric lattice terms; using the even function symmetry reduces the distinct exponential function evaluations to $N+1$. Truncation bounds do not include rounding error in a numerical implementation.

For $F\in L^1(\mathbb{R})$ the periodization is absolutely convergent almost everywhere in a period by integrating the sum of absolute values. Its Fourier series coefficient is $c_n=\tilde{F}(n)/(2\pi )$, where $\tilde{F}(\omega )={\int }_{\mathbb{R}}F(x)e^{-i\omega x}dx$. Absolute integrability permits integrating term by term and changing variables across the disjoint period intervals. This coefficient identity does not by itself assert pointwise convergence of the ordinary Fourier series. Even uniform convergence of the original periodization cannot repair arbitrary point values: $F=1_{\{0\}}$ has zero transform and a nonzero periodization at multiples of $2\pi$. Additional regularity, such as a continuous piecewise continuously differentiable periodization, validates the pointwise Poisson summation formula in this setting.

// Target: fluid-mechanics.bigb

For a Schwartz function on the real numbers, its periodization $P_f$ is a smooth periodic function of period of a function one. Rapid decay gives uniform convergence on $[0,1]$ of the defining series and each differentiated series. Its $k$th Fourier coefficient isThe change of variable $u=x+n$ uses $e^{2\pi ikn}=1$. Rapidly decreasing Fourier coefficients give the absolutely and uniformly convergent Fourier series $P_f(x)={\sum }_k\hat{f}(k)e^{2\pi ikx}$. Evaluation at $x=0$ proves the Poisson summation formula.^[ref]

The same identity explains aliasing in the Nyquist–Shannon sampling theorem: applying the Poisson summation formula to the Schwartz function $x\mapsto f(x)e^{-2\pi i\xi x}$ givesThe sampled spectrum is periodic in frequency, and shifted copies of the continuous Fourier transform add together. When those copies overlap, the samples cannot separate their contributions. The Nyquist–Shannon sampling theorem adds a bandlimited function hypothesis to make reconstruction possible.

Let $\Lambda =B{\mathbb{Z}}^d$ be a full-rank Euclidean lattice, with covolume $m(\Lambda )=|\det B|$ and dual lattice ${\Lambda }^{\vee }=B^{-\mathsf{T}}{\mathbb{Z}}^d$. For a Schwartz function $f\in \mathcal{S}({\mathbb{R}}^d)$ use the Fourier transform conventionThen the Poisson summation formula becomesBoth series have absolute convergence. To see the normalization, apply integer-lattice Poisson summation to $g(u)=f(Bu)$. The change of variables formula gives $\hat{g}(k)=|\det B{|}^{-1}\hat{f}(B^{-\mathsf{T}}k)$, yielding both the reciprocal covolume factor and the dual lattice argument. The factor is the physical lattice-point density. In particular, for spacing $h>0$,For the Gaussian function $f(x)=e^{-\pi a\Vert x{\Vert }^2}$ with $a>0$, the product of the one-dimensional Gaussian Fourier transforms givesHere $\Vert \cdot \Vert$ is the Euclidean norm. As a normalization check, take $\Lambda =\mathrm{diag}(2,3){\mathbb{Z}}^2$. Its covolume is $6$, and its dual lattice is $\mathrm{diag}(1/2,1/3){\mathbb{Z}}^2$, soThis also follows by multiplying two Gaussian Poisson summation identities with parameters $4a$ and $9a$. If one instead uses an angular Fourier transform kernel $e^{-i⟨x,\omega ⟩}$, the reciprocal frequencies are $2\pi {\Lambda }^{\vee }$, while the prefactor remains $m(\Lambda {)}^{-1}$.