0001 function [topcs,pwr,y]=nt_pca_big(x,N,aggregate,normflag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
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
0065 [topcs2,pwr,z]=nt_pca_big(x*topcs,N,aggregate,normflag);
0066 topcs=topcs*topcs2;
0067 else
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
0075 end
0076
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);
0084 pwr=pwr(1:N);
0085 else
0086 topcs(:,end:N)=0;
0087 pwr(end:N)=0;
0088 end
0089
0090
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
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
0114 if 0
0115
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