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

Posted 2022-03-21 fpga&matlab

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

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

（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

1. SNR=0时

1. SNR=10时：

1. SNR=20时

1. SNR=50时

从运行结果上来看，在时间上，基本互相关、PHAT加权、ROTH加权和SCOT加权四种算法的运行时间基本相同；但是从峰度的锐化来说，这四种方式的时延估计的准确性随着信噪比的降低而恶化，互相关函数峰值的尖锐程度随信噪比的降低而降低。对于SCOT加权来说，随着信噪比的降低，性能急剧下降。基本互相关函数和RHOT加权虽然有一定的抗噪能力，但随着信噪比的降低，其波动程度明显加强，特别是对外围的噪声、反射和有限观测数据很敏感，会造成峰值不明显；对于PHAT加权，在较高的信噪比的时候，表现出了波动小、峰值尖锐的特性，在降低信噪比时，也表现出了较强的抗干扰性。