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:

- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_pca0 [topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_pcarot [topcs,eigenvalues]=pcarot(cov,nkeep,threshold,N) - PCA matrix from covariance

This function is called by:

- nt_qpca [squares,quads]=nt_qpca(x,npcs,nsmooth,nquads) - quadratic PCA

## 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