clear;
%	BISP3CUM Auto bispectrum/3rd order cumulant
%
%	[bisp,freq,cum,lag]=bisp3cum(signal,samprate,maxlag,window,scale)
%
%	The maxlag*2+1 x maxlag*2+1 element auto bispectrum and 3rd order cumulant matrices
%	and maxlag*2+1 element frequency and lag vectors are computed from the signal
%	matrix containing samples in rows and records in columns, signal sample rate and
%	maximum lag scalars, and lag window function and scale strings.
%
%	If unspecified, the signal matrix is entered after the prompt from the keyboard,
%	and the default assignments samprate=1 and maxlag=0 are used.  The window and scale
%	strings specify lag window and scale matrix computation, according to:
%
%	window = 'none', 'n', or unspecified does not compute a window
%	       = 'uniform' or 'u' computes the uniform hexagonal window
%	       = 'sasaki' or 's' computes the sasaki window
%	       = 'priestley' or 'p' computes the priestley window
%	       = 'parzen' or 'pa' computes the parzen window
%	       = 'hamming' or 'h' computes the hamming window
%	       = 'gaussian' or 'g' computes the gaussian distribution window
%	       = 'daniell' or 'd' computes the daniell window
%
%	scale  = 'biased', 'b', or unspecified computes the biased estimate
%	       = 'unbiased' or 'u' computes the unbiased estimate

%	Implemented using MATLAB 5.3.1 and additional functions:
%
%	mat=toep(column,row)
%	wind=lagwind(lag,window)
%
%	Implementation:
%
%	cum(k,l) = sum_{n=0}^{N-1} conj(signal(n))*signal(n+k)*signal(n+l)/N
%
%	k,l = {-maxlag,...,-1,0,1,...,maxlag}, n = {0,1,...,N-1}
%
%	bisp=fftshift(fft2(ifftshift(cum.*wind)))
%
%	Example:
%
%	?[b,f,c,l]=bisp3cum([1-i -1+i],1,1)
%
%	b =
%
%	  -5.1962 - 5.1962i        0            -0.0000 - 0.0000i
%	        0                  0                  0          
%	  -0.0000 - 0.0000i   0.0000 + 0.0000i   5.1962 + 5.1962i
%
%	f =
%
%	   -0.5000         0    0.5000
%
%	c =
%
%	  -1.0000 + 1.0000i   1.0000 - 1.0000i        0          
%	   1.0000 - 1.0000i        0            -1.0000 + 1.0000i
%	        0            -1.0000 + 1.0000i   1.0000 - 1.0000i
%
%	l =
%
%	    -1     0     1
%
%	References:
%
%	C. L. Nikias, A. P. Petropulu, Higher-Order Spectra Analysis:  A Nonlinear Signal
%	Processing Framework, PTR Prentice Hall, Englewood Cliffs, NJ, 1993.
%
%	T. S. Rao, M. M. Gabr, An Introduction to Bispectral Analysis and Bilinear Time
%	Series Models, Lecture Notes in Statistics, Volume 24, D. Brillinger, S. Fienberg,
%	J. Gani, J. Hartigan, K. Krickeberg, Editors, Springer-Verlag, New York, NY, 1984.
%
%	Copyright (c) 2000
%	Tom McMurray
%	mcmurray@teamcmi.com

%	assign default input parameters

clear;

nt=128;
f0=10;
fs=4*f0;
nf=128;
for i=1:nt;
% x1(i)=sin(2*pi*f0*i/16);
  x1(i)=sign(sin(2*pi*f0*i/16));
 % x1(i)=sign(randn(1));
end
x=x1+randn(1,nt);

[bisp,freq,cum,lag]=bisp3cum(x,fs,nt-1, 'none','b');
 figure(2)
 mesh(abs(bisp))
figure(3)
mesh(cum)


 %figure(4)
% plot(lag)