0001 function [tosquares,quads,D]=nt_qpca0(x,npcs,nsmooth,nquads)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 if nargin<4; nquads=1; end
0017 if nargin<3; nsmooth=1; end
0018 if nargin<2; npcs=[]; end
0019 [nsamples,nchans,ntrials]=size(x);
0020
0021
0022
0023 THRESH=10^-12;
0024 [topcs,pwr]=nt_pca0(x,[],[],THRESH);
0025 if isempty(npcs) || npcs>size(topcs,2); npcs=size(topcs,2); end
0026 topcs=topcs*diag(1./sqrt(pwr));
0027 topcs=topcs(:,1:npcs);
0028 x=nt_mmat(x,topcs);
0029
0030 nquads=min(nquads,npcs*(npcs+1)/2-1);
0031
0032
0033 c0=zeros(npcs*(npcs+1)/2);
0034 xxx=zeros(nsamples-nsmooth+1,npcs*(npcs+1)/2);
0035 for k=1:ntrials
0036 xx=zeros(nsamples,npcs*(npcs+1)/2);
0037 ii=1;
0038 for jj=1:npcs
0039 for kk=1:jj
0040 xx(:,ii)=x(:,kk,k).*x(:,jj,k);
0041 ii=ii+1;
0042 end
0043 end
0044 xx=filter(ones(nsmooth,1)/nsmooth,1,xx);
0045 xx=xx(nsmooth:end,:,:);
0046 xxx=xxx+xx;
0047 c0=c0+xx'*xx;
0048 end
0049
0050
0051 topcs2=nt_pcarot(c0);
0052 tobest=topcs2(:,2);
0053
0054
0055 ii=1;
0056 A=zeros(npcs);
0057 for k=1:npcs
0058 for j=1:k
0059 if j==k;
0060 A(k,j)=tobest(ii);
0061 else
0062 A(k,j)=tobest(ii)/2;
0063 A(j,k)=tobest(ii)/2;
0064 end
0065 ii=ii+1;
0066 end
0067 end
0068
0069
0070 [V,D]=eig(A);
0071 D=diag(D);
0072 [dummy,idx]=sort(abs(D),'descend');
0073 V=V(:,idx);
0074 D=D(idx);
0075 tosquares=topcs*V;
0076
0077
0078 if nargout>1;
0079 quads=zeros(nsamples-nsmooth+1,nquads+1,ntrials);
0080 for k=1:ntrials
0081 xx=zeros(nsamples,npcs*(npcs+1)/2);
0082 ii=1;
0083 for jj=1:npcs
0084 for kk=1:jj
0085 xx(:,ii)=x(:,kk,k).*x(:,jj,k);
0086 ii=ii+1;
0087 end
0088 end
0089 xx=filter(ones(nsmooth,1)/nsmooth,1,xx);
0090 xx=xx(nsmooth:end,:,:);
0091 quads(:,:,k)=nt_mmat(xx,topcs2(:,1:nquads+1));
0092 end
0093 end
0094