nt_xprod

[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 rem(size(x,1),dsratio)
0028     warning('clipping x to multiple of dsratio');
0029     x=x(1:dsratio*floor(size(x,1)/dsratio), :,:);
0030 end
0031 
0032 if ndims(x)==3
0033     y=nt_fold(nt_xprod(nt_unfold(x),flag,dsratio), size(x,1)/dsratio);
0034 else
0035     [nsamples,nchans]=size(x);
0036     nsamples=floor(nsamples/dsratio);
0037 
0038     y=zeros(nsamples,nchans*(nchans+1)/2);
0039     ind=zeros(nchans*(nchans+1)/2,1);
0040     start=0;
0041     
0042     iProd=1;
0043     for iDiag=start:nchans-1
0044         for kk=1:(nchans-iDiag)
0045             xx=x(:,kk+iDiag).*x(:,kk);
0046             y(:,iProd)=nt_dsample(xx,dsratio);
0047             ind(iProd)=sub2ind([nchans,nchans],kk+iDiag,kk);
0048             iProd=iProd+1;
0049         end
0050     end
0051     if normrow_flag
0052         for iRow=1:size(y,1)
0053             y(iRow,:)=y(iRow,:)/(eps+sum(y(iRow,1:nchans)));
0054         end
0055     end
0056     if strcmp(flag, 'nodiag')
0057         y=y(:,nchans+1:end);
0058     end
0059 
0060     % this could be optimized to save memory:
0061     if strcmp(flag,'full')
0062         y0=y;
0063         y=zeros(nsamples,nchans,nchans);
0064 
0065         [I,J]=ind2sub(nchans,ind);
0066 
0067         for k=1:numel(I)
0068             y(:,I(k),J(k))=y0(:,k);
0069             y(:,J(k),I(k))=y0(:,k);
0070         end
0071 
0072         ind=[];
0073     end
0074     
0075 %     switch order
0076 %         case 'colwise'
0077 %             for iRow=1:nchans
0078 %                 for iCol=1:iRow
0079 %                     xx=x(:,iCol).*x(:,iRow);
0080 %                     y(:,iProd)=nt_dsample(xx,dsratio);
0081 %                     ind(iProd)=sub2ind([size(x,2),size(x,2)],iRow,iCol);
0082 %                     iProd=iProd+1;
0083 %                 end
0084 %             end
0085 %         case 'diagwise'
0086 %             for iDiag=start:nchans-1
0087 %                 for kk=1:(nchans-iDiag)
0088 %                     xx=x(:,kk+iDiag).*x(:,kk);
0089 %                     y(:,iProd)=nt_dsample(xx,dsratio);
0090 %                     ind(iProd)=sub2ind([size(x,2),size(x,2)],kk+iDiag,kk);
0091 %                     iProd=iProd+1;
0092 %                 end
0093 %             end
0094 %         otherwise
0095 %             error('unexpected order flag');
0096 %     end
0097 end
0098 
0099