function     [tdc,wdc,pp,rr]=compcor(x,y,doplot)
%function     [tdc,fdc,wdc,pp,rr,ss]=compcor(x,y,doplot) fdc=tdc & pp=ss are redundant!
% Function that finds the correlation coefficient (peak of the correlation function)
% for two signals 
%
% Syntax: [tdc,wdc] = clascor (x,y,doplot);    
%
% where           x        : 1st signal (with or without noise)
%                 y        : 2nd signal (with or without noise)
%            doplot        : flag for cross correlation plots. 1 = plot , 0 = no plot 
% Returns       tdc        : time domain Correlation coeficient using Cross-correlation of the 2 signals
%               fdc        : fourier domain Correlation coeficient using Cross-correlation of the 2 signals 
%               wdc        : wavelet domain Correlation coeficient using Cross-correlation of the 2 signals
%               pp         : time domain Correlation coeficient using PHAT filter 
%               rr         : time domain Correlation coeficient using ROTH filter 
% Example: 
%                 N = 256;
%                 x = MakeSignal('Doppler',N);
%                 y = x + randn(size(x));
%                 [tdc,fdc,wdc]=clascor(x,y,1)
%
% RJB 10/05/05
% Some code segments ripped off from "transclas()" written by Stefanos Kouteas about 1999

%check for plotting
if nargin < 2
    error('Need at least 2 signal inputs')
elseif nargin <3
    doplot =0;
end
  
%check for Phat and Roth test
if nargout == 4
    DoFlag = 1; 
else
    DoFlag = 0;
end
    
x = x - mean(x);                               %remove mean
y = y - mean(y);
xx=x(:);                                      %make column vectors
yy=y(:);
lx = length(xx);
ly = length(yy);
nfft = 2*max(lx,ly);
w=ones(1,nfft);                                %use a rectangle window (for now)
b=nfft;
%-----------------------------------------------------------
tdc= (xx' * yy)/(norm(xx) .* norm(yy));        %Cross-correlation
%----------------------------------------------------------------
%Wavelet Threshold
%qmf = MakeONFilter('Daubechies',20); 
qmf = MakeONFilter('Coiflet',5);
%qmf = MakeONFilter('Haar');
%qmf = MakeONFilter('Symmlet',8);                   %wdc
wx = FWT_PO(xx,1,qmf);                             %wavelet coeff's  
wy = FWT_PO(yy,1,qmf);
wx = ThreshWave3(wx','S',0,0,sqrt(2*log(lx)),1,qmf);
wy = ThreshWave3(wy','S',0,0,sqrt(2*log(ly)),1,qmf);
wdc = (wx * wy')/(norm(wx) .* (norm(wy)) + eps);
%---------------------------------------------------------
if DoFlag ==1
    %---------------------------------------------------------------
    X=fft(xx,nfft);                                %FFT of signal x
    Y=fft(yy,nfft);                                %FFT of signal y

    XX=X.*conj(X);  %XX=sum(XX,2);                  %Auto-spectrum of x
    YY=Y.*conj(Y);  %YY=sum(YY,2);                  %Auto-spectrum of y 
    XY=Y.*conj(X);  %XY=sum(XY,2);                  %Cross-spectrum of x and y
    %--------------------------------------------------------------
    P=XY./(abs(XY)+eps);                          % Phat in frequency    
    P=P.*w(nfft);                                 %The window (it could be more complicated)
    p=ifft(P,nfft);                             %Phat in time
    p=[p((b/2)+1:end);p(1:(b/2))];                %The PHAT correlation function  
    pp=max(p);
    pp = ckreal(pp);
%---------------------------------------------------------------
    R=XY./(XX +eps)*( (norm(XX))/ (norm(XY)));    %Roth in frequency
    R=R.*w(nfft); 
    r=ifft(R,nfft);     %Roth in time
    r=[r((b/2)+1:end);r(1:(b/2))];
    %rr = max(r)./norm(y);
    rr=max(r);
    rr = ckreal(rr);
end
%------------------------------------------------------------
%Same as PHAT!
%S=XY./(sqrt(XX.*YY)+eps);                      %SCOT in Frequency domain
%S=S.*w(nfft);                                  %The window (it could be more complicated)
%s=ifft(S,nfft);                                %SCOT in Time domain
%b=length(s);                                   %Need to 'fft flip' the time segments 
%s=[s((b/2)+1:end);s(1:(b/2))];                 %The SCOT correlation function 
%ss = max(s);
%ss = ckreal(ss);
%--------------------------------------------------------
%Same as TDC
%c=ifft(XY,nfft);                               %Cross-correlation
%c = c./(norm(x).*norm(y));                     %Need to normalize 
%c=[c((b/2)+1:end);c(1:(b/2))];                 %Need to 'fft flip' the time segments 
%fdc = max(c);
%fdc = ckreal(fdc);                             %fdc
%----------------------------------------------------------------
if doplot ==1
    dt = 1/nfft;                                   %change this if desired
    nfft=nfft-mod(nfft,2);
    n=(-(nfft/2):nfft/2-1)*dt;
    figure(1)
    subplot(221),plot(n,tdc), title('SCOT Cross-Corr.')
    subplot(222),plot(n,wdc), title('PHAT Cross-Corr.')
    subplot(223),plot(n,pp), title('ROTH Cross-Corr.')
    subplot(224),plot(n,rr), title('Cross Corr.')
end
%-------------------------------------------------------------------------------------
function y =ckreal(x)
%check if complex number is actually a real number
%matlab gives ss,pp,rr,cc as complex numbers with zero imaginary parts
%it doesn't effect anything, but it is annoying.

xr = real(x); xi = imag(x); xp = xi./xr;
if (xi < 1e-6) & (xp < 1e-6)
    y = real(x);
else
    y = x;
end
%-----------------------------------------------------------------------------------------