  abridged] A method to rapidly estimate the Fourier power spectrum of a point
distribution is presented. This method relies on a Taylor expansion of the
trigonometric functions. It yields the Fourier modes from a number of FFTs,
which is controlled by the order N of the expansion and by the dimension D of
the system. In three dimensions, for the practical value N=3, the number of
FFTs required is 20. We apply the method to the measurement of the power
spectrum of a periodic point distribution that is a local Poisson realization
of an underlying stationary field. We derive explicit analytic expression for
the spectrum, which allows us to quantify--and correct for--the biases induced
by discreteness and by the truncation of the Taylor expansion, and to bound the
unknown effects of aliasing of the power spectrum. We show that these aliasing
effects decrease rapidly with the order N. The only remaining significant
source of errors is reduced to the unavoidable cosmic/sample variance due to
the finite size of the sample. The analytical calculations are successfully
checked against a cosmological N-body experiment. We also consider the initial
conditions of this simulation, which correspond to a perturbed grid. This
allows us to test a case where the local Poisson assumption is incorrect. Even
in that extreme situation, the third-order Fourier-Taylor estimator behaves
well. We also show how to reach arbitrarily large dynamic range in Fourier
space (i.e., high wavenumber), while keeping statistical errors in control, by
appropriately "folding" the particle distribution.
