互相关延时估计加权函数性能分析(共12页)

1、互相关延时估计加权函数(hnsh)性能分析广义互相关函数法是通过首先求出俩信号之间的互功率谱，然后在频域内给予一定的加权，以此对信号和噪音进行白化处理，从而增强信号中信噪比较高的频率成分，抑制(yzh)噪声的影响，最后再反变换到时域，得到两信号之间的互相关函数，即： （1）其中(qzhng)是广义互相关加权函数。广义互相关加权函数的选择主要基于俩个方面：噪声和反射情况。根据不同的情况选择加权函数，其目的就是使具有比较尖锐的峰值。峰值处就是俩个传感器之间的时延。由于来自同一声源的信号存在一定的相关性，通过计算不同麦克风所接受到的信号之间的相关函数，就可以估计出TDOA值。然而在实际环境中，由于噪

2、声和混响的影响，相关函数的最大峰会被弱化，有时还会出现多个峰值，这些都造成了实际峰值的检测困难。此时就通过加权的方法来锐化峰值，通常我们通过时间、精度来确定算法的合理性。广义互相关函数模拟clear all; clc; close all;N=1024;%长度Fs=500;%采样频率n=0:N-1;t=n/Fs; %时间序列a1=5; %信号幅度a2=5;d=2; %延迟点数x1=a1*cos(2*pi*10*n/Fs); %信号1x1=x1+randn(size(x1); %加噪声x2=a2*cos(2*pi*10*(n+d)/Fs); %信号2x2=x2+randn(size(x2);su

3、bplot(211);plot(t,x1,r);axis(-0.2 1.5 -6 6);hold on;plot(t,x2,:);axis(-0.2 1.5 -6 6);legend(x1信号, x2信号);xlabel(时间/s);ylabel(x1(t) x2(t);title(原始信号);grid on;hold off%互相关函数X1=fft(x1,2*N-1);X2=fft(x2,2*N-1);Sxy=X1.*conj(X2);Cxy=fftshift(ifft(Sxy);%Cxy=fftshift(real(ifft(Sxy);subplot(212);t1=(0:2*N-2)/F

4、s; %注意plot(t1,Cxy,b);title(互相关函数);xlabel(时间/s);ylabel(Rx1x2(t);grid onmax,location=max(Cxy);%求出最大值max,及最大值所在的位置（第几行）location;%d=location-N/2-1 %算出延迟了几个点d=location-NDelay=d/Fs %求得时间延迟运行程序得到(d do)的结果是：d = 2Delay = 0.0040可以看出，通过(tnggu)互相关函数的求解d=2，delay=0.0040，这和我们给出的信号的时延d/Fs=0.0040是一致的。这表明互相关函数可以给出信号的

5、时延估计。PHAT-GCC模拟(mn)clear all; clc; close all;N=1024; %长度(chngd)Fs=500; %采样(ci yn)频率n=0:N-1;t=n/Fs; %时间(shjin)序列a1=5; %信号幅度a2=5;d=9; %延迟点数x1=a1*cos(2*pi*10*n/Fs); %信号1x1=x1+randn(size(x1); %加噪声%x1=x1.*hamming(max(size(x1);%加窗x2=a2*cos(2*pi*10*(n+d)/Fs); %信号2x2=x2+randn(size(x2);%x2=x2.*hamming(max(si

6、ze(x2);%加窗subplot(211);plot(t,x1,r);axis(-0.2 2 -6 6);hold on;plot(t,x2,:);axis(-0.2 2 -6 6);legend(x1信号, x2信号);xlabel(时间/s);ylabel(x1(t) x2(t);title(原始信号);grid on;hold off%互相关函数X1=fft(x1,2*N-1);X2=fft(x2,2*N-1);Sxy=X1.*conj(X2);%Cxy=fftshift(ifft(Sxy);Cxy=fftshift(ifft(Sxy./abs(Sxy);subplot(212);t1

7、=(-N+1:N-1)/Fs;plot(t1,Cxy,b);title(Rx1x2);xlabel(t/s);ylabel(Rx1x2(t);grid onmax,location=max(Cxy);%d=location-N/2-1 d=location-NDelay=d/Fs %求得时间延迟运行程序得到的结果是：d = 1 Delay = 0.0020我们可以看见结果是d=1,delay=0.0020，而实例中给出的时延为d/fs=0.016，这并不表示PHAT-GCC算法是错误的，只是(zhsh)因为，我们在信号中加入了均值为0，方差为1的高斯白噪音，所以才会导致了误差的存在。ROTH-

8、GCC模拟(mn)clear;N=1024;%信号(xnho)长度fs=500;%采样频率n=0:N-1;t=n/fs;%时间序列a1=5;%信号幅度a2=5;%信号幅度d=2;%延迟点数x1=a1*sin(2*pi*10*n/fs);x2=a2*sin(2*pi*10*(n+d)/fs);%x2=awgn(x1./4,-3); %噪声强度大于信号%x2=x2 .* hamming(N);x1=x1+randn(size(x1); %加入噪声x2=x2+randn(size(x2);S1=fft(x1,2*N-1);S2=fft(x2,2*N-1);S12 = S1.* conj(S2);S1

9、1 = S1.* conj(S1);R1 =real(fftshift(ifft(S12./abs(S11);ts=(-N+1:N-1)/fs;plot(ts,R1);xlabel(时间(shjin)/s);ylabel(R1(t);title(互相(h xing)关函数);max,location=max(R1);%d=location-N/2-1 d=location-NDelay=d/fs 运行(ynxng)程序的结果为：d = 4Delay =0.0080SCOT-GCC模拟clear;N=1024;%信号(xnho)长度fs=1000;%采样(ci yn)频率n=0:N-1;t=n/

10、fs;%时间(shjin)序列ts = 1/fs * (-N + (1 : 2*N - 1); %互相关时间序列a1=5;%信号幅度a2=5;%信号幅度d=26;%延迟点数x1=a1*sin(2*pi*10*t)+1.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(x1./4,-3); %噪声强度大于信号%x2=x2 .* hamming(N);x=awgn(x1,20); %加入噪声y=awg

11、n(x2,0.001);S1=fft(x,2*N-1);S2=fft(y,2*N-1);X = S1.* conj(S2);X11 = S1.* conj(S1);X22 = S2.* conj(S2);Y=sqrt(X11.*X22);R1 =real(fftshift(ifft(X./Y);plot(ts,R1);xlabel(时间/s);ylabel(R1(t);title(ifft计算结果)max,location=max(R1);%d=location-N/2-1 d=location-NDelay=d/fs 运行程序的结果是：d = 8Delay =0.0080相同信噪比不同(b

12、tn)算法的比较clear all; clc; close all;N=1024; %长度(chngd)Fs=500; %采样(ci yn)频率n=0:N-1;t=n/Fs; %时间序列a1=30; %信号幅度a2=30;d=9; %延迟点数x1=a1*cos(2*pi*10*n/Fs); %信号1x1=awgn(x1,20); %加噪声%x1=x1.*hamming(max(size(x1);%加窗x2=a2*cos(2*pi*10*(n+d)/Fs); %信号2x2=awgn(x2,20);%x2=x2.*hamming(max(size(x2);%加窗subplot(511);plot(

13、t,x1,r);axis(-0.2 2 -40 40);hold on;plot(t,x2,:);axis(-0.2 2 -40 40);legend(x1信号, x2信号);xlabel(时间(shjin)/s);ylabel(x1(t) x2(t);title(原始(yunsh)信号);grid on;%互相(h xing)关函数ticX1=fft(x1,2*N-1);X2=fft(x2,2*N-1);Sxy=X1.*conj(X2);%Cxy=fftshift(ifft(Sxy);%GCCCxy=fftshift(ifft(Sxy);subplot(512);t1=(-N+1:N-1)/

14、Fs;plot(t1,Cxy,b);title(GCC);xlabel(t/s);ylabel(Cxy);grid on;max1,location1=max(Cxy);%d=location-N/2-1 d1=location1-NDelay1=d1/Fs %求得时间延迟toc%phat-gccticX1=fft(x1,2*N-1);X2=fft(x2,2*N-1);Sxy=X1.*conj(X2);Pxy=fftshift(ifft(Sxy./abs(Sxy);subplot(513);t1=(-N+1:N-1)/Fs;plot(t1,Pxy,b);title(phat-gcc);xlab

15、el(t/s);ylabel(Pxy);grid on;max2,location2=max(Pxy);%d=location-N/2-1 d2=location2-NDelay2=d2/Fs %求得时间延迟toc%rhat-gccticX1=fft(x1,2*N-1);X2=fft(x2,2*N-1);Sxy=X1.*conj(X2);S11 = X1.* conj(X1);Rxy=fftshift(ifft(Sxy./abs(S11);subplot(514);t1=(-N+1:N-1)/Fs;plot(t1,Rxy,b);title(phat-gcc);xlabel(t/s);ylabe

16、l(Rxy);grid on;max3,location3=max(Rxy);%d=location-N/2-1 d3=location3-NDelay3=d3/Fs %求得时间延迟toc%scot-gccticX1=fft(x1,2*N-1);X2=fft(x2,2*N-1);Sxy=X1.*conj(X2);S11 = X1.* conj(X1);S22 = X2.* conj(X2);Y=sqrt(S11.*S22);SCxy=fftshift(ifft(Sxy./Y);subplot(515);t1=(-N+1:N-1)/Fs;plot(t1,SCxy,b);title(scot-gc

17、c);xlabel(t/s);ylabel(SCxy);grid on;max4,location4=max(SCxy);%d=location-N/2-1 d4=location4-NDelay4=d4/Fs %求得时间延迟tocSNR=0时SNR=10时：SNR=20时SNR=50时分析(fnx)：从运行(ynxng)结果上来看，在时间上，基本(jbn)互相关、PHAT加权、ROTH加权和SCOT加权四种算法的运行时间基本相同；但是从峰度的锐化来说，这四种方式的时延估计的准确性随着信噪比的降低而恶化，互相关函数峰值的尖锐程度随信噪比的降低而降低。对于SCOT加权来说，随着信噪比的降低，性能急剧下降。基本互相关函数和RHOT加权虽然有一定的抗噪能力，但随着信噪比的降低，其波动程度明显加强，特别是对外围的噪声、反射和有限观测数据很敏感，会造成峰值不明显；对于PHAT加权，