麦克风阵列声源定位 SRP-PHAT（二）
DOA声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm)，一种是基于空间谱估计如MUSIC等，还有就是基于beamforming的方法，也就是这里要介绍的可控波束响应（steered-response power），
steered-response power
可控波束响应是利用波束形成（beamforming）的方法，对空间不同方向的声音进行增强，得到声音信号最强的方向就被认为是声源的方向。
上一篇中简单介绍了麦克风阵列的背景知识，最简单的SRP就是利用延时-累加(delay-and-sum)的方法，寻找输出能量最大的方向。
其中，语音信号为宽带信号，因此需要做宽带波束形成，这里我们在频域实现
频域宽带波束形成
频域宽带波束形成可以归类为DFT波束形成器，结构如下图
https://img-blog.csdn.net/20180811150554894?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3UwMTA1OTI5OTU=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/00
频域处理也可以看做是子带处理（subband），DFT和IDFT的系数分别对应子带处理中的分析综合滤波器组，关于这一种解释，可参考《传感器阵列波束优化设计与应用》第六章。频域宽带延时累加波束形成的基本过程就是信号分帧加窗->DFT->各频点相位补偿->IDFT
代码实现如下
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 = *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 = ;
% 恢复延时累加的信号
yds(i:i+frameLength-1) = yds(i:i+frameLength-1)+(ifft(Cf));
% 恢复各路对齐后的信号
xf = ;
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/';
= audioread();
s2 = audioread();
s3 = audioread();
s4 = audioread();
s5 = audioread();
s6 = audioread();
signal = ;
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)
= max(P);
figure,plot(0:step:360-step,P/max(P))
ang = (index)*step
程序中用的是圆阵，可以进行二维方向角扫描，不过这里为了简便就固定了俯仰角，只扫描方位角，结果如下
https://img-blog.csdn.net/20180811154104879?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3UwMTA1OTI5OTU=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
结果与预期相同PHAT加权 与GCC-PHAT方法相同，这里也可以对幅度做归一化，只保留相位信息，使得到的峰值更明显，提高在噪声及混响环境下的性能
上面代码中加上这一句
%x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,:)));
测试同样的文件，结果如下
https://img-blog.csdn.net/20180811153854827?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3UwMTA1OTI5OTU=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/0
对比可以看到，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
页:
[1]