0001 function [tosquare,quad,tosquare2,quad2,D]=nt_qca02(x,npcs,nsmooth)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 if nargin<3; nsmooth=1; end
0019 if nargin<2; npcs=[]; end
0020 [nsamples,nchans,ntrials]=size(x);
0021
0022
0023
0024 THRESH=10^-12;
0025 [topcs,pwr]=nt_pca0(x,[],[],THRESH);
0026 if isempty(npcs) || npcs>size(topcs,2); npcs=size(topcs,2); end
0027 topcs=topcs*diag(1/sqrt(pwr));
0028 topcs=topcs(:,1:npcs);
0029 x=nt_mmat(x,topcs);
0030
0031
0032 c0=zeros(npcs*(npcs+1)/2);
0033 xxx=zeros(nsamples-nsmooth+1,npcs*(npcs+1)/2);
0034 for k=1:ntrials
0035 xx=zeros(nsamples,npcs*(npcs+1)/2);
0036 ii=1;
0037 for jj=1:npcs
0038 for kk=1:jj
0039 xx(:,ii)=x(:,kk,k).*x(:,jj,k);
0040 ii=ii+1;
0041 end
0042 end
0043 xx=filter(ones(nsmooth,1)/nsmooth,1,xx);
0044 xx=xx(nsmooth:end,:,:);
0045 xxx=xxx+xx;
0046 c0=c0+xx'*xx;
0047 end
0048 c1=xxx'*xxx;
0049
0050
0051 [todss,pwr0,pwr1]=nt_dss0(c0,c1,[],0);
0052 tobest=todss(:,2);
0053 toworst=todss(:,end);
0054
0055
0056 ii=1;
0057 A=zeros(npcs);
0058 for k=1:npcs
0059 for j=1:k
0060 if j==k;
0061 A(k,j)=tobest(ii);
0062 else
0063 A(k,j)=tobest(ii)/2;
0064 A(j,k)=tobest(ii)/2;
0065 end
0066 ii=ii+1;
0067 end
0068 end
0069
0070
0071 [V,D]=eig(A);
0072 D=diag(D);
0073 [dummy,idx]=sort(abs(D),'descend');
0074 V=V(:,idx);
0075 D=D(idx);
0076 tosquare=topcs*V;
0077
0078
0079 ii=1;
0080 A=zeros(npcs);
0081 for k=1:npcs
0082 for j=1:k
0083 if j==k;
0084 A(k,j)=toworst(ii);
0085 else
0086 A(k,j)=toworst(ii)/2;
0087 A(j,k)=toworst(ii)/2;
0088 end
0089 ii=ii+1;
0090 end
0091 end
0092
0093
0094 [V,D]=eig(A);
0095 D=diag(D);
0096 [dummy,idx]=sort(abs(D),'descend');
0097 V=V(:,idx);
0098 D=D(idx);
0099 tosquare2=topcs*V;
0100
0101
0102
0103
0104 if nargout>1;
0105 quads=zeros(nsamples-nsmooth+1,1,ntrials);
0106 quads2=zeros(nsamples-nsmooth+1,1,ntrials);
0107 for k=1:ntrials
0108 xx=zeros(nsamples,npcs*(npcs+1)/2);
0109 ii=1;
0110 for jj=1:npcs
0111 for kk=1:jj
0112 xx(:,ii)=x(:,kk,k).*x(:,jj,k);
0113 ii=ii+1;
0114 end
0115 end
0116 xx=filter(ones(nsmooth,1)/nsmooth,1,xx);
0117 xx=xx(nsmooth:end,:,:);
0118 quad(:,:,k)=nt_mmat(xx,todss(:,1));
0119 quad2(:,:,k)=nt_mmat(xx,todss(:,end));
0120 end
0121 end
0122