互相關延時估計權重函數性能分析
廣義互相關函數法是通過首先求出倆信号之間的互功率譜，然後在頻域内給予一定的權重，以此對信号和噪音進行白化處理，進而增強信号中信噪比較高的頻率成分，抑制噪聲的影響，最後再反變換到時域，得到兩信号之間的互相關函數，即：
（1）
其中是廣義互相關權重函數。廣義互相關權重函數的選擇主要基于倆個方面：噪聲和反射情況。根據不同的情況選擇權重函數，其目的就是使具有比較尖銳的峰值。峰值處就是倆個傳感器之間的時延。
由于來自同一聲源的信号存在一定的相關性，通過計算不同麥克風所接受到的信号之間的相關函數，就可以估計出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/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，而執行個體中給出的時延為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
- 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
- 相同信噪比不同算法的比較
clear all; clc; close all;
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.2 2 -40 40]);
hold on;
plot(t,x2,':');
axis([-0.2 2 -40 40]);
legend('x1信号', 'x2信号');
xlabel('時間/s');ylabel('x1(t) x2(t)');
title('原始信号');grid on;
%互相關函數
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');grid on;
[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');grid on;
[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(t1,Rxy,'b');
title('phat-gcc');xlabel('t/s');ylabel('Rxy');grid on;
[max3,location3]=max(Rxy);
%d=location-N/2-1
d3=location3-N
Delay3=d3/Fs %求得時間延遲
toc
%scot-gcc
tic
X1=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-gcc');xlabel('t/s');ylabel('SCxy');grid on;
[max4,location4]=max(SCxy);
%d=location-N/2-1
d4=location4-N
Delay4=d4/Fs %求得時間延遲
toc
- SNR=0時
- SNR=10時：
- SNR=20時
- SNR=50時