最近学习了一下《最优阵列处理技术》，应老师要求写一个线性均匀水听器阵列的常规波束形成，由于是初学者，写的可能会有点问题，欢迎大家提出修改建议和指导，写这个主要是记录自己的思考，其次是和初学者进行交流提升。

1. 主要原理阐述

要实现波束形成，首先得了解频率波束响应，波束形成与之息息相关。

1.1 频率波数响应

对于一个水听器阵列，我们要想的得到其对外部信号场的响应，则需要对信号场在空域上进行采样，得到一组信号，表示为：

$f(t,p)=\begin{bmatrix} f(t-p_1)\\ f(t-p_2)\\ ... \\ f(t-p_n) \end{bmatrix}$

对每个阵元的输出用一个线性时不变滤波器进行处理，该滤波器的冲激响应为 $h_n(\tau)$ ，并对所有输出求和，得到阵列的输出 $y(t)$ ，该过程如下图示：

$y(t)=\int h_n(t-\tau)f(\tau,p)d\tau$

针对线性阵列，我们用一张图来说明一下：

设均匀阵元排列间隔为 $d$ ，假如信号平行射入，信号入射角度与端射方向呈 $\theta$ 角，则可以通过计算知道，相邻的阵元接收到信号的间隔为:

$\tau = \frac{dcos\theta}{c}$

设第一个阵元接收到的信号为 $sig(t)$ ，则第二个阵元接收到的信号为 $sig(t+\tau)$ ，以此类推，在理想情况下，第n个阵元接收到的信号为：

$sig(t+(n-1)\tau),(n = 1,2,3,..,N)$

为了能够还原我们想要的信号，我们可以使用一个延时求和波束形成器，其原理如下：

式子表示为：

$h_n(\tau) = \delta(\tau+\tau_n)/N$

也可以写成频域内的简洁形式：

$H_T(\omega)=v_{k}^{H}(k)/N$

其中， $v_k (k)$ 称作流形矢量，定义为：

$v_k(k)= \begin{bmatrix} e^{-jk^TP_0}& e^{-jk^TP_1}& ...& e^{-jk^TP_n} \end{bmatrix}^T$

$v_{k}^H$ 右上角的H为共轭转置的意思。

其中：

$k^TP_n = \omega \tau _n$

$\tau_n = -\frac{\mu ^TP_n}{c}$

$\mu ^T = \begin{bmatrix} sin\theta cos\phi\\sin\theta sin\phi \\ cos\theta \end{bmatrix}$

$P_n = (n-1)d, (n=1,2,3,...,N)$

$\omega$ 为信号的频率， $\mu$ 为方向余弦， $c$ 为介质中的信号传播速度。

为了进一步讲解，我们假设一个基函数，表示如下：

$f(t,p) = e^{j{\omega}t}v_k(k)$

则阵列处理器对一个平面波的响应：

$y(t)=\int h_n(t-\tau)f(\tau,p)d\tau$

可以表示为：

$y(t) = H^T(\omega)v_k(k)e^{jwt}$

在频域内，可以写为：

$Y(w)=H^T(w)v_k(k)$

定义：

$\Upsilon(w) =H^T(w)v_k(k)$

为阵列的频率波束响应。

波束方向图的定义是用入射方向表示频率波数响应函数（在代码中我们使用的角度为60°，可以自行修改）。

2. 代码讲解

2.1 基函数

在代码中，我们沿用基函数：

$f(t,p) = e^{j{\omega}t}v_k(k)$

$v_k(k)= \begin{bmatrix} e^{-jk^TP_0}& e^{-jk^TP_1}& ...& e^{-jk^TP_n} \end{bmatrix}^T$

在我们进行频率波数响应推导的时候，细心的朋友应该发现了，信号接收时间间隔 $\tau$ 的定义即为 $\tau_n$ 的定义， $P_n$ 为位置表示 $(n-1)d$ ，由于是在处理一个平面波信号，方向余弦 $\mu$ 仅需取一个方向，这里为 $cos\theta$ ，故根据 $\tau_n$ 的定义：

$\tau = \frac{dcos\theta}{c}$

$\tau_n =(n-1) \tau = \frac{(n-1)dcos\theta}{c},(n=1,2,3,...,N)$

假设我们的信号为 $s = e^{j\omega_0 t}$ ，其中， $\omega_0$ 为信号声源的频率，那么我们现在可以表示每个阵元的接收信号：

$s_n=e^{j\omega(t+(n-1)\tau)},(n=1,2,3,...,N)$

打开有：

$s_n=e^{j\omega t} * e^{j\omega\tau_n}=f(t,p) = e^{j{\omega}t}v_k(k)$

该部分代码为：

%% 目标信号生成
f0=300;               % 声源频率，Hz
s=exp(1j*2*pi*f0.*t); % 目标复数信号，1 x LEN  "1j为复数，1为数字一，j为负数"  
as=exp(1j*2*pi*ra*f0*cos(thetas_rad)/c); % 目标驾驶向量，N x 1，驾驶向量=exp(-2*pi*f0*τ*1j*)
snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的水听器输出时间波形实数矩阵，Add white Gaussian noise to signal
x = rx.';          % 输出波形矩阵转置，便于二维数组FFT列运算，LEN x N，即接收到的信号

s为源信号，as为流形矢量（也有称驾驶向量的说法），x为得到的输出信号。

注意：正常来讲，我们在进行信号接收时，应该在一定频率内进行扫描，获取不同频率信号再进行筛选；在代码中，我们假设的信号实际是简化了，相当于我们本来就知道信号的频率（设为300Hz），现在只需要在该频率下进行角度扫描，确定信号的来源。

2.2 延时处理

接着，我们刚刚提到了延时求和波束形成器 $H^T(w)$ ，根据频率响应函数的定义，我们需要对信号基函数 $f(t,p) = e^{j{\omega}t}v_k(k)$ 进行延时求和，由此得到波束方向图：

$y(t) = H^T(\omega)v_k(k)e^{jwt}$

我们在1.1频率波数响应函数末尾提到，在代码中我们已经确定了搜索的信号频率，因此仅需要对角度进行扫描，综上所述，代码如下：

%% 目标方位角设置
thetas_deg=60;                 % 目标方位角，度deg
thetas_rad=thetas_deg*deg2rad; % 目标方位角，弧度rad
theta = 0:0.1:180;                  % 任意信号入射角度%% 频域常规波束形成
tic
for i = 1:length(theta)tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延pt_sig = exp(1j*2*pi*f0*tao)' * x.';  %波束，1*LENs_sig(i,:) = pt_sig/N;    %size(s_sig) = i*LEN
end
toc

我们的角度扫描从0到180°，步长为0.1，因此可以得到1801长度的theta，tao为时延补偿，对应 $h_{n}(\tau)$ 中各个阵元的延时补偿。

注意：

1. 频率波束响应公式中 $H^T(w)$ 的转置为共轭转置，在MATLAB中使用单引号，x为输出信号，根据公式不需要转置，这里进行转置主要是为了进行矩阵运算。

2. pt_sig的计算是将每个阵元的对应时延与输出信号相乘得到修正后的信号，然后再相加得到频率响应（频率响应的定义为响应 $h(t)$ 和基函数 $f(t)$ 的卷积累加）：

如果还是不太明白可以跑出代码后对每个矩阵参数的行列数进行细致观察，从而加深理解。

结果：得到了s_sig，其为i*LEN大小的数据，i为theta的长度。矩阵一列的意义为：每个角度下的频率波数响应；矩阵一行的意义为：同一角度下不同采样点的频率波数响应。

3. 结果展示

分析：出现了两个波峰，一个对应60°，一个对应120°，暂且不清楚原因，后续可能会跟进修正，恳请知道原因的朋友赐教，以上。

原因解释：

由于代码基函数部分得到的信号只取了实数部分，后续与修正时延tao的共轭转置相乘时无法通过正交消除120°干扰，如下：

snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的水听器输出时间波形实数矩阵，Add white Gaussian noise to signal
x = rx.';          % 输出波形矩阵转置，便于二维数组FFT列运算，LEN x N，即接收到的信号

修改代码：

%% 频域常规波束形成
tic
xx = (as*s).';
for i = 1:length(theta)tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延pt_sig = exp(1j*2*pi*f0*tao)' * xx.';  %波束，1*LENs_sig(i,:) = pt_sig/N;    %size(s_sig) = i*LEN
end
toc
plot(theta,s_sig(:,1),'b.-');title("频域常规波束形成");

得到图像：

附修改后的MATLAB代码：

clear
close all
clc% 度与弧度常数
deg2rad = pi/180; % 度转换成弧度
rad2deg = 180/pi; % 弧度转换成度%% 水听器阵列参数设置
c = 1500;             % 水体声速，米/秒 m/s
N = 70;               % 阵元数
sqrtN=sqrt(N);        % 驾驶向量归一化常数
T=1;                 % 数据时长，秒sec
FS=24e2;              % 采样频率，赫兹Hz
LEN=T*FS;             % 数据长度，点数
t=(0:LEN-1)./FS;      % 数据时间轴，秒sec，1 x LEN 
fc=300;               % 阵列中心频率，赫兹Hz
lambda = c/fc;        % 波长，米m
lambda2 = lambda / 2; % 半波长，米m   d<(lambda/2)?
d = 0.27;             % 阵元间距，米m
ra=d*(0:N-1).';       % 阵元坐标，N x 1   %% 目标方位角设置
thetas_deg=90;                 % 目标方位角，度deg
thetas_rad=thetas_deg*deg2rad; % 目标方位角，弧度rad
theta = 0:0.1:180;                  % 任意信号入射角度%% 目标信号生成
f0=300;               % 声源频率，Hz
s=exp(1j*2*pi*f0.*t); % 目标复数信号，1 x LEN  "1j为复数，1为数字一，j为负数"  
as=exp(1j*2*pi*ra*f0*cos(thetas_rad)/c); % 目标驾驶向量，N x 1，驾驶向量=exp(-2*pi*f0*τ*1j*)
snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的水听器输出时间波形实数矩阵，Add white Gaussian noise to signal
x = rx.';                                 % 输出波形矩阵转置，便于二维数组FFT列运算，LEN x N，即接收到的信号%% 频域常规波束形成
tic
xx = (as*s).';
for i = 1:length(theta)tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延pt_sig = exp(1j*2*pi*f0*tao)' * xx.';  %波束，1*LENs_sig(i,:) = pt_sig/N;    %size(s_sig) = i*LEN
end
toc
figure(1);
subplot(2,1,1);
plot(theta,s_sig(:,1),'b.-');title("频域常规波束形成");  
subplot(2,1,2);
sig_db = 20*log10(abs(s_sig(:,1).^2));
plot(theta,sig_db,'r.-');title("频域常规波束形成——分贝图");

如果本篇文章对你有帮助，解决了你的疑惑或者问题，请不要吝啬你的点赞哦，谢谢！

clear
close all
clc% 度与弧度常数
deg2rad = pi/180; % 度转换成弧度
rad2deg = 180/pi; % 弧度转换成度%% 阵列参数设置
c = 1500;             % 水体声速，米/秒 m/s
N = 35;               % 阵元数
sqrtN=sqrt(N);        % 驾驶向量归一化常数
T=1;                 % 数据时长，秒sec
FS=24e3;              % 采样频率，赫兹Hz
LEN=T*FS;             % 数据长度，点数
t=(0:LEN-1)./FS;      % 数据时间轴，秒sec，1 x LEN 
fc=300;               % 阵列中心频率，赫兹Hz
lambda = c/fc;        % 波长，米m
lambda2 = lambda / 2; % 半波长，米m
d = 0.27;             % 阵元间距，米m
ra=d*(0:N-1).';       % 阵元坐标，N x 1%% 方位角设置
thetas_deg=90;                 % 目标方位角，度deg
thetas_rad=thetas_deg*deg2rad; % 目标方位角，弧度rad%% 目标信号生成
f0=3000;               % 声源频率，Hz
f1=f0;                % 起始处理频率，Hz
B=10;                 % 处理带宽，Hz
s=exp(1j*2*pi*f0.*t); % 目标复数信号，1 x LEN
sf=fft(s.',LEN);      % 目标信号FFT，LEN x 1
ks=2*pi*f0/c;         % 声源波数，1/m
% 指向目标方位的驾驶向量生成
as=exp(-1j*2*pi*ra*f0*cos(thetas_rad)/c); % 目标驾驶向量，N x 1
rf=sf*(as.');         % 指向目标方向的输出复数信号，LEN x N
% 指向目标方位的输出信号生成
snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的输出时间波形实数矩阵
x = rx.';                                 % 输出波形矩阵转置，便于二维数组FFT列运算，LEN x N
x_fft = fft(x,LEN);   % 对实信号进行FFT%% 基础数据
df = FS / LEN;        % 采样分辨率，即某个采样点对应的频率
figure(1);plot((0:LEN-1)*FS / LEN,abs(x_fft(:,1)));  % 测试图像，观察目标信号傅里叶变换：频率为f0 = 300
place = round(f0/df) + 1;    % 获取信号FFT频率对应位置
sig = x_fft(place,:);    % 获取频率为f0信号的FFT结果  1*N   频率有0频
theta_N = 0:0.1:180;  % 角度遍历
ntheta = length(theta_N);%% 常规频域窄带波束
for i = 1:nthetatao = cos(theta_N(i)*deg2rad)*ra/c;  % 在theta_N(i)角度下的时延  N*1A = exp(1j*2*pi*f0*tao);B = A' * sig.';    % 常规窄带波束  1*N*N*1pt_sig(i) = B;power_sig(i) =  B^2;  % 功率tt(i) = cosd(theta_N(i));
endpt_abs = abs(pt_sig);
pt2one = pt_abs ./ max(pt_abs); % 归一化power_abs = abs(power_sig);
power2one = power_abs ./ max(power_abs);  % 归一化figure(2)
plot(tt,20*log10(pt2one),'r');
% plot(theta_N,(pt2one),'r');
grid onfigure(3)
plot(tt,10*log10(power2one),'b');    % 功率的是10
% plot(theta_N,(power2one),'b');
grid on

运行结果：