宽带波束形成及MATLAB实现

宽带波束形成

宽带波形形成主要分两类，一类是基于频域的宽带波束形成，一类是基于时域的宽带波束形成。

频域宽带波束形成

对于有$M$个阵元的阵列，将各个阵元接收的原始数据进行分段，可以得到如下结果
$$\begin{array}{rl}& x^{\left(n\right)}=\left[\begin{array}{ccccc}{x}^{\left(n\right)}\left(1\right)& \cdots & x^{\left(n\right)}\left(l\right)& \cdots & x^{\left(n\right)}\left(L-1\right)\end{array}\right]\\ & =\left[\begin{array}{ccccc}{x}_1^{\left(n\right)}\left(0\right)& \cdots & x_1^{\left(n\right)}\left(l\right)& \cdots & x_1^{\left(n\right)}\left(L-1\right)\\ ⋮& \cdots & ⋮& \cdots & ⋮\\ x_m^{\left(n\right)}\left(0\right)& \cdots & x_m^{\left(n\right)}\left(l\right)& \cdots & x_m^{\left(n\right)}\left(L-1\right)\\ ⋮& \cdots & ⋮& \cdots & ⋮\\ x_M^{\left(n\right)}\left(0\right)& \cdots & x_M^{\left(n\right)}\left(l\right)& \cdots & x_M^{\left(n\right)}\left(L-1\right)\end{array}\right]\end{array}$$
对上述数据进行FFT，按照频域快拍进行波束形成，其波束形成方法可按照窄带波束形成的方法。

在将时域数据变换到频域的过程中，由于有限长度的数据变换到频域上面对应的频谱应该是无限宽的，而进行波束形成的时候只取了有限个子频带，相当于在频域上进行了滤波，同时由于滤波并不是理想滤波，会存在滤波器的建立时间问题，这样就会导致在不同数据块之间获得的波束输出序列与理想的波束输出序列存在畸变，即会出现不连续的输出结果。这也是频域波束形成存在的问题。此外，频域波束形成需要现将阵列接收的数据进行缓存，也就意味着存在着实时性的问题。

时域宽带波束形成

对于有$M$个阵元的阵列，第$m$号阵元接收的数据可以表示为
$$x_m\left(t\right)=s_0\left[t-{\tau }_m\left({\theta }_0\right)\right]+\sum _{k=1}^K{s}_k\left[t-{\tau }_m\left({\theta }_k\right)\right]+n_m\left(t\right),m=1,2,\cdots ,M$$
其中， { s k ( t ) } k = 0 K \left\{ {{s}_{k}}\left( t \right) \right\}_{k=0}^{K} {sk​(t)}k=0K​表示$K+1$个信号的波形，${\tau }_m\left({\theta }_k\right)$表示第$k$个信号到第$m$个阵元的延时，$n_m\left(t\right)$表示第$m$个阵元接收到的噪声。
如果波束主瓣方向${\theta }_s$刚好指向期望信号方向${\theta }_s$，那么各个阵元的预延时可以表示为
$$T_m=-\tau \left({\theta }_0\right)$$
其作用在于使期望方向信号到达各个阵元具有相同的相位。在实际中，由于该延时往往不是采样周期的整数倍，所以经常采用模拟预延时，在进行采样得到数字信号，因此第$m$个阵元的预延时接收数据样本可以表示为
$$x_m\left(i\right)=x_m\left(t-T_m\right){|}_{t=i{T}_s}$$
各个阵元通过横向的FIR滤波器，设滤波器的长度为$L$，采样周期为$T_s$，各个快拍输入数据可以表示为
$$x_{ml}\left(i\right)=x_m\left[i-\left(l-1\right)\right]=x_m\left[t-T_m-\left(l-1\right)T_s\right]{|}_{t=i{T}_s}m=1,2,\cdots ,M,l=1,2,\cdots ,L$$
对所有数据加权求和就可以得到波束输出的时间序列，对于第$m$号阵元第$l$个快拍的权值为$h_{ml}$，则波束形成输出的时间序列可以表示为
$$y\left(i\right)=\sum _{m=1}^M\sum _{l=1}^L{h}_{ml}{x}_m\left[i-\left(l-1\right)\right]=\sum _{m=1}^M\sum _{l=1}^L{h}_{ml}{x}_{ml}\left(i\right)$$
定义
$$\mathbf{X}\left(i\right)=\left[\begin{array}{cccc}{x}_{11}\left(i\right)& x_{12}\left(i\right)& \cdots & x_{1L}\left(i\right)\\ x_{21}\left(i\right)& x_{22}\left(i\right)& \cdots & x_{2L}\left(i\right)\\ ⋮& ⋮& \ddots & ⋮\\ x_{M1}\left(i\right)& x_{M2}\left(i\right)& \cdots & x_{ML}\left(i\right)\end{array}\right]$$
$$\mathbf{H}=\left[\begin{array}{cccc}{h}_{11}& h_{12}& \cdots & h_{1L}\\ h_{21}& h_{22}& \cdots & h_{2L}\\ ⋮& ⋮& \ddots & ⋮\\ h_{M1}& h_{M2}& \cdots & h_{ML}\end{array}\right]$$
令 $\mathbf{x}\left(i\right)=vec\left(\mathbf{X}\left(i\right)\right)$， $\mathbf{h}=vec\left(\mathbf{H}\right)$
则波束形成器的输出可以表示为
$$y\left(i\right)={\mathbf{h}}^T\mathbf{x}\left(i\right)$$

FIR滤波器设计

长度为$L$的滤波器的冲激响应为$\mathbf{h}={\left[h_1,h_2,\cdots ,h_L\right]}^T$，那么该滤波器的频率响应可以表示为
$$H\left(f\right)=\sum _{l=1}^L{h}_l{e}^{-j2\pi \left(l-1\right)f/f_s}={\mathbf{e}}^T\left(f\right)\mathbf{h}={\mathbf{h}}^T\mathbf{e}\left(f\right)$$
其中，$\mathbf{e}\left(f\right)={\left[1,e^{-j2\pi f/f_s},\cdots ,e^{-j2\pi \left(L-1\right)f/f_s}\right]}^T$，$f_s$为采样频率。
可以采用范数约束优化对其进行求解，若采用$l_2$范数进行优化，可以表示为
$$\begin{array}{rl}& \underset{\mathbf{h}}{\min }\sum _{k=1}^K\left[{\lambda }_k\Vert {\mathbf{e}}^T\left(f_k\right)\mathbf{h}-H_d\left(f_k\right)\Vert \right],f_k\in F_{ML}\\ & s.t.\Vert {\mathbf{e}}^T\left(f_k\right)\mathbf{h}\Vert \le {\xi }_0,f_k\in F_{SL}\end{array}$$
其中，$H_d\left(f_k\right)$表示期望的主瓣响应，${\xi }_0$表示设置的旁瓣值大小，${\zeta }_0$通过控制波束加权向量范数，提高稳健性，其值越小越稳健；${\lambda }_k$表示误差加权系数，一般取${\lambda }_j\text{=}1$。
时域与空域的对应关系如下：

时域	空域
时域滤波	空域滤波
FIR滤波器	波束形成器
频率响应$H(f)$	波束响应$p(\theta )$
滤波器系数$\mathbf{h}$	权系数$\mathbf{w}$
频率响应向量$\mathbf{e}(f)$	导向矢量$\mathbf{a}(\theta )$
滤波器设计	波束形成器设计
谱估计	空间谱估计

仿真结果

下面给出一个时域宽带波束形成的具体的仿真结果，仿真参数如下

仿真参数	参数值
阵元数目	20
期望信号方向	$10^{\circ }$
中心频率	$f_0=100Hz$
采样率	$f_s=5f_0$
滤波器长度	$L=65$

仿真结果如下：

从图中可以看出，基于FIR滤波器设计的宽带波束方向图，能够在主瓣方向形成最大增益，在其他方向形成抑制，也达到了波束形成器的功能，但是从图中也可以看出，不同频率对应的波束宽度也不一样。使用切比雪夫加权，可以得到等旁瓣的波束方向图，结果如下：

这样能更清晰的看出波束方向图的频率偏移性，解决这种问题可以采用凸优化的方法，这里就不再叙述，具体请参考宽带波束形成-----恒定波束宽度设计。三维方向图如下所示：

代码如下：

clear;
close all;
clc;
warning off;

%% 参数设置
M=20;     %阵元数
f0=100;   %中心频率
c=340;     %速度
lamda=c/f0;    %波长
d=lamda/2;     %阵元间距
theta0=[10];     %期望信号方向
fl=f0/2;             %信号的最低频率
fh=2*f0;            %信号的最高频率
fs=5*f0;             %采样率
N=512;              %信号点数
T=512/fs;           %信号周期
t=(0:N-1)/fs;       %时间轴
s=sin(2*pi*(fl+(fh-fl)/(2*T)*t).*t);      %产生宽带信号源

%% 时域FIR滤波器
L=65;                 %滤波器长度
D=(L-1)/2;          %预延迟
tao=(0:M-1).'*d*sind(theta0(1))/c;         %信号延时
f_delta=(fh-fl)/40;         %频率间隔
fd=(fl:f_delta:fh)./fs;             %归一化频率范围
w_win=chebwin(M,30);  %chebwin加权系数，-30dB

for m=1:M
    %% 期望滤波器的响应
    H(m,:)=w_win(m)/sum(w_win)*exp(-1i*2*pi*fd*(D+round(tao(m)*fs)-tao(m)*fs));     %预延迟与小数延迟后的总延迟求得的期望
    %频率响应
    MB=(fl:f_delta:fh)./fs;                                                                %归一化主瓣范围（工作频带范围）
    SB=[(0:f_delta:fl-8*f_delta)./fs (fh+8*f_delta)/fs:f_delta/fs:0.5];    %归一化旁瓣范围
    e_MB=exp(-1i*2*pi*(0:L-1)'*MB);                 %主瓣频率响应向量
    e_SB=exp(-1i*2*pi*(0:L-1)'*SB);                    %旁瓣频率响应向量
    
    %% 优化求解滤波器系数
    cvx_begin
    variable h(L)
    minimize(norm(e_MB.'*h-H(m,:).'));               %旁瓣频率响应约束的条件下
    subject to                                                    %最小化主瓣频率响应与期望响应的误差
    max(abs(e_SB.'*h))<=10^(-40/20);            %求出的滤波器响应
    cvx_end
    h_get(:,m)=h;
end

%% 波束方向图
theta_range=-90:90;           %扫描角度范围
T = -round(tao*fs + D)/fs;   %时间轴
p_total=[];
for i=1:length(MB)
    for j=1:length(theta_range)
        U=exp(-1i*2*pi*MB(i)*fs*((T+(0:M-1).'*d*sind(theta_range(j))/c)*ones(1,L)+ones(M,1)*(0:L-1)/fs));
        P(i,j)=vec(h_get.').'*vec(U);
    end
    figure(1);
    plot(theta_range,20*log10(abs(P(i,:))),'b','linewidth',1);
    hold on;
    ylim([-100 10]);
    xlabel('\theta(\circ)','fontsize',10);
    ylabel('Beampattern/dB','fontsize',10);
    p_total=[p_total;P(i,:)];
end
%% 三维方向图
figure(2);
[x,y]=meshgrid(theta_range,MB);
surf(x,y,20*log10(abs(p_total)));
xlabel('\theta(\circ)','fontsize',10)
ylabel('Normalized frequency','fontsize',10)
zlabel('Beampattern/dB','fontsize',10)
% colorbar;
zlim([-100 10]);

参考文献：
[1]鄢社锋, 马远良. 传感器阵列波束优化设计及应用[M]. 科学出版社, 2009.
[2]Yan S, Ma Y. A unified framework for designing FIR filters with arbitrary magnitude and phase response[J]. Digital Signal Processing, 2004, 14(6): 510-522.