bisp3cum.m

人们熟知的双谱（bispectrum）

function [bisp,freq,cum,lag]=bisp3cum(signal,samprate,maxlag,window,scale)

%	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:
%
%