0001 function [C,idx]=nt_xxcorr(A,B,MAXLAG)
0002
0003
0004
0005
0006
0007
0008
0009
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
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
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
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