0001 function [tocomps,ii]=nt_peaky(c,x,T,nSmooth)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 nt_greetings;
0021
0022 if nargin<1; error('!'); end
0023
0024 topcs=[];
0025 if isempty(c);
0026
0027 if nargin<3; error('!'); end
0028 if nargin<4; nSmooth=1; end
0029 topcs=nt_pca0(x);
0030 x=nt_mmat(x,topcs);
0031 c=nt_xprod(x,'full',T);
0032 if nSmooth>1;
0033 c=filter(ones(nSmooth,1),1,c);
0034 c=c(nSmooth:end,:,:);
0035 end
0036 end
0037
0038 if ndims(c)~=3 || size(c,2)~=size(c,3)
0039 error('c has unexpected shape');
0040 end
0041 nComp=size(c,2);
0042
0043
0044 We iterate from 1 to nComp, updating the tocomps matrix at each step.
0045 At each step we iterate over time intervals. At each time interval we
0046 apply DSS using interval / total power ratio as a bias. We select the
0047 interval that gives the greatest ratio.
0048
0049 We then recalculate the corresponding DSS solution, put the first component
0050 into the tocomps matrix, and recurse on the remainder.
0051
0052 Calculations are done on covariance matrices rather than actual data.
0053
0054
0055 c0=squeeze(mean(c));
0056 if ~isempty(topcs)
0057 tocomps=topcs;
0058 else
0059 tocomps=eye(nComp);
0060 end
0061 iComp=1;
0062 [m n o]= size(c);
0063 while iComp<nComp
0064
0065
0066 ratio=zeros(1,size(c,1));
0067 for iIter=1: size(c,1);
0068 if size(c0,1)>150
0069 N=10;
0070 [V,D]=eigs(squeeze(c(iIter,:,:)),c0,N);
0071 else
0072 [V,D]=eig(squeeze(c(iIter,:,:)),c0);
0073 end
0074 D=sort(diag(D),'descend');
0075 if 0
0076 ratio(iIter)=D(2)/D(1);
0077 else
0078 ratio(iIter)=mean(D(2:end))/D(1);
0079 end
0080 end
0081 [~,idx]=min(ratio);
0082 ii(iComp)=idx;
0083
0084
0085
0086
0087
0088
0089 [todss,pwr0,pwr1]=nt_dss0(c0,squeeze(c(idx,:,:)),[],10.^-10);
0090 profile=pwr1./pwr0;
0091 score(iComp)=profile(1)/profile(2);
0092
0093
0094 tocomps=[tocomps(:,1:iComp-1), tocomps(:,iComp:end)*todss];
0095 nComp=size(tocomps,2);
0096
0097
0098 if iComp<nComp
0099 todss=todss(:,2:end);
0100 c0=todss'*c0*todss;
0101
0102 For speed, we reshape c to apply left and right multiplication to all
0103 matrices at once.
0104
0105 [m n o]=size(c);
0106 c=reshape(permute(c,[2 3 1]), [n o*m]);
0107 c=todss'*c;
0108 n=size(c,1);
0109 c=reshape(permute(reshape(c,[n o m]),[1 3 2]),[m*n o]);
0110 c=c*todss;
0111 o=size(c,2);
0112 c=permute(reshape(c, [n m o]),[2 1 3]) ;
0113 end
0114 iComp=iComp+1;
0115 end
0116
0117
0118
0119