【实验5】——波束形成技术

概念

目的：使多元阵构成的基阵经过适当处理得到在预定方向上的指向性。
意义：对于一个发射系统，具有指向性意味着发射能量可集中在某一方向，这样可以用较小的发射功率他侧更远距离的目标；对于一个接收系统，可以使系统定向接收，从而抑制其他方向的信号个干扰，准确测定目标方位。
方法：对各个阵元进行处理（加权、延时、求和等）形成空间指向性。

自然指向性与任意指向性

　我们设第0个阵元接收到的信号为cos(ωt)cos(\omega t)cos(ωt)，那么第1个阵元接收到的信号为cos(ω(t+τ))cos(\omega (t+\tau))cos(ω(t+τ))，第n个阵元接收到的信号为cos(ω(t+(n−1)τ))cos(\omega (t+(n-1)\tau))cos(ω(t+(n−1)τ))，其中τ=dsin(θ)c\tau=\frac{dsin(\theta)}{c}τ=cdsin(θ)，如果对所有阵元的信号进行一个求和，可以得到指向性因子D(θ)D(\theta)D(θ)：
　D(θ)=sin(Nπdλsin(θ))Nsin(πdλsin(θ))D(\theta)=\frac{sin(\frac{N\pi d}{\lambda } sin(\theta))}{Nsin(\frac{\pi d}{\lambda } sin(\theta))}D(θ)=Nsin(λπdsin(θ))sin(λNπdsin(θ))
　如果我们对每个阵元进行不同的时延，使其接收到的信号分别为cos(ωt)cos(\omega t)cos(ωt)，cos(ω(t+τ−ϕ))cos(\omega (t+\tau-\phi))cos(ω(t+τ−ϕ))，……，cos(ω(t+(n−1)(τ−ϕ)))cos(\omega (t+(n-1)(\tau-\phi)))cos(ω(t+(n−1)(τ−ϕ)))，其中ϕ=dsin(θ0)c\phi=\frac{dsin(\theta_0)}{c}ϕ=cdsin(θ0)，此时的指向性因子：
　D(θ)=sin(Nπdλ(sin(θ)−sin(θ0)))Nsin(πdλ(sin(θ)−sin(θ0)))D(\theta)=\frac{sin(\frac{N\pi d}{\lambda } (sin(\theta)-sin(\theta_0)))}{Nsin(\frac{\pi d}{\lambda } (sin(\theta)-sin(\theta_0)))}D(θ)=Nsin(λπd(sin(θ)−sin(θ0)))sin(λNπd(sin(θ)−sin(θ0)))

空间傅里叶变换角度去看波束形成

如果我们将信号写成复数的形式，那么相差可以用指数表示ejωnτe^{j\omega n\tau}ejωnτ，记第n个阵元的幅值为s((n−1)d)s((n-1)d)s((n−1)d)，指向性因子可以写成下面的形式：
D(θ)=∑n=0N−1s(nd)ejωnτ=∑n=0N−1s(nd)ejωsin(θ)cnd=∑n=0N−1s(nd)ejksin(θ)nd=∑n=0N−1s(nd)ejkxndD(\theta)=\sum_{n=0}^{N-1}s(nd)e^{j\omega n\tau}=\sum_{n=0}^{N-1}s(nd)e^{j\frac{\omega sin(\theta)}{c}nd}=\sum_{n=0}^{N-1}s(nd)e^{j ksin(\theta)nd}=\sum_{n=0}^{N-1}s(nd)e^{j k_xnd}D(θ)=n=0∑N−1s(nd)ejωnτ=n=0∑N−1s(nd)ejcωsin(θ)nd=n=0∑N−1s(nd)ejksin(θ)nd=n=0∑N−1s(nd)ejkxnd
　从上面的推导可以看出，求指向性就是一个空间域的离散序列的傅里叶变换（对水平波数作变换）。
　根据傅里叶变换到的知识，离散序列的谱是周期的：
　
　我们注意到，kxk_xkx的取值是有范围的，为[−kxmax,kxmax][-k_{xmax},k_{xmax}][−kxmax,kxmax]，所以实际只能看到部分的谱，kxmaxk_{xmax}kxmax越大代表能看到的谱就越多，为了不看到第二个峰（即副极大），要使
　kxmax<2πdk_{xmax}<\frac{2\pi}{d}kxmax<d2π
　而在任意指向性时，kx=k(sin(θ)−sin(θ0))k_{x}=k(sin(\theta)-sin(\theta_0))kx=k(sin(θ)−sin(θ0))，所以kxmax=2kk_{xmax}=2kkxmax=2k，故
　2k<2πd2k<\frac{2\pi}{d}2k<d2πd<λ2d<\frac{\lambda}{2}d<2λ
　上式即为空间采样定理，我认为和时间的采样定理还是有一定区别的，时间的采样定理的描述是截止频率<采样率2截止频率<\frac{采样率}{2}截止频率<2采样率，为的是防止频谱混叠；而空间采样定理是为了不看到第二个峰。

主瓣宽度讨论

基于上面两种方法的分析，我们可以从两个角度进行讨论。
　首先从D(θ)D(\theta)D(θ)的表达式入手，通过近似计算，得到主瓣的宽度表达式为：
　θ−θ0=0.44λNdcos(θ0)\theta-\theta_0=\frac{0.44\lambda}{Ndcos(\theta_0)}θ−θ0=Ndcos(θ0)0.44λ

阵长NdNdNd越长，主瓣宽度越小；
声波入射角越偏向端射方向，主瓣宽度越大；
频率越高，主瓣宽度越小；

从傅里叶变换角度考虑：

阵元越长，即空域越宽，波数域越窄，即主瓣宽度约小；
频率改变，波数谱的形状不会变，影响的只是看到的宽度；因为频率越高，波数越大，看到的宽度越大，把看到的谱映射到角度域时，内容更多了，但定义域还是[−90,90][-90,90][−90,90]，所以主瓣的宽度变窄了；

注：以上所有从傅里叶变换角度的分析都是个人理解，可能会存在一定问题。

代码实现

总体分两步：

空间积分，即各阵元的信号求和
时间积分，取信号的能量（信号处理不能通过信号的瞬时值来进行判决，能量正比于幅度的平方）

%阵元的自然指向性
clc;clear all;
fs=1000;%采样率
c=1500;%声速
f0=20;%信号频率 波长75
d=20;%阵元间隔
ts=1;%采样的时间长度
t=1/fs:1/fs:ts;%时间序列
N=40;%阵元个数
target=60;%目标角度
theta=-90:90;
for count=1:length(theta)for n=1:Ntao(n)=(n-1)*d*sind(theta(count))/c;%第n号阵元的时延fai(n)=(n-1)*d*sind(target)/c;%给第n号阵元添加的相位s(n,:)=cos(2*pi*f0*(t+tao(n)));%第n号阵元接收到的信号s1(n,:)=cos(2*pi*f0*(t+tao(n)-fai(n)));%形成具有固定指向性的信号endSt=sum(s);%对所有阵元的信号进行空间求和St1=sum(s1);E(count)=sum(St.^2);%对时间积分，得到能量E1(count)=sum(St1.^2);
end
E=E/max(E);%能量归一化
E1=E1/max(E1);
figure;
subplot(121);
plot(theta,E);xlabel('角度');title('指向性图(直角坐标系)');
subplot(122);
polar(theta*pi/180,E);title('指向性图(极坐标系)');
figure;
subplot(121);
plot(theta,E1);xlabel('角度');title('指向性图(直角坐标系)');
subplot(122);
polar(theta*pi/180,E1);title('指向性图(极坐标系)');

结果：