


Multi-taper spectrum - 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. Defaults to 0.
e.g. For N = 500, if PAD = 0, we pad the FFT
to 512 points; if PAD = 2, we pad the FFT
to 2048 points, etc.
Fs (sampling frequency) - optional. Default 1.
Output:
sigma (nonstationarity index Thomson, 2000)

0001 function sigma = nonst_stat(data,A,sumV,params) 0002 0003 % Multi-taper spectrum - 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. Defaults to 0. 0020 % e.g. For N = 500, if PAD = 0, we pad the FFT 0021 % to 512 points; if PAD = 2, we pad the FFT 0022 % to 2048 points, etc. 0023 % Fs (sampling frequency) - optional. Default 1. 0024 % Output: 0025 % sigma (nonstationarity index Thomson, 2000) 0026 0027 0028 if nargin < 1; error('Need data'); end; 0029 if nargin < 2; params=[]; end; 0030 0031 order = length(A); 0032 N = length(data); 0033 %nfft=2^(nextpow2(N)+pad); 0034 [tapers,pad,Fs]=getparams(params); 0035 tapers=dpsschk(tapers,N,Fs); % check tapers 0036 0037 alpha=zeros(1,order); 0038 for j=1:order 0039 alpha(j) = trace(squeeze(A(:,:,j))*squeeze(A(:,:,j))); 0040 end; 0041 0042 tmp=mtfftc(data,tapers,N,Fs); 0043 %tmp=mtfftc(data,tapers,nfft,Fs); 0044 sigma = zeros(length(data),1); 0045 % Pbar = sum(abs(tmp).^2,2)./sum(weights.^2,2); 0046 Pbar=mean(abs(tmp).^2,2); 0047 for ii=1:order 0048 a0=real(sum(tmp'.*(squeeze(A(:,:,ii))*tmp.')))'/alpha(ii); 0049 sigma=sigma+alpha(ii)*(a0./Pbar-sumV(ii)).^2; 0050 end; 0051