Home > NoiseTools > nt_xprod.m

nt_xprod

PURPOSE ^

[y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) - form all crossproducts

SYNOPSIS ^

function [y,ind]=nt_xprod(x,flag,dsratio,normrow_flag)

DESCRIPTION ^

[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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 28-Nov-2016 20:12:47 by m2html © 2005