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: This function is called by:

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

Generated on Mon 14-May-2018 09:45:18 by m2html © 2005