阅读原文还请移步我的知乎专栏:
https://www.zhihu.com/column/c_1287066237843951616
首先力荐鄢社锋老师的书籍《优化阵列信号处理》,对于麦克风阵列信号处理的讲解非常全面,对阵列的基本概念与优化方法讲的很是透彻,而且最重要事情当然是全中文的啦,而且鄢社锋是本领域的顶尖专家,比很多翻译过来的外文书籍要好很多,值得阅读。个人感觉比陈老师的《Microphone Array Signal Processing》要更容易理解、更容易上手很多,因为这本书里面包含了很多的案例用以对比演示。但是很遗憾,我并没有找到鄢老师的相关代码,所以干脆就自己动手,看书的同时顺带实现了几个。
所以,本文是后续系列文章的第一篇,记录案例实现的过程与代码,对于案例的解释主要为原书内容,自己进行了简化,便于后续理解。本文讲述第二章的案例2.2、2.3与2.4,对比试验了常规波束形成与MVDR最优波束形成的波束图、方位扫描图与阵增益,以及自己实现的代码。
全部实现代码下载位置:
鄢社锋老师的书籍《优化阵列信号处理》前三章关键案例Matlab实现代码-专业指导文档类资源-CSDN下载
例2.2 常规波束图、常规波束扫描方位谱
考虑一个10元标准线列阵,假设一功率为1的单频信号从 方向入射到基阵,不考虑噪声,用如下公式计算数据协方差矩阵:
以基阵在 方向的响应向量作为导向向量,即 ,采用式
计算常规波束加权向量,用式
波束图
计算波束响应。将得到的波束图显示于图 2(a)中。图中可以看见波束主瓣指向方位。这意味着该波束形成器对方向到达的信号响应最大。附代码实现,在图片上方。
%% 2(a) 常规波束形成波束图 代码实现部分
clear; close all; clc;
c=340; %声速
f=1000; %频率
lambda = c / f;
k = 2*pi*f / c;
space = lambda/2; %麦克风间距
M = 10; %麦克风数量
theta_d = 80*pi/180; %入射角度
theta_angle = 0:0.1:180;
theta = theta_angle*pi/180;
Wc = exp(-1i*k*space*[0:M-1]*cos(theta_d)) / M;
B = zeros(size(theta));
%% 计算波束图
for i = 1:length(theta)
p = exp(1i*k*space*[0:M-1]*cos(theta(i)));
B(i) = conj(Wc)*p';
end
B_db = 20*log10(B);
limit_dB = -50;
index = B_db < limit_dB;
B_db(index) = limit_dB;
figure;
plot(theta_angle, B_db, 'linewidth', 1.5);
grid on;xlabel('方位(^o)'); ylabel('20lg(B)/dB');
title('波束图');
2(a) 常规波束形成波束图
波束扫描方位谱
进行波束扫描,即让 在 内变化,取导向向量为 。用下式
计算出波束扫描方位谱,显示于图 2(b)中。从图中可以看出,在 位存在一个主峰值,该峰值表示在该方位存在一个信号。峰值大小为 0dB,表示该方位波束输出功率为1,即该方向到达信号功率估计值为1。同时注意到,虽然只有方向存在信号,其他方位没有信号,但从方位谱上观察,其他方位的波束输出功率并不为0,好像是出现了“能量泄漏”。这是由于波束旁瓣引起的,对于远离信号方向的波束,由于常规波束的旁瓣比较高,即使信号位于这些波束的旁瓣,旁瓣对该信号有一定的抑制,但是波束输出并不为0,因此表现为在方位谱上非信号方向存在输出功率。可以预见,由于常规波束较高的旁瓣,可能会因为强信号的“能量泄漏”而淹没了其他方位的弱信号,这一点在后文的仿真中将会看到。
比较图2(a)与图2(b)发现,指向信号方向的常规波束图与常规波束扫描方位谱曲线完全相同。这是在仅存在单位功率平面波信号时常规波束形成特有的现象,其他波束形成方法不存在此现象。假设基阵接收到功率为 0dB的空间白噪声,信号功率增加到 30dB,即输入信噪比 。采用常规波束扫描得到的方位谱如图 2(c)所示,图中指示在 方向方位谱大约为 30dB,表示该方位信号功率估计值约为 30dB。代码实现参考 2(b),不再赘述。
%% 2(b) 常规波束扫描方位谱,代码实现部分
p = exp(1i*k*space*[0:M-1]*cos(theta_d));
Rx = (p')*conj(p);
B = zeros(size(theta));
for i = 1:length(theta)
Wcc = exp(-1i*k*space*[0:M-1]*cos(theta(i))) / M;
B(i) = conj(Wcc)*Rx*Wcc';
end
B_db = 10*log10(B);
limit_dB = -50;
index = B_db < limit_dB;
B_db(index) = limit_dB;
figure;
plot(theta_angle, B_db, 'linewidth', 1.5);
grid on;
xlabel('方位(^o)'); ylabel('10lg(P)/dB');
title('波束扫描方位谱');
2(b) 常规波束扫描方位谱
2(c) 常规波束扫描方位图
波束阵增益
应用公式:
计算扫描到各方位时波束阵增益,显示于图 2(d)中。图中可见,该曲线形状与图2(a)中波束形状相同,但高出10dB ,即为最高阵增益。从该图可以知,当波束对准信号方向时,波束形成器的阵增益最大;波束观察方位与信号方位存在偏差时,波束阵增益减小。只要方位偏差不太大(在半功率束宽以内),阵增益比最大值下降得并不多。
%% 2(d) 常规波束阵增益 代码实现部分
p_ = exp(1i*k*space*[0:M-1]*cos(theta_d));
B = zeros(size(theta));
for i = 1:length(theta)
p = exp(1i*k*space*[0:M-1]*cos(theta(i)));
B(i) = (conj(p_) * p').^2 / (conj(p_)*p_');
end
B_db = 10*log10(B);
limit_dB = -50;
index = B_db < limit_dB;
B_db(index) = limit_dB;
figure;
plot(theta_angle, B_db, 'linewidth', 1.5);
grid on;
xlabel('方位(^o)'); ylabel('10lg(G)/dB');
title('不同方向波束阵增益');
2(d) 常规波束阵增益
例2.3 MVDR波束形成
考虑一个10元标准线列阵,假设基阵接收噪声是功率为0dB 高斯白噪声,即 。一信噪比为 的单频信号从 方向入射到基阵,用下式计算数据协方差矩阵 。
MVDR波束形成器波束图
采用下式
设计MVDR 波束形成器加权向量。图 3(a)显示了观察方向分别为与 的两个波束图。对于 方向波束形成器,它相当于 的情况, 方向平面波被视作干扰;对于方向波束,方向平面波为期望信号,对应于 的情况。
对于方向的波束图,它在方向的响应为 。由于MVDR波束形成器的特性,为了使波束输出功率最小,该波束在 方向形成零点以最大限度抑制来自该方向的干扰。对于方向的波束图,假想导向向量与真实基阵响应向量相等,又由于,该MVDR波束形成器蜕化为常规波束形成器,所以该波束图与图 2(a)中的常规波束图相同。
%% 3(a) MVDR波束形成器波束图,代码实现
c=340; %声速
f=1000; %频率
lambda = c / f;
k = 2*pi*f / c;
space = lambda/2; %麦克风间距
M = 10; %麦克风数量
theta_angle = 0:0.1:180;
theta = theta_angle*pi/180;
theta_d1 = 80*pi/180; %入射角度
p_1 = exp(1i*k*space*[0:M-1]*cos(theta_d1))';
Rx = conj(p_1)*p_1'*10^(30/10) + eye(M)*10^(0/10);
Wc = Rx\p_1 / (conj(p_1')*inv(Rx)*p_1);
B = zeros(size(theta));
%% 计算波束图 80度
for i = 1:length(theta)
p = exp(1i*k*space*[0:M-1]*cos(theta(i)))';
B(i) = conj(Wc')*p;
end
B_db = 20*log10(B);
limit_dB = -50;
index = B_db < limit_dB;
B_db(index) = limit_dB;
figure;
plot(theta_angle, B_db, 'linewidth', 1.5);
hold on;
3(a) MVDR波束形成器波束图
MVDR波束扫描方位谱
采用MVDR波束扫描估计方位谱,其中扫描方位间隔为 ,得到的方位谱显示于图 3(b)中。与图 2(b)相比较可见,MVDR波束扫描得到的方位谱峰值更尖锐。可见,MVDR波束扫描能提供比常规波束扫描较高的方位分辨能力。
%% 3(b) MVDR波束形成器波束扫描方位谱图,代码实现
theta_angle = 0:1:180;
theta = theta_angle*pi/180;
theta_d1 = 80*pi/180; %入射角度1
p_1 = exp(1i*k*space*[0:M-1]*cos(theta_d1))';
Rx = (p_1)*conj(p_1')*10^(30/10);
B = zeros(size(theta));
for i = 1:length(theta)
p = exp(1i*k*space*[0:M-1]*cos(theta(i)))';
Wcc = Rx\p / (conj(p')*(Rx\p));
B(i) = conj(Wcc')*Rx*Wcc;
end
B_db = 10*log10(B);
limit_dB = -10;
index = B_db < limit_dB;
B_db(index) = limit_dB;
figure;
plot(theta_angle, B_db, 'linewidth', 1.5);
hold on;
plot(80, 30, 'ro', 'linewidth', 1.5);
grid on;
legend('方位谱', '真实信号');
xlabel('方位(^o)'); ylabel('10lg(P)/dB');
title('波束扫描方位谱');
3(b) MVDR波束形成器波束扫描方位谱图
MVDR波束形成器加权向量范数
图 3(v)显示了各扫描方位波束的加权向量范数 。从图中可以看出,在波束观察方向距离信号方向较远时(例如本例中与信号间隔 以上),波束形成器加权向量范数较小,接近于最小值。间隔较近时,随着间隔减小,加权向量范数逐渐增加;而当波束方向恰好等于信号方向时,加权向量范数突然下降到与间隔较远的方位同一大小(因为此时波束形成器蜕化为常规波束形成器)。
%% 3(c) MVDR波束形成器加权向量范数,代码实现
theta_angle = 0:0.1:180;
theta = theta_angle*pi/180;
theta_d1 = 90*pi/180; %入射角度1
p_1 = exp(1i*k*space*[0:M-1]*cos(theta_d1))';
Rx = (p_1)*conj(p_1')*10^(30/10);
B = zeros(size(theta));
for i = 1:length(theta)
p = exp(1i*k*space*[0:M-1]*cos(theta(i)))';
Wcc = Rx\p / (conj(p')*(Rx\p));
B(i) = norm(Wcc)*norm(Wcc);
end
B_db = 10*log10(B);
limit_dB = -50;
index = B_db < limit_dB;
B_db(index) = limit_dB;
figure;
plot(theta_angle, B_db, 'linewidth', 1.5);
grid on;
xlabel('方位(^o)'); ylabel('10lg(w^2)/dB');
title('加权向量范数');
3(c) MVDR波束形成器加权向量范数
例2.4 两信号源情况下MVDR与常规波束形成器的比较
考虑一个 10元标准线列阵,基阵接收噪声是功率为 0dB高斯白噪声,两个信噪比分别为 15dB与 30dB的单频信号分别从 与 方向入射到基阵。
分别采用常规波束形成与MVDR波束形成扫描,得到的方位谱分别如图4(a) 与图4(b) 所示。
4(a) 两信号源常规波束形成器
4(b) 两信号源MVDR波束形成器
参考资料
1、《优化阵列信号处理》,鄢社锋