互相关延时估计加权函数性能分析

1、互相关延时估计加权函数性能分析广义互相关函数法是通过首先求出俩信号之间的互功率谱，然后在频域给予 一定的加权，以此对信号和噪音进行白化处理，从而增强信号噪比较高的频率成 分，抑制噪声的影响，最后再反变换到时域，得到两信号之间的互相关函数，即:凡2 =£ 内2(卬)G?( W)"''"卬(1)其中心2(卬)是广义互相关加权函数。广义互相关加权函数的选择主要基于俩 个方面：噪声和反射情况。根据不同的情况选择加权函数，其目的就是使叫2«) 具有比较尖锐的峰值。与2«)峰值处就是俩个传感器之间的时延。由于来自同一声源的信号存在一定的相关

2、性，通过计算不同麦克风所接受到 的信号之间的相关函数，就可以估计出TDOA值。然而在实际环境中，由于噪声 和混响的影响，相关函数的最大峰会被弱化，有时还会出现多个峰值，这些都造 成了实际峰值的检测困难。此时就通过加权的方法来锐化峰值，通常我们通过时 间、精度来确定算法的合理性。一、广义互相关函数模拟clear all； clc； close all；N=1024；% 长度Fs=500；%采样频率n=0：N-l；t=n/Fs； %时间序列al=5；%信号幅度a2=5；d=2；%延迟点数xl=al*cos(2*pi*10*n/Fs) ；%信号 1xl=xl+randn(size(xl) ；%加噪声

3、x2=a2*cos (2*pi*10*(n+d) /Fs) ； %信号 2x2=x2+randn(size(x2)；subplot(211)；plot(tFxl,frr);axis(-0.2 1.5 -6 6)；hold on；plot(t,x2,'：')；axis(-0.2 1.5 -6 6)；legend ('xl 信号'，'x2 信号')；xlabel ('时间/s') ；ylabel (rxl (t) x2(t)')；title('原始信号')；grid on；hold off%互相关函数Xl=ff

4、t(xl,2*N-l)；X2=fft(x2t2*N-l)；Sxy=Xl.*conj (X2)；Cxy=fftshift(ifft(Sxy)；%Cxy=fftshift(real(ifft(Sxy)；subplot(212)；tl=(0：2*N-2)/Fs；% 注意plot (tl,Cxy,'b')；title('互相关函数')；xlabel ('时间/s') ；ylabel ('Rxlx2(t) ') ；grid onmax, location =max(Cxy);%求出最大值max,及最大值所在的位置(第几行) location

5、；%d=location-N/2-l先算出延迟了几个点d=location-NDelay=d/Fs/求得时间延迟运行程序得到的结果是：d =2Delay =0. 0040原始信号64©NX©Lx20- 2- 4- 6-0.200.20.40.60.811.21.4X1O4互藉具篇数10-1-200,511522.533.5445时间看- exh可以看出，通过互相关函数的求解d=2, delay=0. 0040,这和我们给出的信 号的时延d/Fs=0.0040是一致的。这表明互相关函数可以给出信号的时延估计。 二、 PHAT-GCC 模拟 clear all； clc； cl

6、ose all；N=1024； % 长度Fs=500； %采样频率n=O：N-l；t=n/Fs； %时间序列al=5; %信号幅度a2=5；d=9； %延迟点数xl=al*cos(2*pi*10*n/Fs) ；%信号 1xl=xl+randn (size (xl) ；%加噪声%xl=xl. *hamming(max(size(xl) '加窗 x2=a2*cos(2*pi*10*(n+d) /Fs) ； %信号 2 x2=x2+randn(size(x2)；%x2=x2. *hamming(max(size(x2) ' ；%加窗subplot(211)；plot，');a

7、xis(-0.2 2 -6 6)；hold on；plot(t,x2,': r);axis(-0.2 2 -6 6)；legend ('xl 信号'，'x2 信号')；xlabel ('时间/s') ； ylabel (r xl (t) x2(t)')； title('原始信号')；grid on； hold off%互相关函数Xl=fft(xlT2*N-l)； X2=fft(x2,2*N-l)； Sxy=Xl.*conj(X2)； %Cxy=fftshift(ifft(Sxy)； Cxy=fftshift(iff

8、t(Sxy. /abs(Sxy)； subplot(212);tl=(-N+l：N-l)/Fs； plot(tl,Cxy,rbr)；title ('Rxlx2")；xlabel('t/s1)；ylabel('Rxlx2(t) ') ； grid on max,location=max(Cxy)； %d=location-N/2-l d=location-NDelay=d/Fs%求得时间延迟运行程序得到的结果是： d =1Delay =0. 00206.2 4 2 CD-2B 上-O. 2NXslx4o.o.6.81 .614 di2-2.5.2.1.5

9、-0.500.511.522.5t/s0.050-0.05我们可以看见结果是d=l,delays. 0020 ,而实例中给出的时延为 d/fs=0.016,这并不表示PHAT-GCC算法是错误的，只是因为，我们在信号中加 入了均值为0,方差为1的高斯白噪音，所以才会导致了误差的存在。3、 ROTH-GCC模拟clear；N:1024；%信号长度fs=500；%采样频率n=0：N-l；t二n/fs；%时间序列al=5； %信号幅度a2=5；%信号幅度d=2；%延迟点数xl=al*sin(2*pi*10*n/fs);x2=a2*sin(2*pi*10*(n+d)/fs)；%x2=awgn(xl.

10、/4,-3) ；%噪声强度大于信号%x2=x2 . * hamming(N)；%加入噪声xl=xl+randn(size(xl)； x2=x2+randn(size(x2)；Sl=fft(xl,2*N-l)； S2=fft(x2,2*N-l)；S12 = SI.* con j (S2)；Sil = SI.* conj(Sl)；RI =real(fftshift(ifft(S12. /abs(Sil)；ts=(-N+l：N-l)/fs；plot(ts,Rl);xlabel ('时间/s') ； ylabel (, RI (t),)； title('互相关函数')；

11、max,location=max(RI)； %d=location-N/2-l d=location-NDelay=d/fs 运行程序的结果为： d =4Delay =0. 00804 2 CD-2B-6-06 2NX WK I原始信号 6.6 4 O.2 O,8164112-2.6-2-1.5-1-0.500.511.522.6时间居04、 SC0T-GCC模拟clear；N:1024；%信号长度fs=1000；%采样频率n=O：N-l；t=n/fs；%时间序列ts = 1/fs * (-N + (1 ： 2*N - 1)； %互相关时间序列a 1=5； %信号幅度a2=5；%信号幅度d=2

12、6；%延迟点数xl=al*sin(2*pi*10*t)+l. 9*sin(2*pi*18*t)+2. 8*sin(2*pi*55*t)；x2=a2*sin(2*pi*10*(n+d)/fs)+1. 9*sin(2*pi*18*(n+d)/fs)+2. 8*sin(2*pi*55*(n+d)/fs)；%x2=awgn(xl. /4,-3) ；%噪声强度大于信号%x2=x2 . * hamming(N)；x=awgn(xl,20) ；%加入噪声y=awgn(x2,0. 001)；Sl=fft(x,2*N-l);S2=fft(y,2*N-l)；X = SI.* con j (S2)；Xll = SI

13、.* conj (SI)；X22 = S2. * conj (S2)；Y=sqrt(Xll. *X22)；RI =real(fftshift(ifft(X. /Y)；plot(ts,RI)；xlabel ('时间/s'); ylabel (' RI (t),)；title(* ifft 计算结果')max,location=max(RI)；%d=location-N/2-1d=location-NDelay=d/fs 运行程序的结果是：d =8Delay =0. 0080原始信号eNX果0.2-0.1.5-0.2-1.5-0.6时间抬0.1五、相同信噪比不同算法

14、的比较clear all； clc； close all；N=1024； % 长度Fs=500；/采样频率n=0：N-l；t=n/Fs； %时间序列al=30； %信号幅度a2=30；d=9；%延迟点数xl=al*cos(2*pi*10*n/Fs) ；%信号 1xl=awgn(xlf20) ；%加噪声%xl=xl. *hamming(max(size(xl)'加窗 x2=a2*cos (2*pi*10*(n+d)/Fs) ； %信号 2 x2=awgn(x2,20)；%x2=x2. "hamming(max(size(x2)'%加窗 subplot (511)；plo

15、t (t,xl,'r')；axis (-0.2 2 -40 40)；hold on；plot(t,x2,'axis(-0. 2 2 -40 40)；legend('xl 信号'，'x2 信号')；xlabel ('时间/s') ； ylabel (* xl (t) x2 (t)1); title('原始信号')；grid on；%互相关函数 tic Xl=fft(xl,2*N-l)； X2=fft(x2,2*N-l)； Sxy=Xl.*conj(X2)； %Cxy=fftshift(ifft(Sxy)； %

16、GCCCxy=fftshift(ifft(Sxy)；subplot(512)；tl=(-N+l：N-l)/Fs；plot(tl,Cxy,1b')；title ('GCC')；xlabel('t/s*)；ylabel('Cxy') ； grid on；Imaxl,locationl=max(Cxy)；%d=location-N/2-1 dl=locationl-N Delayl=dl/Fs%求得时间延迟toe %phat-gcc ticXl=fft(xl,2*N-l)；X2=fft(x2,2*N-l)；Sxy=Xl.*conj (X2)；Pxy=f

17、ftshift(ifft(Sxy. /abs(Sxy)； subplot(513)； tl=(-N+l：N-l)/Fs； plot(tl,Pxy,'b')；title('phat-gcc')；xlabel('t/s')；ylabel('Pxy1) ； grid on；max2,location2=max(Pxy)； %d=location-N/2-ld2=location2-NDe 1 ay2=d2/Fs%求得时间延迟toe %rhat-gcc ticXl=fft(xl,2*N-l)； X2=fft(x2,2*N-l)； Sxy=Xl.*

18、conj(X2)；Sil = XL* conj (XI)；Rxy=fftshift(ifft(Sxy. /abs(SI 1)； subplot(514)； tl=(-N+l：N-l)/Fs； plot(tlFRxy,'b')； title(，phat-gcc，)；xlabel('t/s')；ylabel('Rxy') ；grid on；Imax3t1ocation3=max(Rxy)； %d=locati on-N/2T d3=location3-NDe 1 ay3=d3/Fs%求得时间延迟toe%SCOt-gCC tic Xl=fft(xl,2

19、*N-l)； X2=fft(x2,2*N-l)； Sxy=Xl.*conj(X2)； Sil = XL* conj (XI)； S22 = X2. * conj (X2)；Y=sqrt(Sll. *S22);SCxy=fftshift(ifft(Sxy. /Y)； subplot(515)； tl=(-N+l：N-l)/Fs； plot (tl,SCxy,'b')； title(，scot-gcc，)；xlabel('t/s')；ylabel('SCxy') ；grid on；Imax4,location4=max(SCxy)； %d=locat

20、ion-N/21 d4=location4-NDelay4=d4/Fs%求得时间延迟toe SNR=0 时原始信号至 croo oo42 2 4 .(wx (_=xo 石WAxoi011111111111.0111111111111112-1.5-1-0.5 H Cy 0.511.522-1.5-1-0.5 Uy 0.511.5226Illi IlliIlli Illi 111Illi IlliIlliIlli IlliIlli-2-1.6-1-0.500.511.52t/s2.62.5SNR=10 时:原始信号ooooonu 42-2-4 WNX (MX00.20.4 0.6 O.F 11.21.41.61.80-2-1.5-1-0.5 p Dy 0.511.522.50.5杂x n 。u