Truth of the Poisson summation formula The Poisson summation says, roughly, that summing a smooth $L^1$-function of a real variable at integral points is the same as summing its Fourier transform at integral points(after suitable normalization). Here is the wikipedia link.
For many years I have wondered why this formula is true. I have seen more than one proof, I saw the overall outline, and I am sure I could understand each step if I go through them carefully. But still it wouldn't tell me anything about why on earth such a thing should be true.
But this formula is exceedingly important in analytic number theory. For instance, in the book of Iwaniec and Kowalski, it is praised to high heavens. So I wonder what is the rationale of why such a result should be true.
 A: In what follows, I'll use the convention $$ \hat{f}(\xi) = \int_{-\infty}^{\infty} f(x)e^{-2\pi i x \xi}dx,$$ so that $$ f(x) = \int_{-\infty}^{\infty} \hat{f}(\xi)e^{2\pi i x \xi}d\xi.$$
I like the following interpretation of Poisson summation, which also gives a generalization: Consider the Dirac comb distribution $C(x) = \sum_{n\in \mathbb{Z}} \delta(x-n)$. This is a tempered distribution, so it has a Fourier transform. In fact, it is its own Fourier transform. To justify this, I'm going to give a very nonrigorous argument. But if intuition is the main goal, then I think it will help. First, note that $C(x)$ is periodic with period 1. Thus, its "Fourier transform" is actually a Fourier series: its support is in $\mathbb{Z}$. This follows by noting that
$$\begin{align}C(x) &= \int_{-\infty}^{\infty} \hat{C}(\xi)e^{2\pi i x \xi}d\xi; \\
 C(x) &= \sum_{n\in \mathbb{Z}}a_n e^{2\pi i n x};
\end{align}$$
Where the first line is the Fourier inversion formula and the second line is the Fourier series for $C$. It follows by uniqueness that $\hat{C}(\xi) = \sum_{n\in \mathbb{Z}}a_n \delta(\xi - n)$. On the other hand, the (inverse) Fourier transform of $\hat{C}$ is also supported on $\mathbb{Z}$, so $\hat{C}$ is also periodic with period 1. Thus, all the $a_n$ are the same:
$$\hat{C}(\xi) = a\sum_{n\in \mathbb{Z}}\delta(\xi - n),$$ where $a$ is some scalar. It's not hard to see that the scalar has to be 1. 
To derive Poisson summation from this, use the convolution theorem: let $f$ be any function. On the one hand, $$(f*C)(x) = \sum_{n\in \mathbb{Z}} f(x+n).$$ On the other hand, we can use the convolution theorem: $$\widehat{(f*C)}(\xi) = \hat{f}(\xi)\hat{C}(\xi) = \hat{f}(\xi)\sum_{n\in \mathbb{Z}} \delta(\xi-n) = \sum_{n\in \mathbb{Z}} \hat{f}(n)\delta(\xi-n).$$ The last sum gives the Fourier series of the periodic function $f*C$: $$(f*C)(x) = \sum_{n\in \mathbb{Z}} \hat{f}(n)e^{2\pi i n x}.$$ Plugging in $x=0$ gives the Poisson summation formula, QED. But the result for general $x$ is interesting as well: given a function $f$, you can obtain a periodic function $g(x)$ by (a) adding up $f(x+n)$ over all integers $n$, or (b) taking the Fourier transform of $f$ at integer frequencies and making that the Fourier series of $g$. The result is that (a) and (b) give the same function. 
A: The Poisson formula is just the decomposition in Fourier series of the Dirac measure on the circle. 
It happens that the theory of distributions on the circle is far more simpler than the one on the real line since everything is compactly supported. In his book, Laurent Schwartz starts the chapter about the Fourier transform with the circle and shows that any distribution has a convergent Fourier series (VII.1.3 "theorie des distributions").
Moreover the space of distributions is the union of all the negative Sobolev spaces and these spaces can be caracterized in terms of the decay of the Fourier coefficients. I think that the Dirac is in $H^{-1}$. Anyway the usual formula holds and we have
$$
\delta_0(x) = \sum_n \langle \delta_0, e^{inx} \rangle e^{inx} = {1\over 2\pi}\sum_n e^{inx}, \quad x\in {\bf R}/{\bf Z}.
$$
Note that with this choice of normalization, the scalar product is ${1\over 2\pi}\int_0^{2\pi} f(x)\overline{g(x)} \,dx$ so that the $e^{inx}$ are unit vectors.This gives the ${1\over 2\pi}$ on the right hand side.
Now we compose this with the projection $\pi : {\bf R} \rightarrow {\bf R} /2\pi{\bf Z}$ to get the Poisson formula on ${\bf R}$. One must be a bit cautious because although we can always compose a function with another as soon as the domains match, it is not true that composition is always possible for distributions. 
In our particular case, the Dirac measure is the limit of an approximate identity, and so it is easy to guess what is happening just by composing the approximate identity and passing to the limit. Without surprise, the pullback of the Dirac measure at $\{0\}$ is given by the Dirac measure on the pullback of $\{0\}$, which is just the sum of all the Dirac measures on the points $2\pi{\bf Z}$.
$$
\sum_{n\in {\bf Z}} \delta_{2\pi n}(x) = {1\over 2\pi} \sum_{n\in {\bf Z}} e^{inx}, \quad x\in {\bf R}.
$$
We can evaluate this against a function in Schwartz space or use a convolution to obtain the original formula of Poisson. And we got the $2\pi$ at the correct locations without even thinking about it. So that's the truth.
A: It is a special case of the trace formula. Both sides are the trace of the same operator.
A: A variant on @Darsh Ranjan's argument: since $u=\sum_{n\in\mathbb Z}\delta(x-n)$ is invariant under translation by $\mathbb Z$, and is annihilated by multiplication by $e^{2\pi inx}-1$ for all $n\in \mathbb Z$ (this captures the order-zero aspect!) it suffices to show that there is a unique (up to scalar multiples) distribution meeting these conditions, and that the Fourier transform of $u$ also does meet them. The latter is immediate, since Fourier transform interchanges the two conditions.
For uniqueness: annihilation by multiplication by $e^{2\pi ix}-1$ implies that the support is $\mathbb Z$. We know the classification of distributions supported at points, hence, supported on a discrete closed subset: Dirac deltas and derivatives. Using smooth cut-offs to examine the behavior at a given integer, since $e^{2\pi ix}-1$ vanishes just to order $1$, the order of the $n$-th component is $0$. Thus, such a distribution is of the form $\sum_n c_n\cdot \delta(x-n)$. Translation invariance implies that all the coefficients are the same.
