  Fourier extension is an approximation scheme in which a function on an
arbitary bounded domain is approximated using a classical Fourier series on a
bounding box. On the smaller domain the Fourier series exhibits redundancy, and
it has the mathematical structure of a frame rather than a basis. It is not
trivial to construct approximations in this frame using function evaluations in
points that belong to the domain only, but one way to do so is through a
discrete least squares approximation. The corresponding system is extremely
ill-conditioned, due to the redundancy in the frame, yet its solution via a
regularized SVD is known to be accurate to very high (and nearly spectral)
precision. Still, this computation requires ${\mathcal O}(N^3)$ operations. In
this paper we describe an algorithm to compute such Fourier extension frame
approximations in only ${\mathcal O}(N^2 \log^2 N)$ operations for general 2D
domains. The cost improves to ${\mathcal O}(N \log^2N)$ operations for simpler
tensor-product domains. The algorithm exploits a phenomenon called the plunge
region in the analysis of time-frequency localization operators, which
manifests itself here as a sudden drop in the singular values of the least
squares matrix. It is known that the size of the plunge region scales like
${\mathcal O}(\log N)$ in one dimensional problems. In this paper we show that
for most 2D domains in the fully discrete case the plunge region scales like
${\mathcal O}(N \log N)$, proving a discrete equivalent of a result that was
conjectured by Widom for a related continuous problem. The complexity estimate
depends on the Minkowski or box-counting dimension of the domain boundary, and
as such it is larger than ${\mathcal O}(N \log N)$ for domains with fractal
shape.
