[y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) - form all crossproducts y: crossproducts ind: linear index of cross-products x: data (time*channels*trials) flag: 'lower','nodiag','full' [default: 'lower'] dsratio: ratio by which to downsample cross-product [default: 1] 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: linear index 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 [default: 1] 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 y(:,iProd)=nt_dsample(xx,dsratio); 0043 ind(iProd)=sub2ind([nchans,nchans],kk+iDiag,kk); 0044 iProd=iProd+1; 0045 end 0046 end 0047 if normrow_flag 0048 for iRow=1:size(y,1) 0049 y(iRow,:)=y(iRow,:)/(eps+sum(y(iRow,1:nchans))); 0050 end 0051 end 0052 if strcmp(flag, 'nodiag') 0053 y=y(:,nchans+1:end); 0054 end 0055 0056 % this could be optimized to save memory: 0057 if strcmp(flag,'full') 0058 y0=y; 0059 y=zeros(nsamples,nchans,nchans); 0060 0061 [I,J]=ind2sub(nchans,ind); 0062 0063 for k=1:numel(I) 0064 y(:,I(k),J(k))=y0(:,k); 0065 y(:,J(k),I(k))=y0(:,k); 0066 end 0067 0068 ind=[]; 0069 end 0070 0071 % switch order 0072 % case 'colwise' 0073 % for iRow=1:nchans 0074 % for iCol=1:iRow 0075 % xx=x(:,iCol).*x(:,iRow); 0076 % y(:,iProd)=nt_dsample(xx,dsratio); 0077 % ind(iProd)=sub2ind([size(x,2),size(x,2)],iRow,iCol); 0078 % iProd=iProd+1; 0079 % end 0080 % end 0081 % case 'diagwise' 0082 % for iDiag=start:nchans-1 0083 % for kk=1:(nchans-iDiag) 0084 % xx=x(:,kk+iDiag).*x(:,kk); 0085 % y(:,iProd)=nt_dsample(xx,dsratio); 0086 % ind(iProd)=sub2ind([size(x,2),size(x,2)],kk+iDiag,kk); 0087 % iProd=iProd+1; 0088 % end 0089 % end 0090 % otherwise 0091 % error('unexpected order flag'); 0092 % end 0093 end 0094 0095