fixed wideband beamformer

最新推荐文章于 2021-03-18 10:44:26 发布

373955482

最新推荐文章于 2021-03-18 10:44:26 发布

阅读量1.5k

点赞数 2

分类专栏： Microphone Array 数字信号处理音频处理文章标签：麦克风阵列 superdirective beamformer wideband beamformer WGN beamforming

本文链接：https://blog.csdn.net/u010592995/article/details/81699912

版权

数字信号处理同时被 3 个专栏收录

29 篇文章 45 订阅

订阅专栏

音频处理

19 篇文章 19 订阅

订阅专栏

Microphone Array

12 篇文章 36 订阅

订阅专栏

1. fixed beamformer

　　根据滤波器系数是否独立于输入数据（data-independent）,beamformer可以分为固定和自适应两大类，这里先介绍一些固定宽带波束形成方法

2. Evaluation of Beamfomers

2.1 Array-Gain

　　array-gain反映了阵列对输入信号SNR改善程度，定义为阵列输出SNR与输入SNR（参考阵元SNR）之比1，当然，计算时都是用dB单位（ $10log_{10}(*)$ ），这里先省略（for the sake of brevity）

2.2 Directivity

　　方向因子DI（directivity index）通常用来衡量阵列对散射噪声的抑制能力

2.3 WGN

　　与上面相对应的，WGN是描述阵列对空间不相关噪声的抑制能力，像一些环境白噪声、元器件热噪声、以及一些由于DOA误差阵元不匹配引入的误差都可归类为空间白噪声，从这点来看，WGN可以作为用来衡量阵列系统稳定性的指标
　　WGN针对的是空间不相关噪声，因此噪声的互谱密度矩阵为单位矩阵，即 $\Gamma = \mathbf{I}$ ，代入到阵列增益定义式中得

3 superdirective beamfomer

超指向波束形成目的是为了最大化方向因子或阵列增益，在约束目标方向不变的条件下，得到MVDR解形式如下

3.1 performance analysis

以5麦间距5cm阵列，endfire方向为例，与delay-and-sum 对比

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
%   superdirective beamformer 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
% clear all;
close all;
imag=sqrt(-1);
element_num=5;%阵元数
% d_lamda=1/2;%阵元间距d与波长lamda的关系
theta=linspace(-pi/2,pi/2,180);    %扫描范围
theta0=90;                    %注视方向
theta0 = theta0*pi/180;
% w=exp(-1j*2*pi*d_lamda*sin(theta0)*[0:element_num-1]');
k = 64;
y1 = zeros(128,length(theta));
omega = zeros(1,128);
fs = 8000;
N = 512;
d = 0.05;
c = 340;

f = (0:1:N/2+1-1)*fs/N;
f(1) = (1e-14)*fs/N;
Nele = element_num;

% diffuse noise field MSC
Fvv = zeros(N/2+1,Nele,Nele);
for i = 1:Nele
    for j = 1:Nele   
        if i == j
            Fvv(:,i,j) = ones(1,N/2+1);
        else
            dij = abs(i-j)*d;
            Fvv(:,i,j) = sin(2*pi*f*dij*1.01/c)./(2*pi*f*dij*1.01/c);%T(1) = 0.999;%T(2) = 0.996;
        end
    end
end

H = zeros(Nele,N/2+1);
% H = zeros(N/2+1,Nele);
DI_sdb = zeros(N/2+1,1);  % superdirective directive factor
DI_dsb = zeros(N/2+1,1);  % DS directive factor

WGN_sdb = zeros(N/2+1,1);  % superdirective directive factor
WGN_dsb = zeros(N/2+1,1);  % DS directive factor

y2 = zeros(N/2,length(theta));
for k = 1:N/2+1
% for k = 32:32
for  j=1:length(theta)
    omega(k) = 2*pi*(k-1)*fs/N;
    w0=exp(-1j*omega(k)*d*sin(theta0)/c*[0:element_num-1]');
    a =exp(-1j*omega(k)*d*sin(theta(j))/c*[0:element_num-1]');
    y1(k,j)=w0'*a;
    DI_dsb(k) = (abs(w0'*w0))^2 ...
                /(w0'*squeeze(Fvv(k,:,:))*w0);
    WGN_dsb(k) =  (abs(w0'*w0))^2 ...
                /(w0'*w0);       

    % MVDR soulution
    Fvv_k = (squeeze(Fvv(k,:,:))+1e-14*eye(Nele));
    H(:,k) =    Fvv_k\w0 ...
                 ./(w0'/Fvv_k*w0);

    DI_sdb(k) = (abs(H(:,k)'*w0))^2 ...
                /(H(:,k)'*squeeze(Fvv(k,:,:))*H(:,k));

    WGN_sdb(k) = (abs(H(:,k)'*w0))^2 ...
                /(H(:,k)'*H(:,k));

    y2(k,j)=H(:,k)'*(a);

end
end

figure,mesh(theta*180/pi,(0:1:N/2+1-1)*fs/N,pow2db(abs(y1)/max(max(abs(y1)))))
zlim([-40 0]);
figure,mesh(theta*180/pi,(0:1:N/2+1-1)*fs/N,pow2db(abs(y2)/max(max(abs(y2)))))
zlim([-40 0]);

figure,plot(f,pow2db(abs(DI_dsb))),ylim([0 15]),grid on
hold on,plot(f,pow2db(abs(DI_sdb))),ylim([0 15]),grid on
legend('delaysum','superdirective'),title('Directivity Index')

figure,plot(f,pow2db(WGN_dsb)),grid on
hold on,plot(f,pow2db(WGN_sdb))
legend('delaysum','superdirective'),title('WGN')

　　先看波束图对比，从左图可以看出，delay-sum在低频的主瓣非常宽，几乎在700Hz以下就没有方向性了，而对比右图，superdirective beamformer的低频性能就好很多

再看方向因子

　　跟上面的分析相同，superdirective beamformer在低频段有更好的有更好的指向性，其中指向性最高点处可达到 $M^2$ ,而在高频段，DS与superdirctive指向性几乎相同

再看WGN

　　由上图可以看出，DS对空间不相关噪声在所有频率范围内为常数 $（10log(M)）$ ，而SU虽然在低频段指向性很好，但是WGN却很差，也就是说此时SU会显著放大噪声，从这里可以看出，DS的稳定性很好，而MVDR就很容易受到误差的影响。

4 simulation

　　现在就用实际数据测试看看效果，用XMOS 6+1 阵列板录音，先用srp-phat计算得到角度，再经过beamformer处理
测试代码如下

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)),ylim([0.86 1])
ang = (index)*step

% ang = 210;

[ DS0, x1] = DelaySumURA(signal,fs,1024,1024,512,d,ang/180*pi);
% DS = filter(HP_Num,1,DS0);
% audiowrite([path,'xmos-DS0.wav'],real(DS0),fs)
%% diffuse noise field MSC
N = 256;

f = (1:N/2+1)*fs/N;
% f(1) = 1e-8;
Fvv = zeros(N/2+1,M,M);
for i = 1:M
    for j = 1:M   
        if i == j
            Fvv(:,i,j) = ones(1,N/2+1);
        else
            if(abs(i-j)==1)
                dij = d;
            elseif(abs(i-j)==2)
                dij = d*sqrt(3);
            elseif(abs(i-j)==3)
                dij=d*2;
            end
            Fvv(:,i,j) = sin(2*pi*f*dij*1/c)./(2*pi*f*dij*1/c);%T(1) = 0.999;%T(2) = 0.996;
        end
    end
end
[b,a] = sos2tf(SOS);
signal0 = filter(b,a,signal);
[ super1, x1,~,DI] = superdirectiveMVDR(signal0,fs,N,N,N/2,d,ang/180*pi,Fvv);
% audiowrite([path,'xmos-super12.wav'],super1,fs)