互相关延时估计加权函数性能分析教材.docx
《互相关延时估计加权函数性能分析教材.docx》由会员分享，可在线阅读，更多相关《互相关延时估计加权函数性能分析教材.docx（16页珍藏版）》请在冰豆网上搜索。
互相关延时估计加权函数性能分析教材
互相关延时估计加权函数性能分析
广义互相关函数法是通过首先求出俩信号之间的互功率谱，然后在频域内给予一定的加权，以此对信号和噪音进行白化处理，从而增强信号中信噪比较高的频率成分，抑制噪声的影响，最后再反变换到时域，得到两信号之间的互相关函数，即：
（1）
其中
是广义互相关加权函数。
广义互相关加权函数的选择主要基于俩个方面：
噪声和反射情况。
根据不同的情况选择加权函数，其目的就是使
具有比较尖锐的峰值。
峰值处就是俩个传感器之间的时延。
由于来自同一声源的信号存在一定的相关性，通过计算不同麦克风所接受到的信号之间的相关函数，就可以估计出TDOA值。
然而在实际环境中，由于噪声和混响的影响，相关函数的最大峰会被弱化，有时还会出现多个峰值，这些都造成了实际峰值的检测困难。
此时就通过加权的方法来锐化峰值，通常我们通过时间、精度来确定算法的合理性。
一、广义互相关函数模拟
clearall;clc;closeall;
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.21.5-66]）;
holdon;
plot（t,x2,':
'）;
axis（[-0.21.5-66]）;
legend（'x1信号','x2信号'）;
xlabel（'时间/s'）;ylabel（'x1（t）x2（t）'）;
title（'原始信号'）;gridon;
holdoff
%互相关函数
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）'）;gridon
[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模拟
clearall;clc;closeall;
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.22-66]）;
holdon;
plot（t,x2,':
'）;
axis（[-0.22-66]）;
legend（'x1信号','x2信号'）;
xlabel（'时间/s'）;ylabel（'x1（t）x2（t）'）;
title（'原始信号'）;gridon;
holdoff
%互相关函数
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）'）;gridon
[max,location]=max（Cxy）;
%d=location-N/2-1
d=location-N
Delay=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=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./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-N
Delay=d/fs
运行程序的结果为：
d=
4
Delay=
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;%延迟点数
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=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-N
Delay=d/fs
运行程序的结果是：
d=
8
Delay=
0.0080
五、相同信噪比不同算法的比较
clearall;clc;closeall;
N=1024;%长度
Fs=500;%采样频率
n=0:
N-1;
t=n/Fs;%时间序列
a1=30;%信号幅度
a2=30;
d=9;%延迟点数
x1=a1*cos（2*pi*10*n/Fs）;%信号1
x1=awgn（x1,20）;%加噪声
%x1=x1.*hamming（max（size（x1）））';%加窗
x2=a2*cos（2*pi*10*（n+d）/Fs）;%信号2
x2=awgn（x2,20）;
%x2=x2.*hamming（max（size（x2）））';%加窗
subplot（511）;
plot（t,x1,'r'）;
axis（[-0.22-4040]）;
holdon;
plot（t,x2,':
'）;
axis（[-0.22-4040]）;
legend（'x1信号','x2信号'）;
xlabel（'时间/s'）;ylabel（'x1（t）x2（t）'）;
title（'原始信号'）;gridon;
%互相关函数
tic
X1=fft（x1,2*N-1）;
X2=fft（x2,2*N-1）;
Sxy=X1.*conj（X2）;
%Cxy=fftshift（ifft（Sxy））;
%GCC
Cxy=fftshift（ifft（Sxy））;
subplot（512）;
t1=（-N+1:
N-1）/Fs;
plot（t1,Cxy,'b'）;
title（'GCC'）;xlabel（'t/s'）;ylabel（'Cxy'）;gridon;
[max1,location1]=max（Cxy）;
%d=location-N/2-1
d1=location1-N
Delay1=d1/Fs%求得时间延迟
toc
%phat-gcc
tic
X1=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'）;gridon;
[max2,location2]=max（Pxy）;
%d=location-N/2-1
d2=location2-N
Delay2=d2/Fs%求得时间延迟
toc
%rhat-gcc
tic
X1=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
升级会员 