Home > NoiseTools > nt_qpca0.m

nt_qpca0

PURPOSE ^

[tosquares,quads,D]=nt_qpca0(x,npcs,nsmooth,nquads) - quadratic PCA

SYNOPSIS ^

function [tosquares,quads,D]=nt_qpca0(x,npcs,nsmooth,nquads)

DESCRIPTION ^

[tosquares,quads,D]=nt_qpca0(x,npcs,nsmooth,nquads) - quadratic PCA

  tosquares: matrix to linear components closest to largest quadradric component
  quads: largest quadratic component(s)
  D: eigenvalues sorted by absolute value

  x: data (time*channel*trial)
  npcs: maximum number of data PCs to use [default: all]
  nsmooth: square smoothing window to apply to xproducts [default: 1]
  nquads: number of quadratic components to return [default: 1]
 
 NoiseTools.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [tosquares,quads,D]=nt_qpca0(x,npcs,nsmooth,nquads)
0002 %[tosquares,quads,D]=nt_qpca0(x,npcs,nsmooth,nquads) - quadratic PCA
0003 %
0004 %  tosquares: matrix to linear components closest to largest quadradric component
0005 %  quads: largest quadratic component(s)
0006 %  D: eigenvalues sorted by absolute value
0007 %
0008 %  x: data (time*channel*trial)
0009 %  npcs: maximum number of data PCs to use [default: all]
0010 %  nsmooth: square smoothing window to apply to xproducts [default: 1]
0011 %  nquads: number of quadratic components to return [default: 1]
0012 %
0013 % NoiseTools.
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 % PCA & normalize, select PCs to save space
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 % covariance of cross-products
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);                      % lowpass (demodulate)
0045     xx=xx(nsmooth:end,:,:);
0046     xxx=xxx+xx;
0047     c0=c0+xx'*xx;
0048 end
0049 
0050 % DSS to find most repeatable
0051 topcs2=nt_pcarot(c0);
0052 tobest=topcs2(:,2); % first is DC
0053 
0054 % form square matrix
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 % eigenvectors & values
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 % answer best quadratic component(s)
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);                      % lowpass (demodulate)
0090         xx=xx(nsmooth:end,:,:);                                       % chop off onset artifact
0091         quads(:,:,k)=nt_mmat(xx,topcs2(:,1:nquads+1));
0092     end
0093 end
0094

Generated on Sat 29-Apr-2023 17:15:46 by m2html © 2005