高阶谱分析3，适用于阵列信号处理、雷达信号处理等Matlab.zip

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')