function [sg, sl] = glstat(x,cparm,nfft)
%GLSTAT Detection statistics for Hinich's Gaussianity and linearity tests
% [sg, sl] = glstat(x,cparm,nfft)
% x - time series
% cparm - resolution (smoothing) parameter c;
% 0.5 < c < 1.0 [default = 0.51]
% increasing c reduces the variance,
% but increases the bias (window bandwidth).
% nfft - fft length [default = 128]
% sg - statistic for the Gaussianity test
% [observed chi_sq value, dg, PFA]
% sl - statistic for the linearity test
% [R_estimated, lambda, R_theoretical]
%
% Under the assumption of Gaussianity, the test statistics S is
% chi-squared distributed with df degrees of freedom.
% Under the assumption of non-zero skewness and linearity, the test
% statistic R should be approximately equal to the inter-quartile
% range of a chi-squared distribution with two degrees of freedom,
% and non-centrality parameter lambda.
%
% Copyright (c) 1991-98 by United Signals & Systems, Inc. and The MathWorks, Inc.
% $Revision: 1.4 $
% A. Swami January 20, 1993.
% A. Swami, October 14, 1994, Revised.
% A. Swami, February 23, 1995, added `N' printout in linearity test
% 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 (min(size(x)) ~= 1)
error('x: should be a vector')
end
nsamp = length(x);
if (exist('nfft') ~= 1) nfft = 128; end
nrecs = floor(nsamp/nfft);
nrecs = max(nrecs,1);
ksamp = min(nfft,nsamp);
if (nfft > nsamp)
disp([' glstat results are unreliable if zero-padding is severe: '])
disp([' fft length= ',int2str(nfft),' data length=',int2str(nsamp)])
end
if (exist ('cparm') ~= 1) cparm = 0.51; end
if (cparm >= 1.0)
error('cparm cannot be greater than or equal to 1.')
end
if (cparm <=0.5 & nrecs == 1)
error('cparm: for single segments: allowed range is (0.5,1.0) ')
end
M = fix(nfft^cparm); M = M + 1 - rem(M,2);
Msmuth = M;
% ----------------- estimate raw bispectrum ------------------------
mrow = fix(nfft/3)+1; ncol = fix(nfft/2)+1;
F = zeros(mrow,ncol);
S = zeros(nfft,1);
mask = hankel ([1:mrow],[mrow:mrow+ncol-1]);
ind = 1:ksamp;
for k=1:nrecs
y = x(ind);
xf = fft(y-mean(y), nfft);
xc = conj(xf);
S = S + xf .* xc; % power spectrum
F = F + xf(1:mrow) * xf(1:ncol).' .* ...
reshape (xc(mask), mrow, ncol) ;
ind = ind + ksamp;
end
F = F /(nfft*nrecs); % `Raw' bispectrum F as in Eq (2.3)
S = S /(nfft*nrecs); % `Raw' power spectrum
% ---------- zero out area outside principal domain ---------------
ind = (0:(mrow-1))';
F(1:mrow,1:mrow) = triu(F(1:mrow,1:mrow));
Q = ones(mrow, ncol);
Q(1:mrow,1:mrow) = triu(Q(1:mrow,1:mrow)) + eye(mrow);
% the 2f1 + f2 = 1 line:
r = ( rem(nfft,3) == 2); % ans oct 94: also -r next line
for k=mrow+1:ncol-r
j = k-mrow;
Q(mrow-2*j+2:mrow, k) = zeros(2*j-1,1);
F(:, k) = F(:,k) .* Q(:,k);
Q(mrow-2*j+1, k) = 2; % factor 2 on boundary
end
F = F(2:mrow,2:ncol); % in principal domain, no dc terms
Q = Q(2:mrow,2:ncol);
mrow = mrow-1;
ncol = ncol-1;
% -------- smooth the estimated spectrum and bispectrum ------------
% Msmuth x Msmuth box-car smoother
% only every Msmuth term in the smoothed output is retained
% "independent" estimates
M = Msmuth;
mby2 = (M+1)/2;
m1 = rem(mrow,M); m2=rem(ncol,M);
m1 = - m1 + M * (m1 ~= 0);
m2 = - m2 + M * (m2 ~= 0);
F = [F, zeros(mrow,m2); zeros(m1,ncol+m2)];
Q = [Q, zeros(mrow,m2); zeros(m1,ncol+m2)];
k = size(F)/Msmuth;
k1 = k(1); k2 = k(2);
% Apply the box car smoother
% can replace the next ten lines or so with
% B =kron(eye(k1),ones(1,Msmuth)) * F * kron(eye(k2),ones(Msmuth,1))/Msmuth^2;
% Q =kron(eye(k1),ones(1,Msmuth)) * Q * kron(eye(k2),ones(Msmuth,1));
ind = 1:Msmuth;
B = zeros(k1,k2); Q1 = B;
for i=1:k1
for j=1:k2
t = F( (i-1)*Msmuth+ind, (j-1)*Msmuth+ind );
B(i,j) = mean(t(:));
t = Q( (i-1)*Msmuth+ind, (j-1)*Msmuth+ind );
Q1(i,j) = sum(t(:));
end
end
Q = Q1;
% --------------
% At this point B corresponds to B in eq (2.6) of Hinich's paper
% and Q corresponds to the definition following eq (2.8)
% --------------
M = Msmuth;
S = conv(S, ones(M,1))/M;
S = S(M+1:M:M+nfft-1);
S1 = B*0;
% S1 = S(1:k1) * S(1:k2)' .* hankel(S(1:k1),S(k1:k1+k2-1));
S1 = S(1:k1) * S(1:k2)' .* hankel(S(2:k1+1),S(k1+1:k1+k2));
S1 = ones(k1,k2) ./ S1;
ind = find(Q > 0);
Q = Q(ind);
Bic = S1(ind) .* abs(B(ind)).^2; % squared bicoherence
%-------------------------------------------------------------------
%----------------- Gaussianity (non-skewness) test -----------------
%-------------------------------------------------------------------
scale = 2 * (Msmuth^4) / nfft;
Xstat = scale * Bic ./ Q; % 2 |X(m,n)|^2 in (2.9)
df = length(Q) * 2; % degrees of freedom
% asymptotic value is (nfft/M)^2 / 6
chi_val = sum(Xstat); % observed value of chi-square r.v.
% Reference:
% Sankaran's approximation given in [7.9.4], p 197, of
% J.K. Patel and C.B. Read, M. Dekker, New York, 1982.
% `Handbook of the Normal Distribution'
% the cdf of chi(df,lam,x) is approximated by phi(y) [normal distro cdf]
% where y = [ (x/(df+lam))^h - a ] / b;
% and h, a, b are defined below
% For a given df, and chi-sq value, the PFA is the probability that
% we will be wrong in accepting the alternate hypothesis; i.e.,
% data are non-skewed.
h = 1/3; b = sqrt(2/df) / 3; a = (1-b^2);
y = (chi_val/df)^(1/3); y = (y-a)/b;
pfa = 0.5 * erfc(y/sqrt(2));
pfa = round(pfa*10000)/10000; % round to five decimal places
disp(['Test statistic for Gaussianity is ', num2str(chi_val), ...
' with df = ',int2str(df), ', Pfa = ',num2str(pfa)])
sg = [chi_val, df, pfa];
% -----------------------------------------------------------------
% ---------------------- linearity test ---------------------------
% -----------------------------------------------------------------
% ---- find the first and third quartiles of 2|Xmn|^2 = Xstat
ind = find(Q == M^2);
rtest = Xstat(ind);
rtest = sort(rtest);
l1 = length(rtest);
if (l1 < 4)
fprintf('glstat: cparm (%g) is too large or nfft (%g) is too small', ...
cparm,nfft);
fprintf('\n estimated interquartile range (R) set to NaN \n');
Rhat = NaN;
else
lo1 = fix(l1/4); lo2 = fix(l1/4+0.8);
hi1 = fix(3*l1/4); hi2 = fix(3*l1/4+0.8);
r1 = (rtest(lo1)+rtest(lo2))/2;
r3 = (rtest(hi1)+rtest(hi2))/2;
Rhat = r3 - r1; % sample inter-quartile range
end
% ----- estimated lambda: a scaled version of mean FD skewness
% lam = 2/(df * Msmuth^2) * sum ( Q.*Xstat - 2); % Hinich's (4.2)
lam = mean(Bic) * 2 * Msmuth^2 / nfft; % equivalent
%----- find the `theoretical inter-quartile' range of chi-sq(df=2,lam)
% See reference cited earlier.
df = 2; % for our problem, df=2
h = 1/3 + (2/3) * (
MATLAB高阶累积量工具箱.rar
需积分: 46 86 浏览量
2015-07-09
18:53:00
上传
评论
收藏 101KB RAR 举报
Lxs_
- 粉丝: 22
- 资源: 2
最新资源
- 基于Python的图像阴影检测与去除源码(高分期末大作业项目).zip
- 基于C++/Qt实现的井字棋游戏
- 基于 Python 编程语言的 Web 框架Django
- Python和Flask实现的基于体检数据的城市公共健康可视分析系统源码+使用说明.zip
- 基于python实现的华为智慧工地-安全帽检测
- buck-boost_2023-12-16_12-12-13.eprj
- 后端开发关于数据库和API开发的介绍
- 机器学习和数据挖掘课程设计-米其林餐厅数据挖掘管理系统源码+使用文档说明.zip
- html html html展示我与ai的对化
- 数据结构课程设计-全国交通出行咨询模拟系统C语言实现源码.zip
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈