  Measurements of line-of-sight dependent clustering via the galaxy power
spectrum's multipole moments constitute a powerful tool for testing theoretical
models in large-scale structure. Recent work shows that this measurement,
including a moving line-of-sight, can be accelerated using Fast Fourier
Transforms (FFTs) by decomposing the Legendre polynomials into products of
Cartesian vectors. Here, we present a faster, optimal means of using FFTs for
this measurement. We avoid redundancy present in the Cartesian decomposition by
using a spherical harmonic decomposition of the Legendre polynomials.
Consequently, our method is substantially faster: a given multipole of order
$\ell$ requires only $2\ell+1$ FFTs rather than the $(\ell+1)(\ell+2)/2$ FFTs
of the Cartesian approach. For the hexadecapole ($\ell = 4$), this translates
to $40\%$ fewer FFTs, with increased savings for higher $\ell$. The reduction
in wall-clock time enables the calculation of finely-binned wedges in
$P(k,\mu)$, obtained by computing multipoles up to a large $\ell_{\rm max}$ and
combining them. This transformation has a number of advantages. We demonstrate
that by using non-uniform bins in $\mu$, we can isolate plane-of-sky (angular)
systematics to a narrow bin at $\mu \simeq 0$ while eliminating the
contamination from all other bins. We also show that the covariance matrix of
clustering wedges binned uniformly in $\mu$ becomes ill-conditioned when
combining multipoles up to large values of $\ell_{\rm max}$, but that the
problem can be avoided with non-uniform binning. As an example, we present
results using $\ell_{\rm max}=16$, for which our procedure requires a factor of
3.4 fewer FFTs than the Cartesian method, while removing the first $\mu$ bin
leads only to a 7% increase in statistical error on $f \sigma_8$, as compared
to a 54% increase with $\ell_{\rm max}=4$.
