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

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_dsample [y,yy,yyy]=nt_dsample(x,dsr,method) - downsample by averaging neighboring samples
- nt_fold y=fold(x,epochsize) - fold 2D to 3D
- nt_unfold y=nt_fold(x) - unfold 3D to 2D
- nt_xprod [y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) - form all crossproducts

This function is called by:

- nt_bias_cluster [c0,c1,A,todss,pwr0,pwr1]=nt_bias_cluster(x,dsr,flags) - cluster covariance
- nt_cluster_jd [IDX,todss,SCORE,COVS]=nt_cluster_jd(x,dsr,smooth,flags,init,verbose) - cluster with joint diagonalization
- nt_xprod [y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) - form all crossproducts

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