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

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

广义互相关函数法是通过首先求出俩信号之间的互功率谱，然后在频域给予一定的加权，以此对信号和噪音进行白化处理，从而增强信号噪比较高的频率成分，抑制噪声的影响，最后再反变换到时域，得到两信号之间的互相关函数，即：

⎰=π

τψτ0121212)()()(dw e w G w R jw （1） 其中)(12w ψ是广义互相关加权函数。广义互相关加权函数的选择主要基于俩个方面：噪声和反射情况。根据不同的情况选择加权函数，其目的就是使)(12τR 具有比较尖锐的峰值。)(12τR 峰值处就是俩个传感器之间的时延。

由于来自同一声源的信号存在一定的相关性，通过计算不同麦克风所接受到的信号之间的相关函数，就可以估计出TDOA 值。然而在实际环境中，由于噪声和混响的影响，相关函数的最大峰会被弱化，有时还会出现多个峰值，这些都造成了实际峰值的检测困难。此时就通过加权的方法来锐化峰值，通常我们通过时间、精度来确定算法的合理性。

一、 广义互相关函数模拟

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); %信号1

x1=x1+randn(size(x1)); %加噪声

x2=a2*cos(2*pi*10*(n+d)/Fs); %信号2

x2=x2+randn(size(x2));

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)/Fs; %注意

plot(t1,Cxy,'b');

title('互相关函数');xlabel('时间/s');ylabel('Rx1x2(t)');grid on [max,location]=max(Cxy);%求出最大值max,及最大值所在的位置（第几行）location;

%d=location-N/2-1 %算出延迟了几个点

d=location-N

Delay=d/Fs %求得时间延迟

运行程序得到的结果是：

d =

2

Delay =

0.0040

可以看出，通过互相关函数的求解d=2，delay=0.0040，这和我们给出的信号的时延d/Fs=0.0040是一致的。这表明互相关函数可以给出信号的时延估计。

二、PHAT-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); %信号1

x1=x1+randn(size(x1)); %加噪声

%x1=x1.*hamming(max(size(x1)))';%加窗

x2=a2*cos(2*pi*10*(n+d)/Fs); %信号2

x2=x2+randn(size(x2));

%x2=x2.*hamming(max(size(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=(-N+1:N-1)/Fs;

plot(t1,Cxy,'b');

title('Rx1x2');xlabel('t/s');ylabel('Rx1x2(t)');grid on [max,location]=max(Cxy);

%d=location-N/2-1

d=location-N

Delay=d/Fs %求得时间延迟

运行程序得到的结果是：

d =

1

Delay =

0.0020