Home > NoiseTools > nt_qca02.m

# nt_qca02

## PURPOSE

[tosquare,quad,tosquare2,quad2,D]=nt_qca02(x,npcs,nsmooth) - maximize induced power using quadratic component analysis

## SYNOPSIS

function [tosquare,quad,tosquare2,quad2,D]=nt_qca02(x,npcs,nsmooth)

## DESCRIPTION

    [tosquare,quad,tosquare2,quad2,D]=nt_qca02(x,npcs,nsmooth) - maximize induced power using quadratic component analysis

      tosquare: matrix to most reproducible induced component
      quad: most reproducible quadratic component
      tosquare2: matrix to most reproducible induced component
      quad2: most reproducible quadratic component
      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_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- 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

This function is called by:

- nt_qca2 [square,quad,square2,quad2,D]=nt_qca(x,npcs,nsmooth) - maximize induced power using quadratic component analysis

## SOURCE CODE

    0001 function [tosquare,quad,tosquare2,quad2,D]=nt_qca02(x,npcs,nsmooth)
    0002 %[tosquare,quad,tosquare2,quad2,D]=nt_qca02(x,npcs,nsmooth) - maximize induced power using quadratic component analysis
    0003 %
    0004 %  tosquare: matrix to most reproducible induced component
    0005 %  quad: most reproducible quadratic component
    0006 %  tosquare2: matrix to most reproducible induced component
    0007 %  quad2: most reproducible quadratic component
    0008 %  D: eigenvalues sorted by absolute value
    0009 %
    0010 %  x: data (time*channel*trial)
    0011 %  npcs: maximum number of data PCs to use [default: all]
    0012 %  nsmooth: square smoothing window to apply to xproducts [default: 1]
    0013 %  nquads: number of quadratic components to return [default: 1]
    0014 %
    0015 % NoiseTools.
    0016
    0017
    0018 if nargin<3; nsmooth=1; end
    0019 if nargin<2; npcs=[]; end
    0020 [nsamples,nchans,ntrials]=size(x);
    0021
    0022
    0023 % PCA & normalize, select PCs to save space
    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 % covariance of cross-products
    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);                      % lowpass (demodulate)
    0044     xx=xx(nsmooth:end,:,:);
    0045     xxx=xxx+xx;
    0046     c0=c0+xx'*xx;
    0047 end
    0048 c1=xxx'*xxx;
    0049
    0050 % DSS to find most repeatable
    0051 [todss,pwr0,pwr1]=nt_dss0(c0,c1,[],0);
    0052 tobest=todss(:,2); % first is DC
    0053 toworst=todss(:,end);
    0054
    0055 % form square matrix
    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 % eigenvectors & values
    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 % form square matrix
    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 % eigenvectors & values
    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 % best & worst quadratic component(s)
    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);                      % lowpass (demodulate)
    0117         xx=xx(nsmooth:end,:,:);                                       % chop off onset artifact
    0118         quad(:,:,k)=nt_mmat(xx,todss(:,1));
    0119         quad2(:,:,k)=nt_mmat(xx,todss(:,end));
    0120     end
    0121 end
    0122