Interactions among neural signals in different frequency bands have become a focus of strong interest in neuroscience. Bispectral analysis, a type of higher order spectral analysis, provides us with the ability to investigate such nonlinear interactions. Based on the fact that the bispectrum of a linear Gaussian process is zero, a surrogate data method was proposed to test the null hypothesis that the original data were generated by a linear Gaussian process. The method was first tested on two simulation examples. It was then applied to local field potential recordings from a monkey performing a visuomotor task. The analysis reveals nonzero bispectra for beta and gamma band activities in the premotor cortex. The rigorous statistical framework proves essential in establishing these results.