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
共57个文件
m:57个
需积分: 46 34 下载量 155 浏览量
2015-07-09
18:53:00
上传
评论
收藏 101KB RAR 举报
温馨提示
高阶累积量工具箱——添加后可直接使用。课程用过,现在分享给大家。我分也不多,所以少分分享出去啦~
资源推荐
资源详情
资源评论
收起资源包目录
MATLAB高阶累积量工具箱.rar (57个子文件)
MATLAB高阶累积量工具箱
hosa
armaqs.m 4KB
trench.m 2KB
wig2.m 3KB
tde.m 4KB
wig4.m 3KB
cumtrue.m 4KB
tdegen.m 4KB
rivdl.m 5KB
cum4x.m 6KB
tder.m 3KB
nlpow.m 4KB
nlgen.m 2KB
tls.m 1KB
cumest.m 3KB
glstat.m 8KB
hprony.m 3KB
cum2x.m 3KB
hosahelp.m 4KB
bispecdx.m 6KB
bicoherx.m 5KB
tdeb.m 3KB
bispecd.m 6KB
cum3est.m 3KB
ivcal.m 2KB
qpctor.m 5KB
arrcest.m 4KB
maorder.m 3KB
maest.m 3KB
rpiid.m 2KB
nltick.m 5KB
cum4est.m 4KB
wig2c.m 4KB
wig3.m 3KB
Contents.m 5KB
trispect.m 3KB
matul.m 3KB
bicoher.m 4KB
doa.m 6KB
bicepsf.m 5KB
bispect.m 2KB
biceps.m 5KB
armarts.m 3KB
doagen.m 5KB
wig4c.m 4KB
Readme.m 6KB
rivtr.m 3KB
cum2est.m 2KB
bispeci.m 5KB
armasyn.m 4KB
qpcgen.m 6KB
harmest.m 6KB
harmgen.m 4KB
cum3x.m 4KB
use_doa.m 534B
arorder.m 4KB
pickpeak.m 3KB
wig3c.m 4KB
共 57 条
- 1
资源评论
Lxs_
- 粉丝: 22
- 资源: 2
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功