www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/cum3est.m
function y_cum = cum3est (y, maxlag, nsamp, overlap, flag, k1) %CUM3EST Third-order cumulants. % Should be invoked via "CUMEST" for proper parameter checks % y_cum = cum3est (y, maxlag, samp_seg, overlap, flag, k1) % y_cum = cum3est (y, maxlag, samp_seg, overlap, flag, k1) % y: input data vector (column) % maxlag: maximum lag to be computed % samp_seg: samples per segment % overlap: percentage overlap of segments % flag : 'biased', biased estimates are computed [default] % 'unbiased', unbiased estimates are computed. % k1: the fixed lag in c3(m,k1): see below % y_cum: estimated third-order cumulant, % C3(m,k1) -maxlag <= m <= maxlag % Copyright (c) 1991-1999 by United Signals & Systems, Inc. and The Mathworks, Inc. All Rights Reserved. % $Revision: 1.4 $ % A. Swami January 20, 1993 % Modified Jan 20, 94 to handle complex case properly. % c3(i,j) := E x^*(n) x(n+i) x(n+j) (x assumed zero mean) % 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. % c3(i,j) := E x^*(n) x(n+i) x(n+j) (x assumed zero mean) %---------------- Parameter checks done by CUMEST -------------- [n1,n2] = size(y); N = n1*n2; minlag = -maxlag; overlap = fix(overlap/100 * nsamp); nrecord = fix( (N - overlap)/(nsamp - overlap) ); nadvance = nsamp - overlap; y_cum = zeros(maxlag-minlag+1,1); ind = (1:nsamp)'; nlags = 2 * maxlag + 1; zlag = 1 + maxlag; if (flag(1) == 'b' | flag(1) == 'B') scale = ones(nlags,1)/nsamp; else lsamp = nsamp - abs(k1); scale = [lsamp-maxlag:lsamp,lsamp-1:-1:lsamp-maxlag]'; [m2,n2] = size(scale); scale = ones(m2,n2) ./ scale; end for i=1:nrecord x = y(ind); x = x(:) - mean(x); % make sure we have a col vec cx = conj(x); z = x*0; % create the "IV" matrix: offset for second lag if (k1 >= 0) z(1:nsamp-k1) = x(1:nsamp-k1,:) .* cx(k1+1: nsamp,:); else z(-k1+1:nsamp) = x(-k1+1:nsamp) .* cx(1:nsamp+k1); end % compute third-order cumulants y_cum(zlag) = y_cum(zlag) + z' * x; for k = 1:maxlag y_cum(zlag-k) = y_cum(zlag-k) + z([k+1:nsamp])' * x([1:nsamp-k]); y_cum(zlag+k) = y_cum(zlag+k) + z([1:nsamp-k])' * x([k+1:nsamp]); end ind = ind + nadvance; end y_cum = y_cum .* scale / nrecord; return