I meant to submit this post first...
(Both posts refer to BICOHER.m)

Beware that there are multiple ways to normalize the bispectrum, and they are all referred to as the "bicoherence." The formula that this m-file uses is NOT the one that bounds the bicoherence between 0 and 1. If you are looking for that one, check out Kim and Powers' 1979 paper "Digital Bispectral Analysis and Its Applications to Nonlinear Wave Interactions." They provide an explanation of how to code that version, and it has worked reasonably well for me.

By the way, "segsamp" (actually implemented as "nsamp") is the number of samples per segment, that is, the number of points in each chunk of data that you're throwing into the bicoherence. You can use datasets of any length, and the program will typically calculate fft's with a length that is a power of two and is greater than or equal to "segsamp." It does this by adding zeros to the data record. However, you can specify "nfft" to tell the program how many points you want the fft to use. If you set "segsamp" and "nfft" to be equal, then you can avoid zero-padding. Powers of two are calculated faster, but MATLAB can handle fft's of any length.

As for the frequency axis, "waxis" is based on the number of points that are used in the fft. As a result, "0.5" on the axis corresponds to the highest frequency that can be resolved with that number of points, also known as the Nyquist frequency. Exactly how you want to scale the axis will depend on the nature of your data.

On an unrelated note, the demo gives me errors, too.