11.txt

四阶统计求时延,拿来和大家共同研究

%四阶统计求时延
function minnum=FourthTimeDelay(x,y,max_delay,nfft,nsamp,overlap,yk,zk,plotflag) 
% 使用Tuganit的代价函数的方法做基于四阶统计量的时延估计
% max_delay 最大时延数
% nsamp 分段数据的数据量 默认为0
% overlap 每段数据的数据重叠数 (0,99)默认为50 
% nfft fft计算数据长度
% plotflage 图形显示标志位 默认为1
% yk1,zk2 x在计算4阶累积量的变化范围
% ----------------------- Parameter checks ---------------------------
[lx, lrecs] = size(x);
[ly, nrecs] = size(y); 
if (lx ~= ly | lrecs ~= nrecs )
error(' x, y should have identical dimensions')
end

if (ly == 1)
x = x(:); y = y(:); ly = nrecs; nrecs = 1;
end

if (exist('plotflag') ~= 1) plotflag = 1; end
%if (exist('nfft') ~= 1) nfft = 128; end
if (exist('overlap') ~= 1) overlap = 50; end
overlap = min(99, max(overlap,0));
if (nrecs > 1) overlap = 0; end
if (exist('nsamp') ~= 1) nsamp = 0; end
if (nrecs > 1) nsamp = ly; end
if (nrecs == 1 & nsamp <= 0)
nsamp = fix(ly/ (8 - 7 * overlap/100));
end
%if (nfft < nsamp) nfft = 2^nextpow2(nsamp); end
if (exist('plotflag') ~= 1) plotflag = 1; end
if (exist('k1') ~= 1) k1 = 0; end
if (exist('k2') ~= 1) k2 = 0; end
overlap = fix(overlap/100 * nsamp);
nadvance = nsamp - overlap;
nrecs = fix ( (ly*nrecs - overlap) /nadvance);
%------------四阶累积量------------------------------ 
xxxx=zeros(2*max_delay+1,2*max_delay+1,2*max_delay+1);
xyxx=zeros(2*max_delay+1,2*max_delay+1,2*max_delay+1);
locseg = [1:nsamp]';
for krec = 1:nrecs
xseg = x(locseg);
   yseg = y(locseg); 
for k2=-zk:zk
for k1=-yk:yk
xxxx(:,k1+yk+1,k2+zk+1)=xxxx(:,k1+yk+1,k2+zk+1)+cum4est(xseg,max_delay,nsamp,overlap,'biased',k1,k2);
xyxx(:,k1+yk+1,k2+zk+1)=xyxx(:,k1+yk+1,k2+zk+1)+cum4x(xseg,yseg,xseg,xseg,max_delay,nsamp,overlap,'biased',k1,k2);
end
end
locseg = locseg + nadvance;
end
xxxx=xxxx/nrecs;
xyxx=xyxx/nrecs;
save cum4date xxxx xyxx; 
%-------四阶累积谱-------------------------------------------
tirspxxxx=fftn(xxxx);
tirspxyxx=fftn(xyxx);
G=fftshift(tirspxyxx./tirspxxxx);
G=sum(sum(G));%二次积分
ctau=fftshift(real(ifft(G)));
%-----------------时延寻找-----------------------------------
[val, delay] = max(abs(ctau));
delay = delay - max_delay/2-1 ;
disp(['Delay estimated by TDEB is ',int2str(delay)])
for t=1:2*max_delay
gg(t)=ctau(:,:,t);
end
plot(1:2*max_delay ,gg,1:2*max_delay,gg,'o'), grid on 
title(['TDEB: Hologram delay = ',int2str(delay)])
set(gcf,'Name','Hosa TDEB')

作下来误差很大，不稳定。速度超慢（慢在求累积量）。
算法是利用四阶互谱除以四阶自谱然后做两次积分，再求最大值。
我的用法有错误吗？
请大家帮帮！
谢谢