  The anisotropy of the redshift space bispectrum depends upon the orientation
of the triangles formed by three $\vec{k}$ modes with respect to the line of
sight. For a triangle of fixed size ($k_1$) and shape ($\mu,t$), this
orientation dependence can be quantified in terms of angular multipoles
$B_l^m(k_1,\mu,t)$ which contain a wealth of cosmological information. We
propose a fast and efficient FFT-based estimator that computes bispectrum
multipole moments $B_l^m$ of a 3D cosmological field for all possible $l$ and
$m$ (including $m\neq 0$). The time required by the estimator to compute all
multipoles from a gridded data cube of volume $N_g^3$ scales as $\sim N_g^3
\log{(N_g)}$ in contrast to the direct computation technique which requires
time $\sim N_g^6$. Here, we demonstrate the formalism and validate othe
estimator using a simulated non-Gaussian field for which the analytical
expressions for all bispectrum multipoles are known. The estimated results are
found to be in good agreement with the analytical predictions for all $16$
non-zero multipoles (up to $\ell= 6, m=6$). We expect the $m \neq 0$ bispectrum
multipoles to significantly enhance the information available from galaxy
redshift surveys, and future redshifted 21-cm observations.
