bispecd双谱计算

估计高阶统计量的双谱程序，随机过程可用它的各阶矩来完整描述。与矩一样，可用累积量从时域来描述随机过程的统计特性；亦可用高阶谱(包括功率谱)从频域来描述随机过程的统计特性。由于三阶矩与三阶累积量相同，故三阶谱(也称双谱)定义为三阶矩函数(或称三阶自相关函数）的二重傅里叶变换。已测试过可以运行使用

function [Bspecwaxis] = bispecd （y  nfft wind nsamp overlap）
%BISPECD Bispectrum estimation using the direct （fft-based） approach.
%	[Bspecwaxis] = bispecd （y  nfft wind segsamp overlap）
%	y    - data vector or time-series数据向量或时间序列
%	nfft - fft length [default = power of two > segsamp]
%	wind - window specification for frequency-domain smoothing窗口规范
%	       if ‘wind‘ is a scalar it specifies the length of the side
%	          of the square for the Rao-Gabr optimal window  [default=5]
%	       if ‘wind‘ is a vector a 2D window will be calculated via
%	          w2（ij） = wind（i） * wind（j） * wind（i+j）
%	       if ‘wind‘ is a matrix it specifies the 2-D filter directly
%	segsamp - samples per segment [default: such that we have 8 segments]每段的样本
%	        - if y is a matrix segsamp is set to the number of rows行数
%	overlap - percentage overlap [default = 50] 重叠百分比
%	        - if y is a matrix overlap is set to 0.
%
%	Bspec   - estimated bispectrum: an nfft x nfft array with origin
%	          at the center and axes pointing down and to the right.
%	waxis   - vector of frequencies associated with the rows and columns
%	          of Bspec;  sampling frequency is assumed to be 1.

%  Copyright （c） 1991-2001 by United Signals & Systems Inc.
%       $Revision: 1.8 $
%  A. Swami   January 20 1993.

%     RESTRICTED RIGHTS LEGEND
% Use duplication or disclosure by the Government is subject to
% restrictions as set forth in subparagraph （c） （1） （ii） of the
% Rights in Technical Data and Computer Software clause of DFARS
% 252.227-7013.
% Manufacturer: United Signals & Systems Inc. P.O. Box 2374
% Culver City California 90231.
%
%  This material may be reproduced by or for the U.S. Government pursuant
%  to the copyright license under the clause at DFARS 252.227-7013.

% --------------------- parameter checks -----------------------------

    [ly nrecs] = size（y）;
    if （ly == 1） y = y（:）;  ly = nrecs; nrecs = 1; end

    if （exist（‘nfft‘） ~= 1）            nfft = 128; end
    if （exist（‘overlap‘） ~= 1）      overlap = 50;  end
    overlap = min（99max（overlap0））;
    if （nrecs > 1）                  overlap =  0;  end
    if （exist（‘nsamp‘） ~= 1）          nsamp = 0;   end
    if （nrecs > 1）                    nsamp = ly;  end

    if （nrecs == 1 & nsamp <= 0）
       nsamp = fix（ly/ （8 - 7 * overlap/100））;
    end
    if （nfft  < nsamp）   nfft = 2^nextpow2（nsamp）; end

    overlap  = fix（nsamp * overlap / 100）;             % added 2/14
    nadvance = nsamp - overlap;
    nrecs    = fix （ （ly*nrecs - overlap） / nadvance）;


% ------------------- create the 2-D window -------------------------
  if （exist（‘wind‘） ~= 1） wind = 5; end
  [mn] = size（wind）;
  window = wind;
  if （max（mn） == 1）     % scalar: wind is size of Rao-Gabr window
     winsize = wind;
     if （winsize < 0） winsize = 5; end        % the window