www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/bispect.m
function [Bspec,waxis] = bispect (ma, ar, nfft) %BISPECT Theoretical bispectrum of an ARMA model % [Bspec,waxis] = bispect (ma, ar, nfft) % ma - ma parameter vector % ar - ar parameter vector [default: [1]] % nfft - FFT length [default: 512] % Bspec - Bispectrum of ARMA model, an nfft by nfft array % with axes pointing down and right % waxis - frequency axis associated with the rows/columns of Bspec % Copyright (c) 1991-1999 by United Signals & Systems, Inc. and The Mathworks, Inc. All Rights Reserved. % $Revision: 1.7 $ % A. Swami January 20, 1995 % 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. % -------------- if (exist('ma') ~= 1) error('insufficient number of parameters') end q = length(ma)-1; if (exist('ar') ~= 1) ar = 1; end p = length(ar)-1; pq = max(p,q) + 1; if (exist('nfft') ~= 1) nfft = max(2^nextpow2(pq), 512); end nfft = max(nfft, 2^nextpow2(pq) ); Xf = freqz(ma,ar,nfft,'whole') ; Xfc = conj(Xf); Bspec = (Xf * Xf') .* hankel(Xfc,Xfc([nfft,1:nfft-1])); Bspec = fftshift(Bspec); % center the origin; normal axes if (rem(nfft,2) == 0) waxis = [-nfft/2:nfft/2-1]'/nfft; else waxis = [-(nfft-1)/2:(nfft-1)/2]'/nfft; end %contour(abs(Bspec),6,waxis,waxis),grid contour(waxis,waxis,abs(Bspec),6), grid on xlabel('f1'), ylabel('f2') x = [1 1 2/3 1 1/3 0 -1 -1 -2/3 -1 -1/3 0] ; y = [1 0 -1/3 -1 -2/3 -1 -1 0 1/3 1 2/3 1] ; hold on for k=1:12 plot([0 x(k)], [0 y(k)], '--') end hold off set(gcf,'Name','Hosa BISPECT')