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: This function is called by:

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

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