0001 function [todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)
0002
0003
0004
0005
0006
0007
0008
0009
0010 if nargin<1; error('!'); end
0011 if nargin<2||isempty(nIterations); nIterations=5; end
0012 if nargin<3||isempty(exponent); exponent=2; end
0013 if nargin<4; w=[]; end
0014 if nargin<5; smooth=[];end
0015
0016 if ndims(x)>2; x=nt_unfold(x); end
0017
0018 KK={};
0019 K0=kurtosis(x);
0020 if nIterations==1
0021
0022
0023 We choose the channel with highest kurtosis to define a mask to
0024 emphasize high amplitude intervals. DSS then finds components with
0025 variance maximal within those intervals.
0026
0027
0028 K=kurtosis(x);
0029 [~,idx]=max(K);
0030 mask=x(:,idx).^exponent;
0031 mask=filtfilt(ones(5,1),1,mask);
0032 if ~isempty(smooth)&&smooth>1;
0033 mask=filtfilt(ones(1,smooth),1,mask);
0034 end
0035 c0=nt_cov(x,[],w);
0036 c1=nt_cov(nt_demean(x.*bsxfun(@times,x,mask),w),[],w);
0037 [todss]=nt_dss0(c0,c1);
0038 KK{1}=kurtosis(nt_mmat(x,todss));
0039
0040 else
0041
0042
0043 Apply DSS repeatedly, concatenating the DSS matrices. We hope (no
0044 guarantee) that the solution will reveal components with high kurtosis.
0045
0046
0047 z=x;
0048 todss=eye(size(x,2));
0049 for k=1:nIterations
0050 [T]=nt_kurtosis(z,1,exponent,w,smooth);
0051 z=nt_mmat(z,T);
0052 todss=todss*T;
0053 KK{k}=kurtosis(z);
0054 end
0055
0056 end
0057 K=KK{end};
0058
0059 if nargout==0;
0060
0061 figure(101); clf
0062 subplot 121;
0063 plot(K0); hold on
0064 for k=1:nIterations
0065 plot(KK{k});
0066 end
0067 nt_colorlines;
0068 title('kurtosis score for each iteration'); xlabel('component');
0069 legend(num2str((0:nIterations)')); set(gca,'yscale','log')
0070 subplot 222;
0071 z=nt_mmat(x,todss(:,1));
0072 plot(z); title('first DSS component');
0073 subplot 224;
0074 z=nt_mmat(x,todss(:,end));
0075 plot(z); title('last DSS component');
0076 clear c0 c1 todss pwr0 pwr1 K
0077 end
0078