[y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) - form all crossproducts y: crossproducts ind: indices of cross-products x: data (time*channels*trials) flag: 'lower','nodiag','full' [default: 'lower'] dsratio: ratio by which to downsample cross-product normrow_flag: if true, divide each slice by trace (default: false) If flag is 'lower' (default), return lower diagonal terms, including diagonal. The order of cross-products is diagonal first (squares first). If 'nodiag' return lower diagonal terms without diagonal. If 'full', return full array of cross-products as a 3D matrix.
0001 function [y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) 0002 %[y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) - form all crossproducts 0003 % 0004 % y: crossproducts 0005 % ind: indices of cross-products 0006 % 0007 % x: data (time*channels*trials) 0008 % flag: 'lower','nodiag','full' [default: 'lower'] 0009 % dsratio: ratio by which to downsample cross-product 0010 % normrow_flag: if true, divide each slice by trace (default: false) 0011 % 0012 % If flag is 'lower' (default), return lower diagonal terms, including 0013 % diagonal. The order of cross-products is diagonal first (squares first). 0014 % 0015 % If 'nodiag' return lower diagonal terms without diagonal. 0016 % 0017 % If 'full', return full array of cross-products as a 3D matrix. 0018 % 0019 % 0020 0021 if nargin<4 || isempty(normrow_flag); normrow_flag=0; end 0022 if nargin<3 || isempty(dsratio); dsratio=1; end 0023 if nargin<2 || isempty(flag); flag='lower'; end 0024 0025 if ~strcmp(flag,{'lower','nodiag','full'}); error('unexpected flag'); end 0026 0027 if ndims(x)==3 0028 if rem(size(x,1),dsratio); error('size(x,1) must be multiple of dsratio'); end 0029 y=nt_fold(nt_xprod(nt_unfold(x),flag,dsratio), size(x,1)/dsratio); 0030 else 0031 [nsamples,nchans]=size(x); 0032 nsamples=floor(nsamples/dsratio); 0033 0034 y=zeros(nsamples,nchans*(nchans+1)/2); 0035 ind=zeros(nchans*(nchans+1)/2,1); 0036 start=0; 0037 0038 iProd=1; 0039 for iDiag=start:nchans-1 0040 for kk=1:(nchans-iDiag) 0041 xx=x(:,kk+iDiag).*x(:,kk); 0042 % if start==0; 0043 % xx=xx/2; % divide diagonal terms by 2 0044 % end 0045 y(:,iProd)=nt_dsample(xx,dsratio); 0046 ind(iProd)=sub2ind([size(x,2),size(x,2)],kk+iDiag,kk); 0047 iProd=iProd+1; 0048 end 0049 end 0050 if normrow_flag 0051 for iRow=1:size(y,1) 0052 y(iRow,:)=y(iRow,:)/(eps+sum(y(iRow,1:nchans))); 0053 end 0054 end 0055 if strcmp(flag, 'nodiag') 0056 y=y(:,nchans+1:end); 0057 end 0058 0059 % this could be optimized to save memory: 0060 if strcmp(flag,'full') 0061 y0=y; 0062 y=zeros(nsamples,nchans,nchans); 0063 0064 [I,J]=ind2sub(nchans,ind); 0065 0066 for k=1:numel(I) 0067 y(:,I(k),J(k))=y0(:,k); 0068 y(:,J(k),I(k))=y0(:,k); 0069 end 0070 0071 ind=[]; 0072 end 0073 0074 % switch order 0075 % case 'colwise' 0076 % for iRow=1:nchans 0077 % for iCol=1:iRow 0078 % xx=x(:,iCol).*x(:,iRow); 0079 % y(:,iProd)=nt_dsample(xx,dsratio); 0080 % ind(iProd)=sub2ind([size(x,2),size(x,2)],iRow,iCol); 0081 % iProd=iProd+1; 0082 % end 0083 % end 0084 % case 'diagwise' 0085 % for iDiag=start:nchans-1 0086 % for kk=1:(nchans-iDiag) 0087 % xx=x(:,kk+iDiag).*x(:,kk); 0088 % y(:,iProd)=nt_dsample(xx,dsratio); 0089 % ind(iProd)=sub2ind([size(x,2),size(x,2)],kk+iDiag,kk); 0090 % iProd=iProd+1; 0091 % end 0092 % end 0093 % otherwise 0094 % error('unexpected order flag'); 0095 % end 0096 end 0097 0098