0001 function [topcs,pwr,y]=nt_pca_big(x,N,chunksize,normflag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
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
0041 else
0042 [topcs,pwr]=nt_pca0(x);
0043 end
0044 if size(topcs,2)>=N;
0045 topcs=topcs(:,1:N);
0046 pwr=pwr(1:N);
0047 else
0048 topcs(:,end:N)=0;
0049 pwr(end:N)=0;
0050 end
0051
0052
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
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
0076 if 0
0077
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