%
disp(' ');
disp(' cross_correlation.m   ver 1.2  June 19, 2013 ');
disp(' by Tom Irvine  Email: tom@vibrationdata.com ');
%
clear ac1;
clear ac2;
clear ac3;
clear a;
clear b;
clear xc;
clear x;
clear y;
clear THM;
clear length;
clear t;
%
disp(' ');
disp(' Select signal format ');
disp('  1= Common array:  time(sec) & amplitude & amplitude ');
disp('  2= Two separate arrays, each: time(sec) & amplitude ');
%
iform=input(' ');
%
disp(' ');
%
if(iform==1)
    disp(' The signal array must have three columns: ');
    disp('   time(sec) & amplitude & amplitude ');
    disp(' ');  
    THM=input(' Enter the array name: ');
    t=THM(:,1);
    x=THM(:,2);
    y=THM(:,3);
    n=length(t);
    m=n;
    dur=t(n)-t(1);
    dt=dur/(n-1);
%
    disp(' ');
    out1=sprintf(' dt = %8.4g sec   n=%d',dt,n);
    disp(out1);
%
    a=[t x];
    b=[t y];
else
%
    disp(' Each signal must have two columns:  time(sec) & amplitude ');
    disp(' ');
%
    a=input(' Enter the  first signal: ');
    b=input(' Enter the second signal: ');
%
    m=max(size(a));
    n=max(size(b));
%
    dur1=a(m,1)-a(1,1);
    dur2=b(m,1)-b(1,1);
    dur=dur1;
    if(dur2<dur)
        dur=dur2;
    end
%
    dt1=(a(m,1)-a(1,1))/(m-1);
    dt2=(b(n,1)-b(1,1))/(n-1);
%
    disp(' ');
    out1=sprintf(' dt1 = %8.4g sec ',dt1);
    out2=sprintf(' dt2 = %8.4g sec ',dt2);
%
    disp(out1);
    disp(out2);
%
%  Truncate if necessary
%
    if(m<n)
        n=m;
    end
    if(n<m)
        m=n;
    end
%
    pe=abs((dt1-dt2)/dt1);
%
    dt=(dt1+dt2)/2;
%
    if(pe>0.01)
        out1=sprintf('Warning: dt1=%8.4g  dt2=%8.4g ',dt1,dt2);
        disp(out1);
    end    
%
end
%
disp(' ');
disp(' Select calculation method: ')
disp('  1 = Time domain ')
disp('  2 = ifft(fft(x).*conj(fft(y))) ');
method=input(' ');
%
if(method==1)
        xc=zeros(2*n-1,2);
        xc(:,1)=linspace(-dur,dur,(2*n-1));
        n2=2*n;
        progressbar;
%
        clear Z;
        Z=zeros(n,1);
        p=n;
        for i=1:n
            progressbar(i/n);
            clear X;
            clear Y;
            X=a((1:p),2);
            Y=b((i:n),2);
            Z(i)=sum(X.*Y);
            p=p-1;
        end
        Z=Z/n;
%       
        pause(0.1);
        progressbar(1,1);    
%
        sz=size(Z);
        if(sz(2)>sz(1))
            Z=Z';
        end
        Zr=flipud(Z);
        Zr(n)=[];
        xc(:,2)=[ Zr ; Z]; 
%
else
%  
        disp(' ref 1');
        p=2*n-1;
        ac1 = fft(b(:,2),p) .* conj(fft(a(:,2),p));
        ac1=real(ifft(ac1));
        ac1=ac1/n;
        sz=size(ac1);
        if(sz(2)>sz(1))
            ac1=ac1';
        end
        m=floor(length(ac1)/2);
        ac2 = ac1(1:m);
        ac2r=flipud(ac2);
        ac2r(m)=[];
        ac3=[ ac2r ; ac2];
        n=length(ac3);
        t=linspace(0,(n-1)*dt,n);
        t=t-t(n)/2;
        [~,I] = min(abs(t));
        t(I)=0;
        sz=size(t);
        if(sz(2)>sz(1))
            t=t';
        end
        xc=[t,ac3];        
%  
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    clear C;
    clear I;
    clear xcc;
    xcc=[xc(:,1) abs(xc(:,2))];
    [C,I]=max(xcc);
    xmax=xcc(I(2),2);
    tmax=abs(xc(I(2),1));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    figure(1);
    plot(a(:,1),a(:,2));
    grid on;
    xlabel('Time(sec)');
    title(' Signal 1 ');
%
    figure(2);
    plot(b(:,1),b(:,2));
    grid on;
    xlabel('Time(sec)');
    title(' Signal 2 ');
%    
    disp(' ');
    out1=sprintf(' Maximum Absolute occurs at %8.4g sec ',tmax);
    disp(out1);
%
    figure(3);
    plot(xc(:,1),xc(:,2));
    xlabel(' Delay(sec) ');
    out1=sprintf(' Cross-correlation  Max Abs= %8.4g at %8.4g sec',xmax,tmax);
    title(out1);
    grid on;
%