Home > NoiseTools > nt_xxcorr.m

# nt_xxcorr

## PURPOSE

[C,idx]=nt_xxcorr(A,B,MAXLAG) - true unbiased cross-correlation

## SYNOPSIS

function [C,idx]=nt_xxcorr(A,B,MAXLAG)

## DESCRIPTION

    [C,idx]=nt_xxcorr(A,B,MAXLAG) - true unbiased cross-correlation

      C: unbiased cross-correlation function
      idx: index of largest extremum.

      A: first column vector
      B: second column vector
      MAXLAG: lags are between -MAXLAG and +MAXLAG.

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_xxcorr [C,idx]=nt_xxcorr(A,B,MAXLAG) - true unbiased cross-correlation

This function is called by:

- nt_xxcorr [C,idx]=nt_xxcorr(A,B,MAXLAG) - true unbiased cross-correlation

## SOURCE CODE

    0001 function [C,idx]=nt_xxcorr(A,B,MAXLAG)
    0002 %[C,idx]=nt_xxcorr(A,B,MAXLAG) - true unbiased cross-correlation
    0003 %
    0004 %  C: unbiased cross-correlation function
    0005 %  idx: index of largest extremum.
    0006 %
    0007 %  A: first column vector
    0008 %  B: second column vector
    0009 %  MAXLAG: lags are between -MAXLAG and +MAXLAG.
    0010 %
    0011
    0012 if nargin<1; error('!'); end
    0013 if nargin==1; B=A; MAXLAG=[]; end
    0014 if nargin==2;
    0015     if numel(B)==1;
    0016         MAXLAG=B; B=A;
    0017     else
    0018         MAXLAG=[];
    0019     end
    0020 end
    0021 if isempty(MAXLAG); MAXLAG=floor(size(B,1)/4); end
    0022 if size(A,1)==1; A=A(:); B=B(:); end
    0023 if size(B,1)<=2*MAXLAG; error('!'); end
    0024
    0025
    0026
    0027 if size(A,2)==1 && size(B,2)==1;
    0028
    0029     % single channels
    0030     B=B(MAXLAG:end-MAXLAG);
    0031     C=xcorr(A,B);
    0032     C=C(size(A,1)-1+(0:2*MAXLAG));
    0033     [~,idx]=max(abs(C));
    0034
    0035     % power for normalization
    0036     a=cumsum(A.^2);
    0037     a=a(size(B,1)-1:end)-[0;a(1:end-size(B,1)+1)];
    0038     b=sum(B.^2);
    0039     d=sqrt(a*b);
    0040
    0041     C=C./d; C(find(isnan(C)))=0;
    0042
    0043     if nargout==0;
    0044         abscissa=-MAXLAG:MAXLAG;
    0045         plot(abscissa,zeros(size(C)), 'k'); hold on
    0046         plot(abscissa,C);
    0047         plot(abscissa(idx),C(idx),'.k'); hold off
    0048         axis tight; xlabel('lag (samples)');
    0049         lag=idx-MAXLAG-1;
    0050         if C(idx)>1
    0051             if lag>0;
    0052                 title(['X lags Y by ',num2str(lag)]);
    0053             else
    0054                 title(['Y lags X by ',num2str(-lag)]);
    0055             end
    0056         else
    0057             if lag>0;
    0058                 title(['X lags -Y by ',num2str(lag)]);
    0059             else
    0060                 title(['Y lags -X by ',num2str(-lag)]);
    0061             end
    0062         end
    0063         C=[];
    0064     end
    0065
    0066 else
    0067
    0068     % multiple channels
    0069     C=zeros(2*MAXLAG+1,size(A,2),size(B,2));
    0070     idx=zeros(size(A,2),size(B,2));
    0071     for k=1:size(A,2)
    0072         for j=1:size(B,2)
    0073             [a,b]=nt_xxcorr(A(:,k),B(:,j),MAXLAG);
    0074             C(:,k,j)=a;
    0075             idx(k,j)=b;
    0076         end
    0077     end
    0078
    0079     if nargout==0
    0080        imagescc(idx-MAXLAG-1);
    0081        colorbar
    0082        C=[]; idx=[];
    0083     end
    0084
    0085 end