阵列信号处理笔记
文章目录
- 阵列信号处理笔记
- 阵列调向
- 栅瓣
- 调向对方向图的影响
- 调向对HPBW的影响
- 工程相关
- MATLAB代码
- 其他乱七八糟的
阵列调向
阵列调向是阵列信号处理中一个非常基本的问题，通常阵列调向需要做的就是将主瓣对准某一特定位置。这一骚操作可以通过物理手段实现，例如旋转阵列朝向，或者通过电子方式实现。一个比较经典的例子是基站天线的下倾角调节，如图所示。
这种天线通过一个可调节的机械结构改变天线对地的夹角，而为了调节这个机械结构需要相关的操作人员上站调试。现代基站天线基本都具有电子下倾角调节功能，当需要对基站覆盖范围做出某些调整（例如相邻小区干扰严重）时，操作人员在机械倾角的基础上通过合理设置阵列参数，就能达到等效的倾角调节。天线下倾角管理是一个优化问题，而它的基础则依赖于阵列调向理论。
阵列调向理论非常简单，根据前面的章节，任意一个阵列的频域-波数域响应都能用下式表示
γ
(
w
,
k
)
=
W
H
(
w
)
v
(
k
)
(1)
\gamma(w,\mathbf{k})=\mathbf{W}^{\mathrm{H}}(w)\mathbf{v}(\mathbf{k})\tag{1}
γ(w,k)=WH(w)v(k)(1)
阵列调向做的事情在波数域（k域）很简单，形式上就是
γ
(
w
,
k
−
k
0
)
\gamma(w,\mathbf{k-k_0})
γ(w,k−k0)，从数学上则是阵列流形左乘一个调向矩阵
I
s
(
k
0
)
\mathbf{I_s}(\mathbf{k_0})
Is(k0)
γ
(
w
,
k
−
k
0
)
=
W
H
(
w
)
v
(
k
−
k
0
)
=
W
H
(
w
)
I
s
(
k
0
)
v
(
k
)
(2)
\gamma(w,\mathbf{k-k_0})=\mathbf{W}^{\mathrm{H}}(w)\mathbf{v}(\mathbf{k-k_0})=\mathbf{W}^{\mathrm{H}}(w)\mathbf{I_s}(\mathbf{k_0})\mathbf{v}(\mathbf{k})\tag{2}
γ(w,k−k0)=WH(w)v(k−k0)=WH(w)Is(k0)v(k)(2)
其中
I
s
(
k
0
)
=
(
e
j
k
0
T
p
0
0
⋯
0
0
e
j
k
0
T
p
1
⋯
0
0
⋯
⋱
0
0
⋯
0
e
j
k
0
T
p
N
−
1
)
(3)
\mathbf{I_s}(\mathbf{k_0})=\left(\begin{array}{c} e^{j \mathbf{k_0}^\mathrm{T}p_0} & 0 & \cdots &0 \\ 0 &e^{j \mathbf{k_0}^\mathrm{T}p_1} & \cdots &0 \\0&\cdots &\ddots &0 \\ 0 & \cdots &0 &e^{j \mathbf{k_0}^\mathrm{T}p_{N-1}} \\ \end{array}\right)\tag{3}
Is(k0)=
ejk0Tp00000ejk0Tp1⋯⋯⋯⋯⋱0000ejk0TpN−1
(3)
有兴趣可以乘一乘验证一下😀对于u域，
ψ
\psi
ψ域也是如法炮制，但是对
θ
\theta
θ域稍微有点区别。实际上对于ULA来说，调向后的频域-波数域响应就是
γ
(
k
:
k
0
)
=
1
N
v
(
k
0
)
H
v
(
k
)
(4)
\gamma(\mathbf{k}:\mathbf{k_0})=\frac{1}{N}\mathbf{v}(\mathbf{k_0})^{\mathrm{H}}\mathbf{v}(\mathbf{k})\tag{4}
γ(k:k0)=N1v(k0)Hv(k)(4)
证明很简单，令
W
=
[
1
/
N
,
1
/
N
,
⋯
,
1
/
N
]
\displaystyle \mathbf{W}=[1/N,1/N,\cdots, 1/N]
W=[1/N,1/N,⋯,1/N]即可。
栅瓣
栅瓣是波束方向图中间隔出现的，与主瓣高度相同的波瓣。栅瓣的来源可以通过空域采样理论解释，即一个域的离散必定对应于另一个域的周期延拓。在时域中对一个信号 f ( t ) f(t) f(t)进行采样，那么频域必定得到的是以采样角频率为周期重复的 F ( w ) F(w) F(w)。与时域采样类似，空域采样也会存在混叠。关于空域混叠问题在上面一节已经给出了一个描述，这里将给出另一个等价的说法。
阵列调向的方向图可以视为一种移动（在k/u/
ψ
\psi
ψ域），以u域为例，其中
u
0
u_0
u0是移动距离，那么
γ
(
u
:
u
0
)
=
1
N
sin
[
π
N
λ
d
(
u
−
u
0
)
]
sin
[
π
λ
d
(
u
−
u
0
)
]
(5)
\gamma(u:u_0)=\frac{1}{N} \frac{\sin\left[ \frac{\pi}{N\lambda}d(u-u_0) \right]}{\sin \left[ \frac{\pi}{\lambda}d\left( u-u_0 \right) \right]}\tag{5}
γ(u:u0)=N1sin[λπd(u−u0)]sin[Nλπd(u−u0)](5)
其中可视区域为
−
1
⩽
u
⩽
1
-1 \leqslant u \leqslant1
−1⩽u⩽1（即
0
⩽
θ
⩽
18
0
∘
0 \leqslant \theta \leqslant 180 ^{\circ}
0⩽θ⩽180∘），当
u
0
∈
[
−
1
,
1
]
u_0 \in [-1, 1]
u0∈[−1,1]时，可以绘制相应的波束方向图。下面四幅图给出了阵元间距d分别为
λ
\lambda
λ、
0.75
λ
0.75\lambda
0.75λ、
0.5
λ
0.5\lambda
0.5λ、
0.25
λ
0.25\lambda
0.25λ时的动态图。
这是 d = λ d=\lambda d=λ的情况
这是 d = 0.75 λ d=0.75\lambda d=0.75λ的情况
这是 d = 0.5 λ d=0.5\lambda d=0.5λ的情况
这是 d = 0.25 λ d=0.25\lambda d=0.25λ的情况
不难发现，当 d > 0.5 λ d>0.5 \lambda d>0.5λ时，随着 u 0 u_0 u0在 [ − 1 , 1 ] [-1, 1] [−1,1]移动，栅瓣的峰值将落入可视区域内，这种情况被认为是一种峰值响应模糊，也就是说当给我们一张波束方向图（没有任何先验知识），我们无从分辨哪个是主瓣，继而做空间谱估计时也无从分辨哪个峰值对应了真正的波达方向。当 d = 0.5 λ d=0.5 \lambda d=0.5λ时恰好为临界情况，而 d < 0.5 λ d<0.5 \lambda d<0.5λ， ∀ u 0 ∈ [ − 1 , 1 ] \forall u_0 \in [-1,1] ∀u0∈[−1,1]均无栅瓣峰值落入可视区域，并且随着d的减小，栅瓣之间间距也越来越大。
所以这里我们从u域出发得到了相同的结论，即对于均匀线阵，满足无空域混叠的条件是
d
⩽
λ
2
(6)
d \leqslant \frac{\lambda}{2}\tag{6}
d⩽2λ(6)
现在考虑调向范围受限的情形，即
∣
u
0
∣
⩽
Δ
u
m
a
x
，
Δ
u
m
a
x
∈
[
0
,
1
]
\left| u_0 \right| \leqslant \Delta u_{max}，\Delta u_{max} \in [0, 1]
∣u0∣⩽Δumax，Δumax∈[0,1]。为了使得不存在峰值响应模糊，即第一级栅瓣不能落入可视区域，必须满足
d
⩽
1
1
+
Δ
u
m
a
x
λ
(7)
d \leqslant \frac{1}{1+\Delta u_{max}}\lambda\tag{7}
d⩽1+Δumax1λ(7)
简单证明一下。
考虑左侧第一级栅瓣，当主瓣调向达到右侧最大时，第一级栅瓣位于左侧可视区域边界之外（说人话就是可视区域里不能有栅瓣），也就是
π λ d ( − 1 − Δ u m a x ) ⩾ − π \frac{\pi}{\lambda}d(-1-\Delta u_{max})\geqslant-\pi λπd(−1−Δumax)⩾−π
即
d ⩽ 1 1 + Δ u m a x d \leqslant \frac{1}{1+ \Delta u_{max}} d⩽1+Δumax1
根据对称性，右侧第一级栅瓣同理。值得注意的是，当 Δ u m a x = 1 \Delta u_{max}=1 Δumax=1，即 u 0 ∈ [ − 1 , 1 ] u_0 \in [-1, 1] u0∈[−1,1]时，上式简化为
d ⩽ λ 2 d \leqslant \frac{\lambda}{2} d⩽2λ
这与前面MATLAB仿真结果是一致滴。
调向对方向图的影响
阵列调向的方向图在k、u、 ψ \psi ψ域可以视为一种移动（平移）。但是在使用的时候，k、u、 ψ \psi ψ域表述似乎实际意义不怎么明确，因为通常我们没法指着一个地方说握草这地方就是 ψ = π 3 \displaystyle \psi=\frac{\pi}{3} ψ=3π…但是我们能很容易感受 θ = 4 5 ∘ \displaystyle \theta=45 ^{\circ} θ=45∘所对应的大致朝向，但是很不幸的是，在 θ \theta θ域中调向对方向图的影响并非简单的搬移。这种影响要从u域和 θ \theta θ域的关系说起。
我们熟知 u = cos θ \displaystyle u = \cos \theta u=cosθ，需要注意的是这个 θ \theta θ指的是阵列扫描方向与z轴正方向的夹角（见上一篇文章开头的图），或者画在平面图上是这样的。
在很多文献（包括论文、书籍）写的跟我并不一样，他们可能取的是扫描方向与阵列法线的夹角。如果存在一些符号上的差异（例如我是 cos θ \cos \theta cosθ别人是 sin θ \sin \theta sinθ）请仔细甄别。其实原书《Optimum Array Processing》早在前面就已经对这一点做了详细区分，但是很多时候看起来实在有些累。
在u域的阵列调向，通常只需要做u域的平移，即
u
−
u
0
u-u_0
u−u0，这样就能将主瓣对准
u
0
u_0
u0所在的方向。有人会说，那我在
θ
\theta
θ域调向，不是只要做
θ
−
θ
0
\theta-\theta_0
θ−θ0就可以了么？？但是试过的都知道，理想是丰满滴，现实是骨感滴，因为cos并不满足移不变性（shift-invariance property），这件事情天然的就不成立。那为什么在其他域是成立的捏？？（因为是直接由Fouier transform的shift property推导得到）。正确的做法是回到u域，做
cos
θ
−
cos
θ
0
\cos \theta - \cos \theta _0
cosθ−cosθ0，这样一来，对于ULA来说
γ
(
θ
:
θ
0
)
=
1
N
sin
[
N
π
λ
d
(
cos
θ
−
cos
θ
0
)
]
sin
[
π
λ
d
(
cos
θ
−
cos
θ
0
)
]
(8)
\gamma(\theta:\theta_0)=\frac{1}{N} \frac{ \sin \left[ \frac{N\pi}{\lambda}d \left(\cos \theta - \cos \theta_0 \right) \right]}{\sin \left[ \frac{\pi}{\lambda}d \left(\cos \theta - \cos \theta_0 \right) \right]}\tag{8}
γ(θ:θ0)=N1sin[λπd(cosθ−cosθ0)]sin[λNπd(cosθ−cosθ0)](8)
对N=10的标准ULA，我们画出
θ
=
0
∘
\theta =0^{\circ}
θ=0∘、
θ
=
6
0
∘
\theta =60^{\circ}
θ=60∘、
θ
=
13
5
∘
\theta =135^{\circ}
θ=135∘的方向图（这里没有置于极坐标中是为了方便对比）
很显然这再一次印证了 θ \theta θ域的调向并不是简单的平移，因为主瓣、旁瓣形态都发生了一些改变。因此有必要重新研究一下调向过程中HPBW的变化。
如图所示，在u域中，HPBW可以简单表示为
H
P
B
W
=
u
r
−
u
l
=
(
u
0
+
0.45
λ
N
d
)
−
(
u
0
−
0.45
λ
N
d
)
(9)
\mathrm{HPBW}=u_r-u_l=\left( u_0+0.45\frac{\lambda}{Nd} \right)-\left( u_0-0.45\frac{\lambda}{Nd} \right)\tag{9}
HPBW=ur−ul=(u0+0.45Ndλ)−(u0−0.45Ndλ)(9)
利用u域与
θ
\theta
θ域的关系可得（注意arccos在[-1,1]是单调减函数）
θ
H
=
θ
r
−
θ
l
=
arccos
(
u
0
−
0.45
λ
N
d
)
−
arccos
(
u
0
+
0.45
λ
N
d
)
(10)
\theta_{H}=\theta_r-\theta_l=\arccos\left( u_0-0.45\frac{\lambda}{Nd} \right)-\arccos \left( u_0+0.45\frac{\lambda}{Nd} \right)\tag{10}
θH=θr−θl=arccos(u0−0.45Ndλ)−arccos(u0+0.45Ndλ)(10)
这个式子看似挺对的哈，其实不难发现可能存在一个
u
o
′
u_o'
uo′使得
u
0
′
+
0.45
λ
N
d
>
1
(11)
u_0'+0.45 \frac{\lambda}{Nd} > 1\tag{11}
u0′+0.45Ndλ>1(11)
这时候实数域上的反余弦已经没有定义了，很显然这个时候已经找不到
θ
l
\theta_l
θl了，这个极限角度就称为扫描限（Scan Limit）。那么这是什么情况呢…我一开始看教材，这段也挺迷的，然后MATLAB做了一下仿真豁然开朗。
上面这张动图做的是波束从 9 0 ∘ 90 ^{\circ} 90∘调向至 0 ∘ 0^{\circ} 0∘的情形，可以看到 θ H \theta_H θH在逐渐变大，在大概20°左右的时候就超出扫描限了，直观来看就是左半功率点没了，因为它与左侧过来的一个栅瓣融合在一起。当 θ 0 = 0 \theta_0=0 θ0=0时（也就是主瓣沿着阵列方向走）称为端射（End Fire），如下图（偷来的）左侧所示。
由（10）可知， θ H \theta_H θH是与波阵列孔径（就是 N d / λ Nd/\lambda Nd/λ）之比有关的，那么其实可以用下图描述其关系（注意这个图的横坐标和纵坐标都是对数的哈）。
那么这张图说明了一个什么意思呢？我认为大致有2个。
首先可以肯定的是，当调向越接近正向侧（Broadside，也就是 θ 0 = 9 0 ∘ \theta_0=90^\circ θ0=90∘）孔径越大HPBW越小，这意味着大孔径的阵列具有更加优异的角度分辨性能。
其次，当调向越接近端射，HPBW增大越快，也就是说角度分辨率恶化速度越快。这意味着在工程上，阵列调向角度不宜过大，过大的调向角度会导致主瓣宽度增加（另外还有自耦问题和栅瓣迫近的问题）。
调向对HPBW的影响
上面通过一张图定性的阐述一些调向对阵列性能的影响，这里还是要深入探讨一下。当调向角度接近正侧向（Broadside，
θ
=
π
2
\displaystyle \theta = \frac{\pi}{2}
θ=2π）时具有最窄的HPBW。对（10）使用拉格朗日中值定理可得
θ
H
=
0.891
λ
N
d
1
1
−
ξ
2
,
ξ
∈
[
cos
θ
0
−
0.45
λ
N
d
,
cos
θ
0
+
0.45
λ
N
d
]
(12)
\theta_H=0.891\frac{\lambda}{Nd}\frac{1}{\sqrt{1-\xi^2}},\xi\in[\cos \theta_0-0.45\frac{\lambda}{Nd},\cos \theta_0+0.45\frac{\lambda}{Nd}]\tag{12}
θH=0.891Ndλ1−ξ21,ξ∈[cosθ0−0.45Ndλ,cosθ0+0.45Ndλ](12)
显然当
N
d
/
λ
→
0
Nd/\lambda\to 0
Nd/λ→0时，
H
P
B
W
→
0
HPBW\to 0
HPBW→0，不妨取
ξ
\xi
ξ为区间中点，那么
θ
H
≈
0.891
λ
N
d
csc
θ
0
(13)
\theta_H\approx 0.891 \frac{\lambda}{Nd}\csc \theta_0\tag{13}
θH≈0.891Ndλcscθ0(13)
《阵列信号处理》原书这里写的是 csc θ T ‾ \csc\overline{\theta_T} cscθT，应该是一个笔误，如果用余角应该写 sec θ T ‾ \sec \overline{\theta_T} secθT。
那么（13）说滴究竟是什么意思捏？？我觉得需要改写一下这个式子，这样直观意义更加明显
θ
H
≈
0.891
λ
N
(
d
sin
θ
0
)
(14)
\theta _H\approx 0.891 \frac{\lambda}{N\left( d \sin \theta_0 \right)}\tag{14}
θH≈0.891N(dsinθ0)λ(14)
在上一篇文章曾贴了一张图
根据幼儿学学过的三角函数知识不难发现，此时阵元的等效的间距其实变成了 d sin θ 0 d \sin \theta_0 dsinθ0，如果带入上一篇文章讨论过的 θ \theta θ域的HPBW公式不难得到（14）的结论。也就是说，向端射侧调向相当于降低了阵列孔径，这也能解释为什么越接近端射主瓣宽度越宽。这种现象早在上一篇文章介绍空域采样理论解释的时候就已经能够预见。
当
θ
0
=
0
\theta_0=0
θ0=0时，也就是端射阵列时，很显然
θ
H
=
2
arccos
(
1
−
0.45
λ
N
d
)
(15)
\theta_H=2 \arccos \left( 1-0.45 \frac{\lambda}{Nd} \right)\tag{15}
θH=2arccos(1−0.45Ndλ)(15)
改写一下上式，并利用三角升幂与等价无穷小公式，可以得到
θ
H
≈
2
0.891
λ
N
d
(16)
\theta_H \approx 2 \sqrt{0.891 \frac{\lambda}{Nd}}\tag{16}
θH≈20.891Ndλ(16)
简单推导一下，俩改变同时取余弦并移项得到
1 − cos ( θ H 2 ) = 0.45 λ N d 1-\cos\left( \frac{\theta_H}{2} \right)=0.45 \frac{\lambda}{Nd} 1−cos(2θH)=0.45Ndλ
左边用一下升幂可以得到
sin ( θ H 4 ) = 0.45 λ 2 N d \sin \left( \frac{\theta_H}{4} \right)=\sqrt{\frac{0.45\lambda}{2Nd}} sin(4θH)=2Nd0.45λ
当 N d / λ → 0 Nd/\lambda\to 0 Nd/λ→0时， H P B W → 0 HPBW\to 0 HPBW→0，对左边用一下等价无穷小 sin ( θ H 4 ) ∼ θ H 4 \displaystyle \sin\left( \frac{\theta _H}{4} \right)\sim \frac{\theta _H}{4} sin(4θH)∼4θH，简化之可得
θ H ≈ 2 0.891 λ N d \theta_H \approx 2 \sqrt{0.891 \frac{\lambda}{Nd}} θH≈20.891Ndλ
工程相关
根据幼儿园学的信号与系统的相关经验，一个窗的长度越长主瓣宽度就越窄，这种奇妙的时频对偶特性同样体现在空域中，正如上图所示，阵列的孔径波越大，可以预见的是主瓣宽度也会越窄，HPBW也更小，相应的角度分辨率也会越大。然而大孔径有赖于阵列间距d和阵元数量N，根据前面提到的空域采样定理，增加d会导致不满足无混叠条件，增加N可能带来潜在成本增加（特别是硬件的开销和算力需求），并且过小的d又会使得天线尺寸受限而影响性能，这是一种compromise的艺术。
下图是TI公司最新的毫米波雷达SOC AWR2944评估版板载的阵列，其中左边是接收天线，右边是发射天线。该阵列的设计使用了一种叫做虚拟阵列的技术，通过特殊的非均匀阵元排列达到提升角度分辨率的效果。
这种骚操作可以在保证单个天线性能的前提下更好的平衡设计复杂度以及成本。
此外对于特定的应用，不需要如此之高的扫描范围，那么对阵列间距可以适当放宽，根据（7）的结论可以得到
d
⩽
1
1
+
∣
sin
θ
m
a
x
‾
∣
d \leqslant \frac{1}{1+\left| \sin \overline{\theta_{max}} \right|}
d⩽1+
sinθmax
1
其中这个
θ
m
a
x
‾
\overline{\theta_{max}}
θmax是最大调向角度的余角。
MATLAB代码
阵列调向u域平移的动图
clc;
clear;
close all;
lambda = 1;
d = lambda / 2;
for u0 = [-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]
hold off
[B, u] = get_B_u(10, d, lambda, u0);
plot(u, 20 * log10(abs(B)),'linewidth', 1.3);
Graph_Set(gca,gcf, [-4 4],[-30 20],'u','','N=10',[900,400]);
hold on;
fill([-1,-1,1,1],[20,-30,-30,20],'r','FaceAlpha',0.1);
text(-0.65,18,'Visible Region','FontSize',14);
text(-0.35,15,['u = ' num2str(u0, '%.2f')],'FontSize',16);
text(-0.35,12,['\theta = ' num2str(acos(u0) / pi * 180, '%.2f') '^{\circ}'],'FontSize',16);
set(gcf,'position',[500,500,900,400]);
drawnow;
pause(1);
end
function [B, u] = get_B_u(N, d, lambda, u0)
u = -4 : 0.001 : 4;
B = 1 / N * sin(N * pi / lambda * d * (u-u0)) ./ sin(pi / lambda * d * (u-u0));
end
function Graph_Set(ax, fig, XLim, YLim, Xlabel, Ylabel, Title, size)
ax.FontSize = 11;
ax.Title.String = Title;
ax.XLabel.String = Xlabel;
ax.YLabel.String = Ylabel;
ax.XLim = XLim;
ax.YLim = YLim;
ax.XGrid = 'on';
ax.YGrid = 'on';
ax.XMinorTick = 'on';
ax.YMinorTick = 'on';
set(fig,'position',[500,300, size]);
end
θ \theta θ域调向， θ = 0 ∘ \theta =0^{\circ} θ=0∘、 θ = 6 0 ∘ \theta =60^{\circ} θ=60∘、 θ = 13 5 ∘ \theta =135^{\circ} θ=135∘的方向图
clc;
clear;
close all;
lambda = 1;
[B, u] = get_B_theta(10, lambda * 0.5, lambda, 0);
B_db = 20 * log10(abs(B));
subplot(3,1,1);
hold on
plot(u, B_db,'linewidth', 1.3);
fill([-3,0,3],[5,1,5],'r');
hold off
Graph_Set(gca,gcf,[-180 180],[-30 20],'\theta','','\theta_0=0^{\circ}',[900,600]);
[B, u] = get_B_theta(10, lambda * 0.5, lambda, 60);
B_db = 20 * log10(abs(B));
subplot(3,1,2);
hold on
plot(u, B_db,'linewidth', 1.3);
fill([57,60,63],[5,1,5],'r');
hold off
Graph_Set(gca,gcf,[-180 180],[-30 20],'\theta','','\theta_0=60^{\circ}',[900,600]);
[B, u] = get_B_theta(10, lambda * 0.5, lambda, 135);
B_db = 20 * log10(abs(B));
subplot(3,1,3);
hold on
plot(u, B_db,'linewidth', 1.3);
fill([132,135,138],[5,1,5],'r');
hold off
Graph_Set(gca,gcf,[-180 180],[-30 20],'\theta','','\theta_0=135^{\circ}',[900,600]);
function [B, theta] = get_B_theta(N, d, lambda, theta0)
theta = -180 : 0.01 : 180;
B = 1 / N * sin(N * pi / lambda * d * (cos(theta/180*pi)-cos(theta0/180*pi))) ./ sin(pi / lambda * d * (cos(theta/180*pi)-cos(theta0/180*pi)));
end
function Graph_Set(ax, fig, XLim, YLim, Xlabel, Ylabel, Title, size)
ax.FontSize = 11;
ax.Title.String = Title;
ax.XLabel.String = Xlabel;
ax.YLabel.String = Ylabel;
ax.XLim = XLim;
ax.YLim = YLim;
ax.XGrid = 'on';
ax.YGrid = 'on';
ax.XMinorTick = 'on';
ax.YMinorTick = 'on';
set(fig,'position',[500,300, size]);
end
阵列调向 θ \theta θ域变化的动图
clc;
clear;
close all;
lambda = 1;
theta0 = 20;
N = 10;
d = lambda * 0.5;
for theta0 = 90 : -5 : 0
hold off
[B, theta] = get_B_theta(N, d, lambda, theta0);
B_db = 20 * log10(abs(B));
plot(theta, B_db,'linewidth', 1.3);
Graph_Set(gca,gcf,[-180 180],[-30 20],'\theta','','',[900,400]);
hold on;
theta_hpbw_r = acos(cos(theta0 / 180 * pi) - 0.45 * lambda /(N * d)) / pi * 180;
theta_hpbw_l = acos(cos(theta0 / 180 * pi) + 0.45 * lambda /(N * d)) / pi * 180;
text(-70,16,['\theta = ' num2str(theta0, '%.2f') '^{\circ}'],'FontSize',20);
if(real(theta_hpbw_l) ~= 0)
fill([theta_hpbw_l,theta_hpbw_l,theta_hpbw_r,theta_hpbw_r],[20,-30,-30,20],'r','FaceAlpha',0.1);
text(-70,11,['\theta_H = ' num2str(theta_hpbw_r - theta_hpbw_l, '%.2f') '^{\circ}'],'FontSize',20)
else
text(-120,11,'Exceed Scan Limit','FontSize',20);
end
drawnow;
pause(0.5);
end
function [B, theta] = get_B_theta(N, d, lambda, theta0)
theta = -180 : 0.01 : 180;
B = 1 / N * sin(N * pi / lambda * d * (cos(theta/180*pi)-cos(theta0/180*pi))) ./ sin(pi / lambda * d * (cos(theta/180*pi)-cos(theta0/180*pi)));
end
function Graph_Set(ax, fig, XLim, YLim, Xlabel, Ylabel, Title, size)
ax.FontSize = 11;
ax.Title.String = Title;
ax.XLabel.String = Xlabel;
ax.YLabel.String = Ylabel;
ax.XLim = XLim;
ax.YLim = YLim;
ax.XGrid = 'on';
ax.YGrid = 'on';
ax.XMinorTick = 'on';
ax.YMinorTick = 'on';
set(fig,'position',[500,200,size]);
end
HPBW与阵列孔径、调向角度的关系图
这是原书自带的代码，加了点东西美化一下（面子工程）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.24
% HPBW versus steering angle
% Lillian Xu
% Modified by Xin Zhang
% Last updated by L. Xu 11/30/00, K. Bell 7/22/01, 10/4/01, 10/17/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
Nd_set=[1000:-5:700 700:-0.5:200 200:-0.05:20 15:-0.005:1];
theta_deg=[2.5 5 10 20 30 45 90];
%theta_deg=[2.5 5 10 20 30 45 60 90];
theta_set=theta_deg*pi/180;
figure
for num=1:size(theta_set,2)
theta=theta_set(num);
bw=0;
for num1=1:size(Nd_set,2)
Nd=Nd_set(num1);
r=cos(theta)+0.443/Nd;
if abs(r)<=1
bw(num1)=acos(r)-acos(cos(theta)-0.443/Nd);
end
end
loglog(Nd_set(1:size(bw,2)),abs(bw*180/pi),'linewidth', 1.3);
hold on
rightpoint=abs(bw(1)*180/pi);
text(1040,rightpoint,[num2str(theta_deg(num)),'^\circ'],'fontsize',10); % used to be 7
end
axis([1 1000 0.04 100])
xlabel('{\itNd}/\lambda','Fontsize',16)
ylabel('3-dB beamwidth in degrees','Fontsize',16)
scanlimit=-acos(1-2*0.443./Nd_set);
loglog(Nd_set,abs(scanlimit*180/pi),'--','linewidth', 1.3)
Endfire = 2*acos(1-0.433./Nd_set);
loglog(Nd_set,abs(Endfire*180/pi),'linewidth', 1.3)
rightpoint=abs(Endfire(1)*180/pi);
text(1040,rightpoint,'\theta=0^\circ','fontsize',10); % used to be 7
%title('HPBW versus steering angle: standard linear array with uniform weighting')
text(12,35,'Endfire','Fontsize',12)
text(12,14.4,'Scan limit','Fontsize',12)
text(12,1.8,'Broadside','Fontsize',12)
plot([1 1000],[0.1 0.1],':')
plot([1 1000],[1 1],':')
plot([1 1000],[10 10],':')
plot([1 1000],[100 100],':')
plot([10 10],[0.04 200],':')
plot([100 100],[0.04 200],':')
axis([1 1000 0.04 200])
set(gcf,'position',[500,200,1200,600]);
gca.XGrid = 'on';
gca.YGrid = 'on';
gca.XMinorTick = 'on';
gca.YMinorTick = 'on';
grid on
其他乱七八糟的
这篇文章里用了几个动图，感觉效果还挺不错的。其实本来MATLAB就可以做gif，不过出来的文件实在有些大啊，这里还是采用了保存图片+第三方软件生成gif。
MATLAB保存gif的语句是
print(gcf,'-dpng','example.png');
这里强烈安利一个图像方面大名鼎鼎的开源软件ImageMagick，这个软件有非常丰富的说明文档，包括详细的技术，特别值得一读。我曾经做过一些图像方面的研究，代码参考了很多ImageMagick源码。
例如将1.png、2.png、3.png合并成out.gif，每一张图片300ms，循环播放，可以使用如下命令
convert.exe -delay 30 1.png 2.png 3.png -loop 0 out.gif
其中-delay 30表示一帧300ms，-loop 0表示无限次循环播放。需要注意的是，一定要把ImageMagick添加到环境变量中，不然convert就找不到命令了。
