互相關延時估計加權函數性能分析

? ? ? ?廣義互相關函數法是通過首先求出倆信號之間的互功率譜，然后在頻域內給予一定的加權，以此對信號和噪音進行白化處理，從而增強信號中信噪比較高的頻率成分，抑制噪聲的影響，最后再反變換到時域，得到兩信號之間的互相關函數，即：

（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時

從運行結果上來看，在時間上，基本互相關、PHAT加權、ROTH加權和SCOT加權四種算法的運行時間基本相同；但是從峰度的銳化來說，這四種方式的時延估計的準確性隨著信噪比的降低而惡化，互相關函數峰值的尖銳程度隨信噪比的降低而降低。對于SCOT加權來說，隨著信噪比的降低，性能急劇下降。基本互相關函數和RHOT加權雖然有一定的抗噪能力，但隨著信噪比的降低，其波動程度明顯加強，特別是對外圍的噪聲、反射和有限觀測數據很敏感，會造成峰值不明顯；對于PHAT加權，在較高的信噪比的時候，表現出了波動小、峰值尖銳的特性，在降低信噪比時，也表現出了較強的抗干擾性。

總結

以上是生活随笔為你收集整理的互相关延时估计加权函数性能分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。