nt_xxcorr

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

  C: normalized unbiased cross-correlation function
  idx: index of largest extremum.
  
  A: first column vector
  B: second column vector (nsamples < A)
  centerflag: if true pad B on both sides [default: false]

 NoiseTools

0001 function [C,idx]=nt_xxcorr(A,B,centerflag)
0002 %[C,idx]=nt_xxcorr(A,B,centerflag) - true normalized unbiased cross-correlation function
0003 %
0004 %  C: normalized unbiased cross-correlation function
0005 %  idx: index of largest extremum.
0006 %
0007 %  A: first column vector
0008 %  B: second column vector (nsamples < A)
0009 %  centerflag: if true pad B on both sides [default: false]
0010 %
0011 % NoiseTools
0012 
0013 if nargin<2; error('!'); end
0014 if nargin<3||isempty(centerflag); centerflag=0; end
0015 if size(A,1)==size(B,1); warning('same nrows, output=1'); end
0016 if size(A,1)<size(B,1); error('should be nrows(A)>nrows(B)'); end
0017 
0018 
0019 if size(A,2)==1 && size(B,2)==1; 
0020     
0021     % single channels
0022     nA=size(A,1);
0023     nB=size(B,1);
0024     tmp=nA-nB;
0025     if ~centerflag
0026         B=[zeros(tmp,1);B];
0027     else
0028         B=[zeros(floor(tmp/2),1);B;zeros(tmp-floor(tmp/2),1)];
0029     end
0030     C=xcorr(A,B);
0031     C=C(nB-1+(1:tmp));
0032     N=cumsum(A.^2);
0033     N(nB+1:end) = N(nB+1:end) - N(1:end-nB);
0034     N=N*sum(B.^2); 
0035     N=N(nB+1:end);
0036     C=C./sqrt(N);
0037     C(find(isnan(C)))=0;
0038     [~,idx]=max(C);
0039     
0040     if nargout==0;
0041         if ~centerflag; 
0042             abscissa=0:size(C,1)-1;
0043         else
0044             abscissa=(0:size(C,1)-1) - ceil(size(C,1)/2);
0045         end
0046         plot(abscissa,C); hold on
0047         plot(abscissa(idx),C(idx),'.k'); hold off
0048         axis tight; xlabel('lag (samples)');
0049         C=[];
0050     end
0051     
0052 else
0053     
0054     % multiple channels
0055     C=zeros(size(A,1)-size(B,1),size(A,2),size(B,2));
0056     idx=zeros(size(A,2),size(B,2));
0057     for k=1:size(A,2)        
0058         for j=1:size(B,2)
0059             [a,b]=nt_xxcorr(A(:,k),B(:,j),centerflag);
0060             C(:,k,j)=a;
0061             idx(k,j)=b;
0062         end
0063     end
0064     
0065     if nargout==0
0066        imagescc(idx-MAXLAG-1);
0067        colorbar
0068        C=[]; idx=[];
0069     end
0070     
0071 end
0072 
0073 % tests
0074 if 0
0075     x=ones(1000,1);
0076     nt_xxcorr(x,x);
0077 end
0078 
0079 if 0 
0080     x=ones(1000,1);
0081     y=ones(500,1);
0082     nt_xxcorr(x,y);
0083 end  
0084 
0085 if 0 
0086     x=randn(1000,1);
0087     y=x(1:500);
0088     x=x.*(1:1000)'.^2;
0089     nt_xxcorr(x,y);
0090 end