Home > NoiseTools > nt_pca_big.m

nt_pca_big

PURPOSE ^

[topcs,pwr,y]=nt_pca_big(x,N,chunksize,normflag) - approximate first PCs of big data

SYNOPSIS ^

function [topcs,pwr,y]=nt_pca_big(x,N,chunksize,normflag)

DESCRIPTION ^

[topcs,pwr,y]=nt_pca_big(x,N,chunksize,normflag) - approximate first PCs of big data

  r: matrix to project to PCs
  pwr: power per PC
  y: PCs

  x: data
  N: desired number of components [default: 10]
  chunksize: number of channels per chunk
  normflag: normalize at each level [default: don't]

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [topcs,pwr,y]=nt_pca_big(x,N,chunksize,normflag)
0002 %[topcs,pwr,y]=nt_pca_big(x,N,chunksize,normflag) - approximate first PCs of big data
0003 %
0004 %  r: matrix to project to PCs
0005 %  pwr: power per PC
0006 %  y: PCs
0007 %
0008 %  x: data
0009 %  N: desired number of components [default: 10]
0010 %  chunksize: number of channels per chunk
0011 %  normflag: normalize at each level [default: don't]
0012 %
0013 
0014 
0015 if nargin<1; error('!'); end
0016 if nargin<2||isempty(N); 
0017     N=10;
0018     disp('defaulting to N=10');
0019 end
0020 if nargin<3||isempty(chunksize); chunksize=500; end
0021 if nargin<4||isempty(normflag); normflag=0; end
0022 if chunksize<2; error('!'); end
0023 
0024 [nsamples,nchans,ntrials]=size(x); 
0025 x=nt_unfold(x);
0026 
0027 if normflag;
0028     tonorm=1./sqrt(mean(x.^2)); 
0029     tonorm(find(isnan(tonorm)))=0;
0030     x=bsxfun(@times,x,tonorm);
0031 end
0032 
0033 if nchans>chunksize; 
0034     topcs=zeros(nchans,N*2);
0035     a=round(nchans/2);
0036     [topcs(1:a,1:N),p,z]=nt_pca_big(x(:,1:a),N,chunksize,normflag);
0037     [topcs(a+1:end,N+1:2*N),p,z]=nt_pca_big(x(:,a+1:end),N,chunksize,normflag);
0038     [topcs2,pwr,z]=nt_pca_big(x*topcs,N,chunksize,normflag);
0039     topcs=topcs*topcs2;
0040     %figure(1); plot(xx*tp(:,1)); drawnow; pause
0041 else
0042     [topcs,pwr]=nt_pca0(x);
0043 end
0044 if size(topcs,2)>=N; 
0045     topcs=topcs(:,1:N); % truncate to save space
0046     pwr=pwr(1:N);
0047 else
0048     topcs(:,end:N)=0; % should rarely happen
0049     pwr(end:N)=0;
0050 end
0051    
0052 %disp([size(topcs)])
0053 
0054 if normflag; topcs=bsxfun(@times,topcs,tonorm'); end
0055 
0056 if nargout>1; y=x*topcs; y=nt_fold(y,nsamples);end
0057 
0058 if ~nargout % don't return, just plot something
0059     y=nt_unfold(y);
0060     semilogy(mean(y.^2), '.-');
0061     xlabel('pc'); ylabel('power');
0062     figure(2); clf
0063     plot(y); xlabel('sample');
0064     figure(3); clf
0065     nt_imagescc(nt_normcol(x'*y)'); 
0066     ylabel('pc'); xlabel('channel');
0067     clear topcs y
0068 end
0069 
0070 
0071 
0072 
0073 
0074 
0075 %% test code
0076 if 0
0077     % basic
0078     x=randn(10000,300);
0079     [topcs,pwr,y]=nt_pca_big(x);
0080     figure(1); clf; plot(pwr);
0081 end
0082 if 0 
0083     x=sin(2*pi*3*(1:1000)'/1000)*randn(1,1000);
0084     x=2*x + randn(size(x));
0085      [topcs,pwr,y]=nt_pca_big(x);
0086     figure(1); plot(pwr);
0087     figure(2); subplot 121; plot(x); subplot 122; plot(x*topcs);
0088 end   
0089 if 0 
0090     nchans=100000;
0091     nnoise=900;
0092     x=repmat(sin(2*pi*3*(1:1000)'/1000), 1, nchans);
0093     x= 0.06*nt_normcol(x) + nt_normcol(randn(size(x,1),nnoise)*randn(nnoise,nchans));
0094     [topcs,pwr,y]=nt_pca_big(x,10);
0095     figure(1); plot(pwr);
0096     figure(2); subplot 121; plot(x(:,1:100)); subplot 122; plot(y(:,1));
0097 end   
0098

Generated on Mon 28-Nov-2016 20:12:47 by m2html © 2005