You start out with the Helmholtz equation that holds true in a homogeneous lossless medium for any rectangular component of both the $E$ and $H$ fields:
$$\nabla^2 U + \kappa^2 U = 0 \tag{1}\label{1}$$
That is $U$ may stand for any of the $E_x, E_y,...,H_z$. A general solution of $\eqref{1}$ is the plane wave $e^{\mathfrak {j} \kappa \hat {\mathbf{k}} \cdot \mathbf{r}}$ where $\hat {\mathbf{k}}$ is the unit vector in the direction of propagation and $\kappa = \omega/c =2\pi/\lambda$.
Now assume that a wave, say $U(x,y,z)$ propagates *mostly* along the $z$ axis (paraxial approximation) and is incident upon the $z=0$ plane. Decompose this field at an arbitrary $z$ into plane wave components using the 2D Fourier transform:
$$U(x,y,z)=\int \int dudv A_z(u,v) e^{\mathfrak j \kappa(xu+yv)}\tag{2}\label{2}$$
Here $$dU(x,y,z)= A_z(u,v) e^{\mathfrak j \kappa (xu+yv)} dudv$$ is a differential plane wave moving in the direction represented by the direction cosines $u,v,w$ where $u^2+v^2+w^2=1$, i.e., $w=\sqrt{1-u^2-v^2}$.
Next we wish to find the relationship between the field (component) $U(x,y,z)$ at $z=0$ and at arbitrary $z$. From linearity $ A_z(u,v) = A_0(u,v) e^{\mathfrak{j}\kappa w}$. The inverse Fourier transform gives:
$$A_0(x,y)=b_1\int \int dx'dy' U(x',y',0) e^{-\mathfrak j \kappa(x'u+y'v)}\tag{3}\label{3}$$ for some constant $b$ (in fact $b_1=1/(2\pi)^2$
Again because of the linear shift invariant nature of hte optical (EM) system the *output* field at $z$ is a convolution of the *input* field at $z=0$ and kernel, say, $G(x,y)$, that is
$$U(x,y,z) =\int \int dx'dy' U(x',y',0) G(x-x',y-y')\tag{4}\label{4}$$
and
$$G(x,y) = \int \int du dv e^{\mathfrak j \kappa w} e^{\mathfrak j \kappa (ux+vy)} \tag{5}\label{5}$$ where $w=\sqrt{1-u^2-v^2}$.
When polar coordinates are introduced $u=R cos\theta$ and $v=R sin \theta$ after much pain we get (see the Erdelyi Higher Transcendental Functions, 4.52)
$$G(R cos\theta,R sin \theta) = \frac{1}{\mathfrak j \lambda}\frac{e^{\mathfrak j \kappa \sqrt{R^2+z^2}}}{\sqrt{R^2 + z^2}} \frac{z}{\sqrt{R^2 + z^2}}\left( 1-\frac{1}{\mathfrak j \kappa \sqrt{R^2 + z^2}}\right) \tag{6}\label{6}$$
Now in the far field (Fraunhofer) we have
- when $z \gg R$ the term in parentheses is approximated as $1$
- Define $cos\phi = \frac{z}{\sqrt{R^2+z^2}}$, this angle $\phi$ is the direction the elementary wavelet has with the axis $z$ as it propagates, and $cos\phi$ is the obliquity factor introduced first by Fresnel. Taking $\phi \approx 0$ that is $cos\phi \approx 1$ is the paraxial approximation.
- The term $\sqrt{R^2+z^2}$ is the distance between the point $x',y'$ in the $z=0$ plane and the field point $x,y,z$, and in the paraxial regime this can be approximated as $z$.
Therefore
$$G(R cos\theta,R sin \theta)\approx \tilde G(R) = \frac{1}{\mathfrak j \lambda}\frac{e^{\mathfrak j \kappa \sqrt{R^2+z^2}}}{z} \tag{7}\label{7}$$
This gives us the field at $z$ approximately by calcuating the convolution integral as
$$ U(x,y,z)\approx \frac{b_2}{z}\int \int dx'dy'U(x',y',0)e^{\mathfrak j \kappa \sqrt{R^2+z^2}} \tag{8}\label{8}$$
Finally, in the Fraunhofer far-field we also have the approximation
$$\sqrt{R^2+z^2} = \sqrt{z^2+(x-x')^2+(y-y')^2} \\
\approx z+ \frac{1}{2z}(x-x')^2+\frac{1}{2z}(y-y')^2,$$ and thus
$$ U(x,y,z)\approx \frac{b_2 e^{\mathfrak j \kappa z}}{z}\int \int dx'dy'U(x',y',0)e^{\mathfrak j \kappa (\frac{1}{2z}(x-x')^2+\frac{1}{2z}(y-y')^2)} \tag{9}\label{9}$$
The result $\eqref{9}$ is the diffraction integral, the input (incident) field $U(x,y,0)$is at $z=0$ represents any component of $E$ or $H$ and is supposed to be known. It propagates in the free half-space $z>0$ and via this convolution integral one can calculate the same field component at any point $x,y,z$.
If you wish so you can also write the $\eqref{9}$ convolution integral as a Fourier transform by expanding the quadratic terms in the exponent:
$$ U(x,y,z)\approx \frac{b_2 e^{\mathfrak j \kappa (z+x^2/(2z) +y^2/(2z))}}{z}\int \int dx'dy' e^{\mathfrak j \tfrac{\kappa}{2z} (x'^2 +y'^2)} U(x',y',0)e^{-\mathfrak j \kappa (xx'+yy')/z} \tag{10}\label{10}$$
but notice that $\eqref{10}$ is not the Fourier transform of $U(x,y,0)$ but instead there is a quadratic phase factor also involved, so that aside from some unit modulus phase factor outside the integral this is the Fourier transform of $e^{\mathfrak j \kappa(x^2 +y^2)/(2z)} U(x,y,0)$. Now if the extent of the diffracting object at $z=0$ is such that $\tfrac{\kappa}{2z}(x^2 +y^2) \ll 1 $ or $x^2 +y^2 \ll z \lambda $ then the quadratic phase factor can be neglected and we do have the Fourier transform of the incident field.
[1] Stark: Application of Optical Fourier Transforms, Academic Press 1982
Best Answer
To give an unmathematical catchy answer, let's look at Fraunhofer diffraction in double slit experiment. Interference at the observation plane depends on slit parameter $d$. What is the frequency of slits? E.g. $1\,\text{mm}\frac{1}{d}$: number of slits per length. Concluding frequiency in the setup. The following argumentation links this frequency to the fourier transform. The physical significance is in the real optics setup. The setup is easier described, when transformed in fourier space.
The double slit link on high school level above gives all the math without integrals. After visualizing the following concept, you see that the integrals are just math to convert the diffraction pattern in fourier space via fourier transformation.
Using trigonometry first compute phase difference $\Delta\phi(\theta,d)$. Go deeper in this concept using a sketch to visualize phase difference of $n\cdot\lambda$, $n\in\mathbf{N}$ as bright maxima in diffraction pattern. There is no magic in the next step. It's just a another point of view: Try to grasp $\frac{1}{d}$ as a parameter on its own: $\Delta\phi(\theta,\frac{1}{d})$.
Fourier space is a synonym for frequency domain. Acoustics examples are given in Eichenlaubs SE answer and ptomatoes optics explanation.
No literature: Calculate and understand yourself with the links above.