查看: 3267|回复: 0

麦克风阵列声源定位 SRP-PHAT（二）

[复制链接]

字体大小: 正常放大

浅夏110

542 主题	15 听众	1万积分

TA的每日心情

	开心 2020-11-14 17:15

签到天数: 74 天

[LV.6]常住居民II

群组: 2019美赛冲刺课程

群组: 站长地区赛培训

群组: 2019考研数学桃子老师

群组: 2018教师培训（呼伦贝

群组: 2019考研数学站长系列

电梯直达

1^#

发表于 2020-5-15 15:04 |只看该作者 |倒序浏览

|招呼Ta 关注Ta |邮箱已经成功绑定

DOA
　　声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm)，一种是基于空间谱估计如MUSIC等，还有就是基于beamforming的方法，也就是这里要介绍的可控波束响应（steered-response power），

steered-response power
　　可控波束响应是利用波束形成（beamforming）的方法，对空间不同方向的声音进行增强，得到声音信号最强的方向就被认为是声源的方向。
　　上一篇中简单介绍了麦克风阵列的背景知识，最简单的SRP就是利用延时-累加(delay-and-sum)的方法，寻找输出能量最大的方向。
　　其中，语音信号为宽带信号，因此需要做宽带波束形成，这里我们在频域实现

频域宽带波束形成
　　频域宽带波束形成可以归类为DFT波束形成器，结构如下图

频域处理也可以看做是子带处理（subband），DFT和IDFT的系数分别对应子带处理中的分析综合滤波器组，关于这一种解释，可参考《传感器阵列波束优化设计与应用》第六章。

频域宽带延时累加波束形成的基本过程就是信号分帧加窗->DFT->各频点相位补偿->IDFT 5 T$ B4 b8 q# @: E- Q; a6 j$ [. D, N
代码实现如下

nction [ DS, x1] = DelaySumURA( x,fs,N,frameLength,inc,r,angle)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%frequency-domain delay-sum beamformer using circular array
%
%    input :
%       x : input signal ,samples * channel
%       fs: sample rate
%       N : fft length,frequency bin number
%frameLength : frame length,usually same as N
%       inc : step increment
%       r : array element radius
%    angle : incident angle
%
%    output :
%       DS : delay-sum output
%       x1 : presteered signal,same size as x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c = 340;
Nele = size(x,2);
omega = zeros(frameLength,1);
H = ones(N/2+1,Nele);

theta = 90*pi/180; %固定一个俯仰角
gamma = [30 90 150 210 270 330]*pi/180;%麦克风位置
tao = r*sin(theta)*cos(angle(1)-gamma)/c;    %方位角 0 < angle <360
yds = zeros(length(x(:,1)),1);
x1 = zeros(size(x));

% frequency bin weights
% for k = 2:1:N/2+1
for k = 1:1:5000*N/fs
omega(k) = 2*pi*(k-1)*fs/N;
% steering vector
H(k,

= exp(-1j*omega(k)*tao);
end

for i = 1:inc:length(x(:,1))-frameLength

d = fft(bsxfun(@times, x(i:i+frameLength-1,

,hamming(frameLength)));

x_fft=bsxfun(@times, d(1:N/2+1,

,H);

% phase transformed
%x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,

));
yf = sum(x_fft,2);
Cf = [yf;conj(flipud(yf(2:N/2)))];

% 恢复延时累加的信号
yds(i:i+frameLength-1) = yds(i:i+frameLength-1)+(ifft(Cf));

% 恢复各路对齐后的信号
xf = [x_fft;conj(flipud(x_fft(2:N/2,

))];
x1(i:i+frameLength-1,

= x1(i:i+frameLength-1,

+(ifft(xf));
end
DS = yds/Nele;

end
然后遍历各个角度重复调用这个函数，测试实际录音数据，代码如下

%% SRP Estimate of Direction of Arrival at Microphone Array
% Frequency-domain delay-and-sum test
%
%%

% x = filter(Num,1,x0);
c = 340.0;

% XMOS circular microphone array radius
d = 0.0420;
% more test audio file in ../../TestAudio/ folder
path = '../../TestAudio/XMOS/room_mic5-2/';
[s1,fs] = audioread([path,'音轨-2.wav']);
s2 = audioread([path,'音轨-3.wav']);
s3 = audioread([path,'音轨-4.wav']);
s4 = audioread([path,'音轨-5.wav']);
s5 = audioread([path,'音轨-6.wav']);
s6 = audioread([path,'音轨-7.wav']);
signal = [s1,s2,s3,s4,s5,s6];
M = size(signal,2);
%%
t = 0;

% minimal searching grid
step = 1;

P = zeros(1,length(0:step:360-step));
tic
h = waitbar(0,'Please wait...');
for i = 0:step:360-step
% Delay-and-sum beamforming
[ DS, x1] = DelaySumURA(signal,fs,512,512,256,d,i/180*pi);
t = t+1;
%beamformed output energy
P(t) = DS'*DS;
waitbar(i / length(step:360-step))
end
toc
close(h)
[m,index] = max(P);
figure,plot(0:step:360-step,P/max(P))
ang = (index)*step

程序中用的是圆阵，可以进行二维方向角扫描，不过这里为了简便就固定了俯仰角，只扫描方位角，结果如下

结果与预期相同

PHAT加权

　　与GCC-PHAT方法相同，这里也可以对幅度做归一化，只保留相位信息，使得到的峰值更明显，提高在噪声及混响环境下的性能
5 f- k+ q) x% F2 x5 i　　上面代码中加上这一句

%x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,

));
测试同样的文件，结果如下

对比可以看到，PHAT加权的方法性能更好
参考
1.《SRP-PHAT-A High-Accuracy, Low-Latency Technique for Talker Localization in Reverberant Environments Using Microphone Arrays》
2. 《传感器阵列波束优化设计与应用》
————————————————
版权声明：本文为CSDN博主「373955482」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/u010592995/article/details/81586504

zan