互相关延时估计加权函数性能分析教材.docx

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

3、subplot(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)

4、/Fs; %注意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 = 2Delay = 0.0040可以看出，通过互相关函数的求解d=2，delay=0.0040，这和我们给出的信号的时延d/Fs=0.0040是一致的。这表明互相关函数可以给出信号的时延估计。二、 PHA

5、T-GCC模拟clear all; clc; close all;N=1024; %长度Fs=500; %采样频率n=0:N-1;t=n/Fs; %时间序列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(size(x2);%加窗subplot(211);plot(t,x1,

6、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=(-N+1:N-1)/Fs;plot(t1,Cxy,b);tit

7、le(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算法是错误的，只是因为，我们在信号中加入了均值为0，方差为1的高斯白噪音，所以才会导致了误差的存在。三、 ROTH-GCC模拟clear;N=1024;%信号长度fs=500;%采样频率n

8、=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);S11 = S1.* conj(S1);R1 =real(fftshift(ifft(S12./

9、abs(S11);ts=(-N+1:N-1)/fs;plot(ts,R1);xlabel(时间/s);ylabel(R1(t);title(互相关函数);max,location=max(R1);%d=location-N/2-1 d=location-NDelay=d/fs 运行程序的结果为：d = 4Delay =0.0080四、 SCOT-GCC模拟clear;N=1024;%信号长度fs=1000;%采样频率n=0:N-1;t=n/fs;%时间序列ts = 1/fs * (-N + (1 : 2*N - 1); %互相关时间序列a1=5;%信号幅度a2=5;%信号幅度d=26;%延迟点

10、数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=awgn(x2,0.001);S1=fft(x,2*N-1);S2=fft(y,2*N-1);X = S1.* conj(S2);X11 = S1.* conj(S1);X22

11、 = 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五、 相同信噪比不同算法的比较clear all; clc; close all;N=1024; %长度Fs=500; %采样频率n=0:N-1;t=n/Fs; %时间序列a1=30; %信

12、号幅度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(t,x1,r);axis(-0.2 2 -40 40);hold on;plot(t,x2,:);axis(-0.2 2 -40 40);legend(x1信号, x2信号);xlabel(时间/s);y

13、label(x1(t) x2(t);title(原始信号);grid on;%互相关函数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)/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=

14、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);xlabel(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(t