function [delay,ctau] = tdeb (x,y, max_delay,  nfft, wind, nsamp, overlap)
%TDEB	Time Delay Estimation using conventional bispectrum method.
%	[delay,ctau] = tdeb (x,y,  max_delay, nfft, wind, segsamp, overlap)
%	x         - signal at sensor 1
%	y         - signal at sensor 2
%	max_delay - maximum expected delay in samples
%	nfft - fft length [default = power of two >= segsamp]
%	wind - window specification for frequency-domain smoothing
%	       if 'wind' is a scalar, it specifies the length of the side
%	          of the square for the Rao-Gabr optimal window  [default=5]
%	       if 'wind' is a vector, a 2D window will be calculated via
%	          w2(i,j) = wind(i) * wind(j) * wind(i+j)
%	       if 'wind' is a matrix, it specifies the 2-D filter directly
%	segsamp - samples per segment [default: 4*max_delay+1]
%	        - if x is a matrix, segsamp is set to the number of rows
%	overlap - percentage overlap, allowed range [0,99]. [default = 50];
%	        - if x is a matrix, overlap is set to 0.
%	delay     - estimated delay (in samples)
%       ctau      - estimated third-order hologram.

%  Copyright (c) 1991-2001 by United Signals & Systems, Inc. 
%       $Revision: 1.9 $
%  A. Swami   January 20, 1995

%     RESTRICTED RIGHTS LEGEND
% Use, duplication, or disclosure by the Government is subject to
% restrictions as set forth in subparagraph (c) (1) (ii) of the
% Rights in Technical Data and Computer Software clause of DFARS
% 252.227-7013.
% Manufacturer: United Signals & Systems, Inc., P.O. Box 2374,
% Culver City, California 90231.
%
%  This material may be reproduced by or for the U.S. Government pursuant
%  to the copyright license under the clause at DFARS 252.227-7013.

% parameter checks ---------------------------------------------------------

[lx,nrecs] = size(x);
[ly,ny] = size(y);
if (ly ~= lx | ny ~= nrecs)
    error('matrices x and y should have the same dimensions')
end
if (lx == 1), lx = nrecs; nrecs = 1; x = x(:); y = y(:); end

if (exist('max_delay') ~= 1)
   error('max_delay must be specified');
end

if (exist('nsamp') ~= 1) nsamp =  2^nextpow2(4*max_delay+1); end
if (exist('overlap') ~= 1) overlap = 50; end
overlap = min(99,max(overlap,0));

if (exist('wind') ~= 1)       wind = 5; end

if (nrecs > 1)
   overlap = 0;      nsamp   = lx;
end
if (exist('nfft') ~= 1)  nfft  =   0; end
if (nfft < nsamp)
   nfft = 2 ^ nextpow2(nsamp) ;
end

plotflag = 0;
Bxyx = bispecdx (x, y, x, nfft, wind, nsamp, overlap, plotflag) ;
Bxxx = bispecdx (x, x, x, nfft, wind, nsamp, overlap, plotflag) ;

Bxyx = fftshift(Bxyx);     % odd nfft ?
Bxxx = fftshift(Bxxx);

 coher = sum(Bxyx ./ Bxxx);     % sum integrates over omega_2
 ctau  = fftshift(real(ifft(coher)));
 [val, delay] = max(abs(ctau));
 delay = delay - nfft/2 - 1;
 disp(['Delay estimated by TDEB is ',int2str(delay)])
 plot(-nfft/2:nfft/2-1,ctau, -nfft/2:nfft/2-1,ctau,'o'), grid on 
 title(['TDEB: Hologram delay = ',int2str(delay)])
 set(gcf,'Name','Hosa TDEB')