gccstandard.m

广义的互相提取时延算法,此程序是基于matlab的源程序,适用于非高斯信号和高斯白噪声的情况

%吴逢吉,无线定位TOA/TDOA参估计方法的研究,仿真问题(非高斯噪声白信号+不相关白噪声)(13-17)
clear;clc;
t=-1023:1023;
sx=-1+2.*rand(1,1024);%产生均值0,方差1/1.732,[-1,1]均匀分布的随机序列.
sx1=1000.*sx';
sy1=1.5.*circshift(sx1,1004);%延时20单位sy1=6.*circshift(sx1,-1004);%延时20单位
sy=sy1';
x1=awgn(sx,0,'measured');%加入SNR为零的正态分布白噪声随机序列,'measured'表测量信号噪声.x1=awgn(sx,0,'measured','db' )
y1=awgn(sy,0,'measured' );%第二种改法y1=awgn(sy,0,'measured','db' );
rxx=xcorr(x1,x1,'biased');%求信号x1自相关函数,度应是1024
ryy=xcorr(y1,y1,'biased');%求信号y1自相关函数
rxy=xcorr(x1,y1,'biased');%求信号x1,y1互相关函数[rxx,lags]=xcorr(x1,x1,'biased')
sxx=fft(rxx,2047);%求x1功率谱密度
syy=fft(ryy,2047);%求y1功率谱密度
sxy=fft(rxy,2047);%求互相关谱密度,
cxy=(sxy.* sxy)/(sxx.*syy);
Hs=cxy./((abs(sxy)).*(1-cxy));
sxyf=Hs.*sxy;
TDOA=ifft(sxyf,2047);%注2047必与t的长度一致,否则会报错.
GCCA=(TDOA.*conj(TDOA)).^0.5;
t1=-30:0.01:30;
final=interp1(t,GCCA,t1,'spline');%插值
plot(t1,final)
axis([-40,40,-0.1,0.2])
grid on