


Nonstationarity test - continuous process
Usage:
sigma=nonst_test(data,A,sumV,params)
Input:
Note units have to be consistent. See chronux.m for more information.
data (1d array in samples) -- required
A quadratic coefficient matrix - (Compute this separately since
the computation is time consuming - [A,sumV]=quadcof(N,NW,order). order
has to < 4NW.)
sumV sum of the quadratic inverse basis vectors
params: structure with fields tapers, pad, Fs, fpass, err, trialave
-optional
tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional. If not
specified, use [NW K]=[3 5]
pad (padding factor for the FFT) - optional (can take values -1,0,1,2...).
-1 corresponds to no padding, 0 corresponds to padding
to the next highest power of 2 etc.
e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
to 512 points, if pad=1, we pad to 1024 points etc.
Defaults to 0.
Fs (sampling frequency) - optional. Default 1.
Output:
sigma (nonstationarity index Thomson, 2000)

0001 function sigma = nonst_stat(data,A,sumV,params) 0002 0003 % Nonstationarity test - continuous process 0004 % 0005 % Usage: 0006 % 0007 % sigma=nonst_test(data,A,sumV,params) 0008 % Input: 0009 % Note units have to be consistent. See chronux.m for more information. 0010 % data (1d array in samples) -- required 0011 % A quadratic coefficient matrix - (Compute this separately since 0012 % the computation is time consuming - [A,sumV]=quadcof(N,NW,order). order 0013 % has to < 4NW.) 0014 % sumV sum of the quadratic inverse basis vectors 0015 % params: structure with fields tapers, pad, Fs, fpass, err, trialave 0016 % -optional 0017 % tapers (precalculated tapers from dpss, or in the form [NW K] e.g [3 5]) -- optional. If not 0018 % specified, use [NW K]=[3 5] 0019 % pad (padding factor for the FFT) - optional (can take values -1,0,1,2...). 0020 % -1 corresponds to no padding, 0 corresponds to padding 0021 % to the next highest power of 2 etc. 0022 % e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT 0023 % to 512 points, if pad=1, we pad to 1024 points etc. 0024 % Defaults to 0. 0025 % Fs (sampling frequency) - optional. Default 1. 0026 % Output: 0027 % sigma (nonstationarity index Thomson, 2000) 0028 0029 0030 if nargin < 1; error('Need data'); end; 0031 if nargin < 2; params=[]; end; 0032 0033 order = length(A); 0034 N = length(data); 0035 %nfft=max(2^(nextpow2(N)+pad),N); 0036 [tapers,pad,Fs]=getparams(params); 0037 tapers=dpsschk(tapers,N,Fs); % check tapers 0038 0039 alpha=zeros(1,order); 0040 for j=1:order 0041 alpha(j) = trace(squeeze(A(:,:,j))*squeeze(A(:,:,j))); 0042 end; 0043 0044 tmp=mtfftc(data,tapers,N,Fs); 0045 %tmp=mtfftc(data,tapers,nfft,Fs); 0046 sigma = zeros(length(data),1); 0047 % Pbar = sum(abs(tmp).^2,2)./sum(weights.^2,2); 0048 Pbar=mean(abs(tmp).^2,2); 0049 for ii=1:order 0050 a0=real(sum(tmp'.*(squeeze(A(:,:,ii))*tmp.')))'/alpha(ii); 0051 sigma=sigma+alpha(ii)*(a0./Pbar-sumV(ii)).^2; 0052 end; 0053