www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/maorder.m
function qopt = maorder(x,qmin,qmax,pfa, flag) %MAORDER MA order determination % qopt = maorder(x, qmin, qmax,pfa, flag) % x - time series (must be a vector) % qmin - minimum MA order [default = 0] % qmax - maximum MA order [default = 10] % pfa - probability of false alarm [default = 0.05] % flag - if non-zero (default is 1), various values are displayed % qopt - estimated MA order % Copyright (c) 1991-1999 by United Signals & Systems, Inc. and The Mathworks, Inc. All Rights Reserved. % $Revision: 1.4 $ % 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. % -------------------- Parameter checks -------------------------- if (exist('x') ~= 1) error(' vector x must be specified') end [m,n] = size(x); if (min(m,n) ~= 1) error(' x must be a vector, not a matrix') end if (m==1) x = x(:); m=n; n=1; end nsamp = m; if (exist('qmin') ~= 1) qmin = 0; end if (exist('qmax') ~= 1) qmax = 10; end if (qmax < qmin) qmax = qmin + 1; end if (exist('pfa') ~= 1) pfa = 0.05; end if (exist('flag') ~= 1) flag = 1; end cvec = cumest (x, 3, qmax+1,nsamp,0,'unbiased'); % the c(m,0) slice cvec = cvec(qmax+3:2*qmax+3); % c(m,0) m=1, ... ,qmax+1 var = zeros(qmax+1,1); ord = zeros(qmax+1,1); thres = zeros(qmax+1,1); xsq = x.^2; z = zeros(nsamp,1); inveps = erfinv(1-pfa); % ---------------------------------------------------------------- for q=qmin:qmax m=q+1; z = xsq(1:nsamp-m) .* x(m+1:nsamp); % z(n) = y^2(n) * y(n+m) rzz = cumest (z,2,2*qmax,nsamp-m,0,'unbiased'); % covariance of z(n) rzz = rzz + cvec(m)^2; % correlation of z(n) var(m) = sum(rzz) / (nsamp-m); % variance of c(q+1,0) thres(m) = inveps * sqrt(2*var(m)); % (1-pfa) sigma bounds ord(m) = (abs(cvec(m)) < thres(m)); % less ? then MA(q) end q = max(find(ord == 0)); % max m for which MA(m) hypothesis fails ind = find( ord(q+1:qmax+1) > 0); % now find where hypothesis holds qopt = min(ind) + q - 1; % fix the index if (isempty(qopt)) qopt = qmax+1; end if (flag) disp([' q var(cqk) cqk thres result ']) for q=qmin:qmax m = q+1; fprintf(' %3g %15.5e %15.5e',q, var(m), cvec(m)) fprintf(' %15.5e %1g \n', thres(m), ord(m)) end disp([' Estimated MA order is ', int2str(qopt) ]) end return