I've seen the standard proof that the Fourier transform is self-inverse (up to an overall factor determined by conventions), which is essentially equivalent to $\int_{-\infty}^\infty e^{2 \pi i \omega x} dx = \delta(\omega)$. That isn't too hard to see from definitions.
But intuitively, the Fourier transform of a function gives the amplitude of each component frequency in an obvious sense. Is it intuitively obvious that this operation is self-inverse or is there a deep reason why it must be so?
This won't be the most rigorous answer, but just a representation my intuition with regards to the Fourier transform.
In nonrelativistic quantum mechanics there are two important operators which act on $L^2(\mathbb{R})$.
The position operator,
$$ (\hat{x} \psi)(x) = x\psi(x),$$
and the momentum operator,
$$ (\hat{p} \psi)(x) = -i \frac{\partial \psi}{\partial x} (x).$$
Now the "eigenvectors" of $\hat{x}$ are $a_{x'}(x) = \delta(x-x')$ and the "eigenvectors" of $\hat{p}$ are $b_p(x) = e^{ipx}$.
Now if you study linear algebra at all you know that the choice of basis is arbitrary. Also you would know that the eigenvectors of a Hermitian operator always form a complete basis. Although our operators our Hermitian, their eigen-functions certainly don't live in $L^2(\mathbb{R})$. I'm going to completely ignore that fact and pretend that I can use them as basis.
For instance the expansion of $\psi(x)$ in terms of $a_x$'s will be,
$$ \psi(x) = \int \mathrm dx'\ c(x') a_{x'}(x),$$
we can use what we know about delta functions to conclude that the appropiate expansion coefficients, $c(x')$, are in fact $c(x')=\psi(x')$.
If we want to expand $\psi$ in the momentum basis (the $b_p$'s) then we write the superposition and try to determine the coefficients.
$$\boxed{ \psi(x) = \int \mathrm dp \ c(p) b_{p}(x)},$$
the coefficient function, $c(p)$, is a vector in $L^2(\mathbb{R})$ as well. In fact what we have writtend down here is an isomorphism between the x-basis and the p-basis. However the bahavior of the operators on functions expressed in this new bases is different.
$$ \hat{x} \psi(x) = \int \mathrm dp \ c(p) x e^{ipx} = \int \mathrm dp \ c(p) \left(-i\frac{\partial}{\partial p} \right)e^{ipx} = \int \mathrm dp \ i \frac{\partial c(p)}{\partial p} e^{ipx} $$
$$ \boxed{ \hat{x} = i \frac{\partial}{\partial p}} $$
$$ \hat{p} \psi(x) = \int \mathrm dp \ c(p) \left(-i \frac{\partial}{\partial x} \right) e^{ipx} = \int \mathrm dp \ c(p) p e^{ipx} = \int \mathrm dp \ i \frac{\partial c(p)}{\partial p} e^{ipx} $$
$$ \boxed{ \hat{p} = p} $$
This means that in the p-basis the eigenvalues of $\hat{x}$ are $a_x(p) = e^{-ipx}$. I can write $c(p)$ in terms of the x-basis then by the integral,
$$ c(p) = \int \mathrm dx \ \phi(x) e^{-ipx} $$
I've used $\phi(x)$ to represent the coefficients in this expansion, but in practice we actually know the vector in the x-basis which is isomorphic to $c$, that would be $\psi(x)$.
$$ \boxed{ c(p) = \int \mathrm dx \ \psi(x) e^{-ipx} }$$
One issue is the normalization of the transforms are off. This is because the don't have finite norms. For that you need to actually prove the Fourier inversion theorem.
Consider the time-frequency $(t, f)$ plane interpretation. Suppose we start with $f(t)$ along time axis. Define the rotation matrix (called Elliptic group in this context) such that
$$ M_\gamma = \begin{pmatrix}\cos(\gamma) & -\sin(\gamma)\\ \sin(\gamma) & \cos(\gamma)\end{pmatrix}, \qquad \qquad \gamma = m \pi/2 $$
Thus, applying Fourier once rotates it toward frequency axes, apply it again, moves it in the opposite side of time axes (As you can verify through integration).
In general, to apply the Fourier transform $m$ times is to rotate the signal with $\gamma = \frac{m\pi}{2}$ degree (with time and frequency axes being orthogonal).
$$\mathcal{F}^m = \begin{cases}f(t)& \mbox{ if } m \equiv 0 \mod 4 \\\mathcal{F} & \mbox{ if } m \equiv 1 \mod 4\\ f(-t) & \mbox{ if } m \equiv 2 \mod 4\\ \mathcal{F}^{-1} & \mbox{ if } m \equiv 3 \mod 4\end{cases} $$
Hence, it is not self-inverse, but rather enjoys a 4-fold periodicity property.
It should also be noted that rotation by arbitrary degrees, (i.e. $m$ not integer) is at the heart of the Fractional Fourier transform.
EDIT ($M_\gamma$ and integral transform)
Briefly, multiplication of $f(t)$ by $t$ in time domain (let that be through an operator, $Q$) corresponds to differentiation (with rotation by -i) in frequency domain (through an operator $P$). Define the (linear) operator $C_M$ such that $C_M Q = CQC^-1 = dQ - bP$ and $CPC^{-1} = -cQ + aP$ parametrized by
$$M = \begin{bmatrix}a & b\\ c& d\end{bmatrix} \in SL(2, R)$$
When $M$ is a rotation matrix, i.e. $M = M_\gamma$, and specifically $\gamma = \frac{\pi}{2}$ (so that $a, d = \cos(\pi/2) = 0, b = -c = 1$), we get the Fourier transform (applied once). For example, its action on $tf(t)$ is $C_M tf(t) = CQC^{-1}F(f) = -P F(f) = i\frac{d F}{df}$.
More importantly, composition of $C_{M_1}$ and $C_{M_2}$ corresponds to multiplication of $M_1$ and $M_2$, which explains the four-periodicity. The argument for $Q$ (with similar one for $P$) is as follows:
\begin{align}C_{M_2}Q'C_{M_2}^{-1} &= C_{M_2}C_{M_1}Q(C_{M_2}C_{M_1})^{-1} \\&= C_{M_2}(d_1Q - b_1P)C_{M_2}^{-1} \\&= c_2b_1 + d_2d_1Q - (a_2b_1 + b_2d_1)P \\&= C_{M_1M_2}QC_{M_1M_2}^{-1}\end{align}
In terms of kernel of the integral transform, $\int df F(f) C_M(t,f)$ is defined to be $C_M(t, f) = \theta_M e^{i(at^2 - 2tf + d f^2)/2b}$, which is $e^{-itf}$ for Fourier transform case, up to a constant $\theta_M$.
All the interesting properties of integral transforms are encoded through the matrix $M$, see chapter 9 in Wolf's Integral Transforms in Science and Engineering.
I'm sure the Quantum mechanics interpretation/motivation will be valuable, especially for how these came about in the first place, as hinted by the post of @Spencer here.
It is not selfinverse: $F^4=I$ but $F^2 \neq I$.
$$ f=\sum_{n}(f,e_n)e_n,\;\;\;\; e_n(x)=\frac{1}{\sqrt{2\pi}}e^{inx} \\ g=\int_{s}(f,e_s)e_s,\;\;\;\; e_s(x)=\frac{1}{\sqrt{2\pi}}e^{isx} $$ See any similarities? The first is the expansion of a function in terms of the Fourier series on $[-\pi,\pi]$. The second is the expansion for the Fourier integral. The discrete Fourier transform is $\hat{f}(n)=(f,e_n)$, which is the coefficient function required to reconstruct $f$ in terms of the orthonormal basis $\{e_n \}_{n=-\infty}^{\infty}$. The Fourier transform looks more symmetric because both the coefficient function $\hat{f}(s)=(f,e_{s})$ and the reconstruction using this coefficient function are both integrals. But notice there's a difference, and the sign change is not just superficial.
The roles of the forward transform and the reverse integral reconstruction are quite different, even though they end up being related in this case. It's more of a coincidence that the computation of the coefficient function $\hat{f}(s)=(f,e_s)$ and the inverse reconstruction $\int_{s} \hat{f}(s)e_{s}$ are related. The original context definitely distinguished between the two, and they are distinguished by a sign change, too. Similarity between the forward transform and the reconstruction integral is interesting, but it seems to me that one loses context in viewing it that way. There are many other such integral eigenfunction expansions where it is important to distinguish between the coefficient computation and the integral reconstruction.
