EDIT: Upon reading the question, I leapt to the conclusion that the asker was curious about the Fourier transforming properties of lenses. I am embarrassed to say that I did not immediately read the article. Upon doing so, I realize that the question is indeed about an all-optical FFT for temporal, not spatial frequencies. I've left my old answer below, but I've added a more appropriate answer here.
The article you linked is, like most popular reporting about science, very vague. I've done some checking, and it appears to me that the journal article being reported on is this one:
D. Hillerkuss, M. Winter, M. Teschke, A. Marculescu, J. Li, G. Sigurdsson, K. Worms, S. Ben Ezra, N. Narkiss, W. Freude, and J. Leuthold, "Simple all-optical FFT scheme enabling Tbit/s real-time signal processing," Opt. Express 18, 9324-9340 (2010)
The authors describe a method for rapidly isolating signals sent over optical fibers on different wavelengths. This basic idea has been used for quite a long time, and the simpler approach is called Wavelength Division Multiplexing (WDM). The idea is simple enough -- encode optical signals on light of different wavelengths, send them over a fiber, and use a prism or diffraction grating to separate the signals at the receiving end.
However, there is a limited range of wavelengths that are useful for transmission over long distances in optical fibers, so the more tightly spaced the signal wavelengths are, the more data can ultimately be sent. Modern techniques use the very closely spaced spectral lines produced by mode-locked lasers. This technique is called Orthogonal Frequency Division Multiplexing (at least that is what I gather. I am not an expert on telecom technology). At this level, it becomes impractical to isolate the signals with a diffraction grating, as the spatial separation between them is too small. Current technology apparently deals with this by digitizing the optical signals and performing Fourier analysis off line (i.e. not in real time) which limits the data rate.
The optical FFT developed by the authors is a solution to this problem. They have developed an optical apparatus which takes the incoming signal and splits it multiple times. Each copy of the signal is delayed by a specific amount, and then fed through a series of fiber-optic junctions which combine or phase shift as needed to implement a shift-and-add FFT algorithm. The key figure from the paper is reproduced below (this is allowed, I hope?)

The all-optical computation is carried out in real-time, allowing very high data rates.
Old Answer
Under the conditions that the medium of propagation is isotropic, and the light is monochromatic and coherent, the electric and magnetic fields along any chosen vector will differ by at most a proportionality constant, so the field can be expressed as a complex scalar field which specifies the magnitude and phase of the EM wave. This is called the "scalar approximation" and provides a mathematical basis for an approximate description of wave optics which applies in the vast majority of situations.
Using this scalar approximation, an optical field at one plane can be used to compute the field after some propagation distance. In general, this can be done using the Rayleigh-Sommerfeld Diffraction Integral:
$$ U(P_0) = \frac{1}{i \lambda} \int\int_\Sigma U(P_1) \frac{\mathrm{exp}(i k r_{01})}{r_{01}} \cos{\theta} \;\; \mathrm{ d}s$$
Where $U(P_0)$ is the optical field at plane $P_0$, $\lambda$ is the wavelength, $k$ is the waveneumber, $r_{01}$ is the vector from the point being integrated over in $P_1$ to the point being calculated in $P_0$. The region of integration, $\Sigma$, is the aperture in the input plane.
In most cases this integral cannot be solved analytically, and direct numerical solutions are very time consuming. However, under certain conditions (which are easily obtained in the real world) this integral can be approximated by Fraunhofer Diffraction or Fresnel Diffraction. I will limit myself here to the latter for simplicity, but the two situations are computationally very similar, so the answer to your question would be almost the same if I presented Fraunhofer Diffraction instead.
Fresnel Diffraction is obtained by expanding $r_{01}$ by a Taylor series and keeping terms up to second order. It is accurate when the propagation distance is sufficiently larger than the aperture. The exact definition of this condition varies from author to author, and results do not depend strongly on where this limitation is chosen. Under this assumption, the diffraction integral is reduced to:
$$ U(x,y) \propto \int\int^{\infty}_{-\infty} U(\xi, \eta) \mathrm{exp}\left[\frac{i \pi}{\lambda z}((x-\xi)^2 + (y-\eta)^2)\right] \;\; \mathrm{d}\xi \; \mathrm{d}\eta $$
This is in the form of a convolution, so convolution theorem allows this to be expressed as a Fourier transform as well:
$$ U(x,y) \propto \int\int^{\infty}_{-\infty} U(\xi, \eta) \mathrm{exp}\left[\frac{i k}{\lambda z} (\xi^2 +\eta^2)\right] \mathrm{exp}\left[\frac{-i 2 \pi}{\lambda z} (x\xi + y\eta)\right] \;\; \mathrm{d}\xi \; \mathrm{d}\eta $$
This would be an equality rather than a proportionality equation, but I have left of some complex factors of unity magnitude for simplicity. My notation has changed here to a less compact and more practical form because this integral is usually intended to be used for actual computation, as opposed to mathematical derivations; $x$ and $y$ are output coordinates, and $\xi$ and $\eta$ are input plane coordinates.
As you can see, this is simply a Fourier transform of $U(\xi, \eta)$, with some appropriately chosen scaling applied (specifically, $f_X \rightarrow \frac{x}{\lambda z}$, and $f_Y \rightarrow \frac{y}{\lambda z}$. $f_X$ and $f_Y$ are called "spatial frequencies" to avoid confusion with the more common temporal frequencies that we talk about in the context of Fourier analysis of time series.
As might be obvious by the scaling I described, this integral is easiest to compute in terms of the angular position in the output plane ($\frac{x}{z}$) when z is very large. In other words, if an optical field is allowed to propagate to infinity, its angular distribution will be the Fourier transform of the initial field.
In the case of a lens of focal length $f$, let us modify the initial field $U(\xi, \eta)$ by a quadratic phase term $\mathrm{exp}\left[\frac{-i\pi}{\lambda f}(\xi^2 + \eta^2)\right]$. We are interested in the field at the focal plane of the lens, so we will propagate a distance $z = f$.
Immediately, by looking at the Fourier transform formulation of the diffraction integral, it is obvious that the quadratic phase term inside the integral is exactly opposed by the phase term due to the lens (note that $k=\frac{2\pi}{\lambda}$). The resulting field is simply a Fourier transform of the input field, expressed in spatial coordinates rather than angular coordinates.