Home > NoiseTools > nt_pca_big_old.m

nt_pca_big_old

PURPOSE ^

[topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag) - approximate PCA of big data

SYNOPSIS ^

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

DESCRIPTION ^

[topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag) - approximate PCA of big data

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

  x: data
  N: desired number of PCs
  aggregate: 'local' [default], 'global', or 'globmean'
  normflag: normalize at each level [default: don't]

 Data are processed recursively in chunks of 2*N columns.  Each chunk is
 submitted to PCA and either N PCs are retained ('local') or 1 PC is
 retained ('global').  If 'globmean' the data are averaged over chunks 
 of nchans/N channels, then PCA

 NoiseTools

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag)
0002 %[topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag) - approximate PCA 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 PCs
0010 %  aggregate: 'local' [default], 'global', or 'globmean'
0011 %  normflag: normalize at each level [default: don't]
0012 %
0013 % Data are processed recursively in chunks of 2*N columns.  Each chunk is
0014 % submitted to PCA and either N PCs are retained ('local') or 1 PC is
0015 % retained ('global').  If 'globmean' the data are averaged over chunks
0016 % of nchans/N channels, then PCA
0017 %
0018 % NoiseTools
0019 nt_greetings();
0020 
0021 if nargin<1; error('!'); end
0022 if nargin<2||isempty(N); 
0023     N=100;
0024     disp('defaulting to N=30');
0025 end
0026 if nargin<3||isempty(aggregate); aggregate='local'; end
0027 if nargin<4||isempty(normflag); normflag=0; end
0028 
0029 [nsamples,nchans,ntrials]=size(x); 
0030 x=nt_unfold(x);
0031 
0032 if normflag;
0033     tonorm=1./sqrt(mean(x.^2)); 
0034     tonorm(find(isnan(tonorm)))=0;
0035     x=bsxfun(@times,x,tonorm);
0036 end
0037 if ~ (strcmp(aggregate,'local') || strcmp(aggregate,'global') || strcmp(aggregate,'globmean'))
0038     error('!');
0039 end
0040 
0041 
0042 
0043 if nchans<=2*N; 
0044     [topcs,pwr]=nt_pca0(x);
0045 else
0046     if strcmp(aggregate,'local')
0047         topcs=zeros(nchans,N*2);
0048         a=round(nchans/2);
0049         [topcs(1:a,1:N)]=nt_pca_big(x(:,1:a),N,aggregate,normflag);
0050         [topcs(a+1:end,N+1:2*N)]=nt_pca_big(x(:,a+1:end),N,aggregate,normflag);
0051         [topcs2,pwr,z]=nt_pca_big(x*topcs,N,aggregate,normflag);
0052         topcs=topcs*topcs2;
0053     elseif strcmp(aggregate,'global') 
0054         nChunks=ceil(nchans/(2*N));
0055         chunkSize=ceil(nchans/nChunks);
0056         nKeep=ceil(N/nChunks);
0057         topcs=zeros(nchans,nChunks*nKeep);
0058         for iChunk=1:nChunks
0059             a=(iChunk-1)*chunkSize;
0060             idx=(a+1):min(a+chunkSize,nchans);
0061             tp=nt_pca_big(x(:,idx),N,aggregate,normflag);
0062             topcs(idx,(iChunk-1)*nKeep+(1:nKeep))=tp(:,1:nKeep);
0063         end
0064         %figure(1); clf; nt_imagescc(topcs);
0065         [topcs2,pwr,z]=nt_pca_big(x*topcs,N,aggregate,normflag);
0066         topcs=topcs*topcs2;
0067     else % globmean
0068         chunkSize=ceil(nchans/(2*N));
0069         topcs=zeros(nchans,2*N);
0070         for iChunk=1:2*N
0071             a=(iChunk-1)*chunkSize;
0072             idx=(a+1):min(a+chunkSize,nchans);
0073             topcs(idx,iChunk)=1/numel(idx);
0074             %disp(idx)
0075         end
0076         %figure(1); clf; nt_imagescc(topcs);
0077         [topcs2,pwr,z]=nt_pca_big(x*topcs,N,aggregate,normflag);
0078         topcs=topcs*topcs2;
0079     end
0080         
0081 end
0082 if size(topcs,2)>=N; 
0083     topcs=topcs(:,1:N); % truncate to save space
0084     pwr=pwr(1:N);
0085 else
0086     topcs(:,end:N)=0; % should rarely happen
0087     pwr(end:N)=0;
0088 end
0089    
0090 %disp([size(topcs)])
0091 
0092 if normflag; topcs=bsxfun(@times,topcs,tonorm'); end
0093 
0094 if nargout>1; y=x*topcs; y=nt_fold(y,nsamples);end
0095 
0096 if ~nargout % don't return, just plot something
0097     y=nt_unfold(y);
0098     semilogy(mean(y.^2), '.-');
0099     xlabel('pc'); ylabel('power');
0100     figure(2); clf
0101     plot(y); xlabel('sample');
0102     figure(3); clf
0103     nt_imagescc(nt_normcol(x'*y)'); 
0104     ylabel('pc'); xlabel('channel');
0105     clear topcs y
0106 end
0107 
0108 
0109 
0110 
0111 
0112 
0113 %% test code
0114 if 0
0115     % basic
0116     x=randn(10000,3000);
0117     [topcs,pwr,y]=nt_pca_big(x,10,'local');
0118     figure(1); clf; plot(pwr);
0119 end
0120 if 0 
0121     x=sin(2*pi*3*(1:1000)'/1000)*randn(1,1000);
0122     x=2*x + randn(size(x));
0123      [topcs,pwr,y]=nt_pca_big(x);
0124     figure(1); plot(pwr);
0125     figure(2); subplot 121; plot(x); subplot 122; plot(x*topcs);
0126 end   
0127 if 0 
0128     nchans=100000;
0129     nnoise=900;
0130     x=repmat(sin(2*pi*3*(1:1000)'/1000), 1, nchans);
0131     x= 0.06*nt_normcol(x) + nt_normcol(randn(size(x,1),nnoise)*randn(nnoise,nchans));
0132     [topcs,pwr,y]=nt_pca_big(x,100,'global');
0133     figure(1); plot(pwr);
0134     figure(2); subplot 121; plot(x(:,1:100)); subplot 122; plot(y(:,1));
0135 end   
0136 if 0 
0137     nchans=100000;
0138     nnoise=900;
0139     x=repmat(sin(2*pi*3*(1:1000)'/1000), 1, nchans);
0140     x= 0.01*nt_normcol(x) + nt_normcol(randn(size(x,1),nnoise)*randn(nnoise,nchans));
0141     [topcs,pwr,y]=nt_pca_big(x,100,'globmean');
0142     figure(1); plot(pwr);
0143     figure(2); subplot 121; plot(x(:,1:100)); subplot 122; plot(y(:,1));
0144 end   
0145

Generated on Fri 02-Dec-2016 14:27:04 by m2html © 2005