We present an efficient separable approach to the estimation and reconstruction of the bispectrum and the trispectrum from observational (or simulated) large scale structure data. This is developed from general CMB (poly-)spectra methods which exploit the fact that the bispectrum and trispectrum in the literature can be represented by a separable mode expansion which converges rapidly (with $n_textrm{max}={cal{O}}(30)$ terms). With an effective grid resolution $l_textrm{max}$ (number of particles/grid points $N=l_textrm{max}^3$), we present a bispectrum estimator which requires only ${cal O}(n_textrm{max} times l_textrm{max}^3)$ operations, along with a corresponding method for direct bispectrum reconstruction. This method is extended to the trispectrum revealing an estimator which requires only ${cal O}(n_textrm{max}^{4/3} times l_textrm{max}^3)$ operations. The complexity in calculating the trispectrum in this method is now involved in the original decomposition and orthogonalisation process which need only be performed once for each model. However, for non-diagonal trispectra these processes present little extra difficulty and may be performed in ${cal O}(l_textrm{max}^4)$ operations. A discussion of how the methodology may be applied to the quadspectrum is also given. An efficient algorithm for the generation of arbitrary nonGaussian initial conditions for use in N-body codes using this separable approach is described. This prescription allows for the production of nonGaussian initial conditions for arbitrary bispectra and trispectra. A brief outline of the key issues involved in parameter estimation, particularly in the non-linear regime, is also given.