The proof I came up with for the $L^2$ isometry of the Fourier transform on $L^1(\Bbb R)\cap L^2(\Bbb R)$ hinges on the $L^2$ completeness of the Hermite-Gauss functions. It is well-known that these are eigenfunctions of the Fourier transform however the proofs using the Hermite-Gauss functions immediately jump to an $L^2$ theory and neglect the integral operator aspect of the Fourier transform, i.e. the $L^1(\Bbb R)$ aspect. The convention for the Fourier transform I am using is
$$\mathcal{F}f(k) = \frac{1}{\sqrt{2\pi}}\int_{\Bbb R} e^{-ikx}f(x)\,dx.$$
The Hermite polynomials $H_n$ are given by the recursion relation
$$ H_n(x) = (-1)^n e^{x^2}\frac{d^n}{dx^n}e^{-x^2}.$$
The Hermite-Gauss functions are given by $H_n(x)e^{-\frac{x^2}{2}}$. It is not hard to prove that these are eigenfunctions of the Fourier transform with the following recursion relation:
$$ \left(x-\frac{d}{dx}\right)\left(H_n(x)e^{-\frac{x^2}{2}}\right) = H_{n+1}(x)e^{-\frac{x^2}{2}}.$$
We proceed by induction on $n$. It is well-known that the Gaussian is an eigenfunction of the Fourier transform with eigenvalue $1$ and this takes care of the base case.
Suppose then that $H_n(x)e^{-\frac{x^2}{2}}$ is an eigenfunction with eigenvalue $(-i)^n$, then we wish to show that $H_{n+1}(x)e^{-\frac{x^2}{2}}$ is also an eigenfunction with eigenvalue $(-i)^{n+1}$. Taking the Fourier transform of this, we get
\begin{align}
\int_{\Bbb R} e^{-ikx} H_{n+1}(x)e^{-\frac{x^2}{2}}\,dx &= \int_{\Bbb R} e^{-ikx} \left(x-\frac{d}{dx}\right)\left(H_n(x)e^{-\frac{x^2}{2}}\right)\,dx \\
&= \left(i\frac{d}{dk}-ik\right)\int_{\Bbb R} e^{-ikx}H_n(x)e^{-\frac{x^2}{2}}\,dx \\
&= -i(-i)^n\left(k-\frac{d}{dk}\right)\left(H_n(k)e^{-\frac{k^2}{2}}\right) \\
&= (-i)^{n+1} H_{n+1}(k)e^{-\frac{k^2}{2}}.
\end{align}
Thus the Hermite-Gauss functions are eigenfunctions of the Fourier transform with eigenvalues $\pm 1,\pm i$. It remains to see that they are mutually orthogonal. Suppose then that $m < n$.
$$\int_{\Bbb R}H_m(x)e^{-\frac{x^2}{2}} H_n(x)e^{-\frac{x^2}{2}}\,dx = \int_{\Bbb R} H_m(x) \frac{d^n}{dx^n} e^{-x^2}\,dx.$$
Since $H_m$ is an $m$th degree polynomial, an integration by parts $m$ times will leave a constant multiple of an $n-m$th derivative of $e^{-\frac{x^2}{2}}$ and so the integral is zero, i.e. if $m\neq n$, $H_m(x)e^{-\frac{x^2}{2}}$ and $H_n(x)e^{-\frac{x^2}{2}}$ are orthogonal.
Denoting the (normalized) Hermite-Gauss functions by $\psi_n$, from the above property, $\|\mathcal{F}\psi_n\| = \|(-i)^n \psi_n\| = \|\psi_n\|$. As noted above, the Hermite-Gauss functions are complete in $L^2(\Bbb R)$ which can be proved without the Fourier transform (though most modern proofs do use the Fourier transform).
Consider then $f\in L^1(\Bbb R)\cap L^2(\Bbb R)$. On this space, the Fourier transform still acts as an integral operator. Moreover,
$$\int_{\Bbb R^2}|e^{-ikx}f(x)\psi_n(k)|\,dx\,dk < \infty$$
since $f,\psi_n\in L^1(\Bbb R)$. Making use of Fubini's theorem, we have that
$$\int_{\Bbb R} \mathcal{F}f(k)\psi_n(k)\,dx = \int_{\Bbb R} f(x)\mathcal{F}\psi_n(x) = (-i)^n \int_{\Bbb R} f(x)\psi_n(x)\,dx.$$
Thus
$$\sum_{n=0}^{\infty} |\langle \mathcal{F}f,\psi_n\rangle|^2 = \sum_{n=0}^{\infty} |(-i)^n\langle f,\psi_n\rangle|^2 = \|f\|^2 < \infty$$
and so $\mathcal{F}f\in L^2(\Bbb R)$ and moreover $\|\mathcal{F}f\| = \|f\|.$ Thus $\mathcal{F}(L^1(\Bbb R)\cap L^2(\Bbb R))\subseteq L^2(\Bbb R)$ as claimed.