octave_packages/tsa-4.2.4/bispec.m

   1 function [BISPEC,BIACF,ACF] = bispec(Z,N);
   2 % Calculates Bispectrum
   3 % [BISPEC] = bispec(Z,N);
   4 %
   5 % Input:        Z    Signal
   6 %               N  # of coefficients
   7 % Output:       BiACF  bi-autocorrelation function = 3rd order cumulant
   8 %               BISPEC Bi-spectrum
   9 %
  10 % Reference(s):
  11 % C.L. Nikias and A.P. Petropulu "Higher-Order Spectra Analysis" Prentice Hall, 1993.
  12 % M.B. Priestley, "Non-linear and Non-stationary Time series Analysis", Academic Press, London, 1988.
  13
  14 %       $Id: bispec.m 5090 2008-06-05 08:12:04Z schloegl $
  15 %       Copyright (C) 1997-2003,2008 by Alois Schloegl <a.schloegl@ieee.org>
  16 %
  17 %    This program is free software: you can redistribute it and/or modify
  18 %    it under the terms of the GNU General Public License as published by
  19 %    the Free Software Foundation, either version 3 of the License, or
  20 %    (at your option) any later version.
  21 %
  22 %    This program is distributed in the hope that it will be useful,
  23 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25 %    GNU General Public License for more details.
  26 %
  27 %    You should have received a copy of the GNU General Public License
  28 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  29
  30
  31 P=N+1;
  32 ACF=zeros(1,N+1);
  33 BIACF=zeros(2*N+1,2*N+1);
  34
  35 Z=Z(:);
  36 M=size(Z,1);
  37 M1=sum(Z)/M;
  38 Z=Z-M1*ones(size(Z));
  39
  40 for K=0:N,
  41         jc2=Z(1:M-K).*Z(1+K:M);
  42         ACF(K+1)=sum(jc2)/M;
  43         for L = K:N,
  44                 jc3 = sum(jc2(1:M-L).*Z(1+L:M))/M;
  45                 BIACF(K+P,  L+P)  =jc3;
  46                 BIACF(L+P,  K+P)  =jc3;
  47                 BIACF(L-K+P, -K+P)=jc3;
  48                 BIACF(-K+P, L-K+P)=jc3;
  49                 BIACF(K-L+P, -L+P)=jc3;
  50                 BIACF(-L+P, K-L+P)=jc3;
  51         end;
  52 end;
  53
  54 BISPEC=fft2(BIACF,128,128);