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

DOA

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

steered-response power

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

频域宽带波束形成

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

频域处理也可以看做是子带处理（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 = [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/fsomega(k) = 2*pi*(k-1)*fs/N;   % steering vectorH(k,:) = exp(-1j*omega(k)*tao);
endfor i = 1:inc:length(x(:,1))-frameLengthd = 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 energyP(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方法相同，这里也可以对幅度做归一化，只保留相位信息，使得到的峰值更明显，提高在噪声及混响环境下的性能
　　上面代码中加上这一句
　　

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

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

对比可以看到，PHAT加权的方法性能更好

代码及测试文件在github

参考

1.《SRP-PHAT-A High-Accuracy, Low-Latency Technique for Talker Localization in Reverberant Environments Using Microphone Arrays》
2. 《传感器阵列波束优化设计与应用》