function [t,ccf,scot,phat,roth] = gccf(x,y)
% Routine for computing generalized correlation functions.
% Returns: standard ccf, SCOT, PHAT, and Roth functions.
%
% by D. Holmgren, Aug 5 94. holmgren@phobos.astro.uwo.ca

nl = input('Enter FFT length: ');
wx = zeros(nl,1); wy = wx;
nx = length(x); ny = length(y);
wx(1:nx) = x; wy(1:ny) = y;
% generate delay axis...
t = nl/2 - (0:1:nl-1)';
% compute transforms...
X = fft(wx);
Y = fft(wy);
% shift and smooth spectra...
Xs = fftshift(X);
Ys = fftshift(Y);
Xsm = smooth(Xs);
Ysm = smooth(Ys);
% compute various power spectra...
Gxy = Xsm .* conj(Ysm);
Gxym = sqrt(real(Gxy).^2 + imag(Gxy).^2);
Gxx = Xsm .* conj(Xsm);
Gyy = Ysm .* conj(Ysm);
% compute generalised correlations...
ccf = ifft(X .* conj(Y)) ./ sqrt(2048);
ccf = fftshift(ccf);
scot = ifft(fftshift(Gxy ./ sqrt(Gxx .* Gyy)));
scot = fftshift(scot);
phat = ifft(fftshift(Gxy ./ Gxym));
phat = fftshift(phat);
roth = ifft(fftshift(Gxy ./ Gxx));
roth = fftshift(roth);