閱讀原文還請移步我的知乎專欄：
https://www.zhihu.com/column/c_1287066237843951616

首先力薦鄢社鋒老師的書籍《優化陣列信號處理》，對于麥克風陣列信號處理的講解非常全面，對陣列的基本概念與優化方法講的很是透徹，而且最重要事情當然是全中文的啦，而且鄢社鋒是本領域的頂尖專家，比很多翻譯過來的外文書籍要好很多，值得閱讀，個人感覺比陳老師的《Microphone Array Signal Processing》要更容易理解、更容易上手很多，因為這本書里面包含了很多的案例用以對比演示，但是很遺憾，我并沒有找到鄢老師的相關代碼，所以干脆就自己動手，看書的同時順帶實作了幾個，

所以，本文是后續系列文章的第一篇，記錄案例實作的程序與代碼，對于案例的解釋主要為原書內容，自己進行了簡化，便于后續理解，本文講述第二章的案例2.2、2.3與2.4，對比試驗了常規波束形成與MVDR最優波束形成的波束圖、方位掃描圖與陣增益，以及自己實作的代碼，

例2.2 常規波束圖、常規波束掃描方位譜

考慮一個10元標準線列陣，假設一功率為1的單頻信號從 $80^\circ$ 方向入射到基陣，不考慮噪聲，用如下公式計算資料協方差矩陣：

$\bold R_x=E\left[ \bold x(n) \bold x^H(n) \right]$

以基陣在 $10^\circ$ 方向的回應向量作為導向向量，即 $\bar {\bold {a}}_s=\bold a\left( 80^\circ\right)$ ，采用式

$\bold w_c=\bar {\bold a}_s/M$

計算常規波束加權向量，用式

$p\left( \theta \right)=\bold w^H \bold a\left( \theta \right),\theta\in\Theta$

波束圖

計算波束回應，將得到的波束圖顯示于圖 2（a）中，圖中可以看見波束主瓣指向 $80^\circ$ 方位，這意味著該波束形成器對 $80^\circ$ 方向到達的信號回應最大，附代碼實作，在圖片上方，

%% 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) 常規波束形成波束圖

波束掃描方位譜

進行波束掃描，即讓 $\theta$ 在 $\left[ -90^\circ,90^\circ \right]$ 內變化，取導向向量為 $\bar {\bold {a}}_s=\bold a\left( \theta\right)$ ，用下式 $P\left( \theta \right)=\bold w^H\left( \theta \right)\bold R_x\bold w\left( \theta \right),\theta\in \Theta$

計算出波束掃描方位譜，顯示于圖 2（b）中，從圖中可以看出，在 $80^\circ$ 位存在一個主峰值，該峰值表示在該方位存在一個信號，峰值大小為 0dB，表示該方位波束輸出功率為1，即該方向到達信號功率估計值為1，同時注意到，雖然只有 $80^\circ$ 方向存在信號，其他方位沒有信號，但從方位譜上觀察，其他方位的波束輸出功率并不為0，好像是出現了“能量泄漏”，這是由于波束旁瓣引起的，對于遠離信號方向的波束，由于常規波束的旁瓣比較高，即使信號位于這些波束的旁瓣，旁瓣對該信號有一定的抑制，但是波束輸出并不為0，因此表現為在方位譜上非信號方向存在輸出功率，可以預見，由于常規波束較高的旁瓣，可能會因為強信號的“能量泄漏”而淹沒了其他方位的弱信號，這一點在后文的仿真中將會看到，

比較圖2（a）與圖2（b）發現，指向信號方向的常規波束圖與常規波束掃描方位譜曲線完全相同，這是在僅存在單位功率平面波信號時常規波束形成特有的現象，其他波束形成方法不存在此現象，假設基陣接收到功率為 0dB的空間白噪聲，信號功率增加到 30dB，即輸入信噪比 $SNR_{in}=30dB$ ，采用常規波束掃描得到的方位譜如圖 2（c）所示，圖中指示在 $80^\circ$ 方向方位譜大約為 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) 常規波束掃描方位圖

波束陣增益

應用公式： $G_c=\frac{\left| \bar{\bold a} ^H_s \bold a_s \right|^2}{\bar{\bold a} ^H_s \bold \rho _n\bold a_s}$

計算掃描到各方位時波束陣增益，顯示于圖 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 高斯白噪聲，即 $\bold R_n=\bold I$ ，一信噪比為 $SNR=30dB$ 的單頻信號從 $\theta_0=80^ \circ$ 方向入射到基陣，用下式計算資料協方差矩陣 $\bold R_x=E\left[ \bold x(n) \bold x^H(n) \right]$ ，

MVDR波束形成器波束圖

采用下式

$\bold w_{MVDR}=\frac{\bold R_x^{-1} \bar{\bold a}_s}{\bar{\bold a}_s^H\bold R_x^{-1} \bar{\bold a}_s}$

設計MVDR 波束形成器加權向量，圖 3（a）顯示了觀察方向分別為 $\theta_s=80^\circ$ 與 $\theta_s=100^\circ$ 的兩個波束圖，對于 $80^\circ$ 方向波束形成器，它相當于 $\beta=0$ 的情況， $100^\circ$ 方向平面波被視作干擾；對于 $100^\circ$ 方向波束， $100^\circ$ 方向平面波為期望信號，對應于 $\beta=1$ 的情況，

對于 $80^\circ$ 方向的波束圖，它在 $80^\circ$ 方向的回應為 $1(0dB)$ ，由于MVDR波束形成器的特性，為了使波束輸出功率最小，該波束在 $100^\circ$ 方向形成零點以最大限度抑制來自該方向的干擾，對于 $100^\circ$ 方向的波束圖，假想導向向量與真實基陣回應向量相等，又由于 $\bold R_n=\bold I$ ，該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波束掃描估計方位譜，其中掃描方位間隔為 $0.1^\circ$ ，得到的方位譜顯示于圖 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）顯示了各掃描方位波束的加權向量范數 $10lg\left( \left| \left| \bold w \right| \right|^2 \right)$ ，從圖中可以看出，在波束觀察方向距離信號方向較遠時（例如本例中與信號間隔 $10^\circ$ 以上）,波束形成器加權向量范數較小，接近于最小值，間隔較近時，隨著間隔減小，加權向量范數逐漸增加；而當波束方向恰好等于信號方向時，加權向量范數突然下降到與間隔較遠的方位同一大小（因為此時波束形成器蛻化為常規波束形成器），

%% 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('加權向量范數');