当前位置：首页 > news >正文

信号处理之声源定位

news 来源：原创 2024/11/18 10:41:05

1. 什么是声源定位？

本文讨论，通过麦克风阵列，判断声音的方位（距离、方位角azimuth和俯仰角elevation）。

2. 声源定位算法

远场模型定义如下图：

信号源s到观察信号y1和y2的时间差为

可以采用如下方法计算时间差 $\tau_{12}$

（1）互相关方法（cross correlation）

计算观察信号 $y_1(k)$ 和 $y_2(k)$ 之间的互相关函数。

互相关函数取最大值时的p值，即对应的时延，此时 $y_1(k)$ 和 $y_2(k+p)$ 重叠。

（2）广义互相关方法（generalized cross correlation）

由于噪声和混响的存在，导致互相关函数的峰值不明显，导致估计不准，考虑采用广义互相关方法。即先将时域转换为频域，在频域进行归一化操作，达到降噪的效果，再傅里叶逆变换至时域。

根据Wiener-Khichine定理，互功率谱(一般用大写字母 $\Phi$ 或 $P$ 表示)

$\Phi _{y_1y_2}(f)=Fourier(r_{y_1y_2}(p))$

$r_{y_1y_2}(p)=Fourier^{-1}(\Phi {y_1y_2}(f))$

这里可以对功率谱进行归一化

$r_{y_1y_2}(p)=Fourier^{-1}(\frac{1}{|\Phi {y_1y_2}(f)|}\Phi {y_1y_2}(f))$

接下来同CC方法一样

附上Matlab代码

%% 利用基于时延估计的方法（CC，GCC，GCC-PHAT）
% 作者：https://blog.csdn.net/pk296256948/article/details/116294055
clc
clear
close all

% 加载matlab自带的敲锣声
load gong;
% 采样频率
Fs = 8192;
% 采样周期
dt = 1/Fs;
% 定义声源
music_src = y;

% 设置两个麦克风的坐标
mic_d = 2;
mic_x = [-mic_d mic_d];
mic_y = [0 0];
plot(mic_x,mic_y,'x'); % 观察两麦克风的位置
grid on
axis([-5 5 -5 5])
hold on;
quiver(-5,0,10,0,1,'color','black'); % 画x,y轴
quiver(0,-5,0,10,1,'color','black');
xlabel('X')
ylabel('Y')
title("声源到两麦克风距离");
% 设置声源位置
s_x = 20;
s_y = 10;
plot(s_x,s_y,'o');
quiver(s_x,s_y,-s_x-mic_d,-s_y,1); % 假设两麦克风接收声信号路线
quiver(s_x,s_y,-s_x+mic_d,-s_y,1);

% 求出两麦克风分别到声源的距离
dis_s1 = sqrt((mic_x(1)-s_x).^2+(mic_y(1)-s_y).^2); % 麦克风1到声源距离
dis_s2 = sqrt((mic_x(2)-s_x).^2+(mic_y(2)-s_y).^2); % 麦克风2到声源距离
c = 343; % 设置声速
delay = abs((dis_s1-dis_s2)./c); % 计算两麦克风之间的时延

% 设置延时
music_delay = delayseq(music_src,delay,Fs);
figure(2);
subplot(2,1,1);
plot(music_src);
axis([0 length(music_src) -2 2]);
subplot(2,1,2);
plot(music_delay);
axis([0 length(music_delay) -2 2]);

% 使用matlab自带的gcc-phat算法
[tau,R,lag] = gccphat(music_delay,music_src,Fs);
disp("使用gcc-phat时延为");
disp(tau);
figure(3);
t=1:length(tau);
plot(lag,real(R(:,1)));

% 使用cc算法
[rcc,lag] = xcorr(music_delay,music_src);
figure(4);
plot(lag/Fs,rcc);
[M,I] = max(abs(rcc));
lagDiff = lag(I);
timeDiff = lagDiff/Fs;
disp("使用自相关的时延为");
disp(timeDiff);

% 使用公式gcc-phat算法
RGCC = fft(rcc);
rgcc = ifft(RGCC*1./abs(RGCC));
figure(5);
plot(lag/Fs,rgcc);
[M,I] = max(abs(rgcc));
lagDiff = lag(I);
timeDiff = lagDiff/Fs;
disp("使用gcc-phat的时延为");
disp(timeDiff);

% 计算角度，这里假设为平面波
dis_r = tau*c;
angel = acos(tau*c./(mic_d*2))*180/pi;
if dis_s1<dis_s2
    angel = 180-angel;
end
disp("计算声源的角度为");
disp(angel);

参考资料

[1] Generalized cross-correlation - MATLAB gccphat- MathWorks 中国

[2] 声源定位之GCC-PHAT算法_抽屉疯了的博客-CSDN博客

[3] 声源定位 - 知乎