基于MATLAB的均匀面阵MUSIC算法DOA估计仿真
基于MATLAB的均匀面阵MUSIC算法DOA估计仿真
文章目录
- 前言
- 一、二维MUSIC算法原理
- 二、二维MUSIC算法MATLAB仿真
- 三、MATLAB源代码
- 总结
前言
\;\;\;\;\; 在波达角估计算法中,MUSIC 算法与ESPRIT算法属于特征结构子空间算法,是波达角估计算法中的基石。在前面的文章 一文读懂MUSIC算法DOA估计的数学原理并仿真 中详细介绍了一维MUSIC算法即线阵MUSIC算法DOA估计的原理及仿真,本文将介绍二维MUSIC算法即均匀面阵的MUSIC算法DOA估计原理及MATLAB仿真。
提示:以下是本篇文章正文内容,尊重版权,引用请附上链接。
一、二维MUSIC算法原理
下图为面阵入射信号模型,
\;\;\;\;\;
假设从远场有
K
K
K 个互不相关的窄带信号,入射到一个阵元个数为
M
×
N
M×N
M×N 的平面阵列上。记第
i
i
i个入射信号的方位角和俯仰角分别为
θ
i
\theta_i
θi和
φ
i
\varphi_i
φi ,则阵列接收信号可以表示为:
z
(
t
)
=
A
s
(
t
)
+
n
(
t
)
\boldsymbol{z}(t)=\boldsymbol A \boldsymbol s(t)+\boldsymbol n(t)
z(t)=As(t)+n(t)其中
A
\boldsymbol A
A是维度为(MN×K)的均匀矩形阵列的阵列流形,可以表示为如下所示的式子:
A
=
[
a
(
θ
k
,
φ
1
)
,
a
(
θ
2
,
φ
2
)
,
⋯
,
a
(
θ
K
,
φ
K
)
]
T
\mathbf{A}=\begin{bmatrix}\boldsymbol{a}(\theta_k,\varphi_1),\boldsymbol{a}(\theta_2,\varphi_2),\cdots,\boldsymbol{a}(\theta_K,\varphi_K)\end{bmatrix}^T
A=[a(θk,φ1),a(θ2,φ2),⋯,a(θK,φK)]T
a
(
θ
k
,
φ
k
)
\boldsymbol{a}(\theta_k,\varphi_k)
a(θk,φk)为第k个入射信号的导向矢量,仅仅由阵列的阵元排布和参考阵元的选择所决定,用公式可以表示为:
a
(
θ
k
,
φ
k
)
=
a
x
(
θ
k
,
φ
k
)
⊗
a
y
(
θ
k
,
φ
k
)
∈
C
M
N
×
1
\boldsymbol{a}(\theta_k,\varphi_k)=\boldsymbol{a}_x(\theta_k,\varphi_k)\otimes\boldsymbol{a}_y(\theta_k,\varphi_k)\in C^{MN\times1}
a(θk,φk)=ax(θk,φk)⊗ay(θk,φk)∈CMN×1 其中
⊗
\otimes
⊗表示的是克罗内克内积(Kronecker Product),
a
x
(
θ
k
,
φ
k
)
\boldsymbol{a}_x(\theta_k,\varphi_k)
ax(θk,φk)表示x轴方向上均匀线阵接收信号的方向矢量,
a
y
(
θ
k
,
φ
k
)
\boldsymbol{a}_y(\theta_k,\varphi_k)
ay(θk,φk)表示y轴方向上均匀线阵接收信号的方向矢量,可分别写为如下数学表达式:
a
x
(
θ
k
,
φ
k
)
=
[
a
x
,
0
(
θ
k
,
φ
k
)
,
a
x
,
1
(
θ
k
,
φ
k
)
,
⋯
,
a
x
,
M
−
1
(
θ
k
,
φ
k
)
]
T
\boldsymbol{a}_x(\theta_k,\varphi_k)=\begin{bmatrix}a_{x,0}(\theta_k,\varphi_k),a_{x,1}(\theta_k,\varphi_k),\cdots,a_{x,M-1}(\theta_k,\varphi_k)\end{bmatrix}^T
ax(θk,φk)=[ax,0(θk,φk),ax,1(θk,φk),⋯,ax,M−1(θk,φk)]T
a
y
(
θ
k
,
φ
k
)
=
[
a
y
,
0
(
θ
k
,
φ
k
)
,
a
y
,
1
(
θ
k
,
φ
k
)
,
⋯
,
a
y
,
N
−
1
(
θ
k
,
φ
k
)
]
T
\boldsymbol{a}_y(\theta_k,\varphi_k)=\begin{bmatrix}a_{y,0}(\theta_k,\varphi_k),a_{y,1}(\theta_k,\varphi_k),\cdots,a_{y,N-1}(\theta_k,\varphi_k)\end{bmatrix}^T
ay(θk,φk)=[ay,0(θk,φk),ay,1(θk,φk),⋯,ay,N−1(θk,φk)]T 式中的
s
(
t
)
\mathbf{s}(t)
s(t)是信号源矢量,
n
(
t
)
\mathbf{n}(t)
n(t)为高斯白噪声矢量,服从
N
(
0
,
σ
2
)
N(0,\sigma^2)
N(0,σ2)分布,可以分别表示如下式子:
s
(
t
)
=
[
s
0
(
t
)
,
s
1
(
t
)
,
⋯
,
s
K
−
1
(
t
)
]
T
\mathbf{s}(t)=\left[\mathbf{s}_0(t),\mathbf{s}_1(t),\cdots,\mathbf{s}_{K-1}(t)\right]^T
s(t)=[s0(t),s1(t),⋯,sK−1(t)]T
n
(
t
)
=
[
n
0
(
t
)
,
n
1
(
t
)
,
⋯
,
n
M
N
(
t
)
]
T
\mathbf{n}(t)=\left[\mathbf{n}_0(t),\mathbf{n}_1(t),\cdots,\mathbf{n}_{MN}(t)\right]^T
n(t)=[n0(t),n1(t),⋯,nMN(t)]T
\;\;\;\;\;
阵列接收信号的协方差矩阵可以表示为:
R
=
E
[
z
z
H
]
\mathbf{R} = \mathbb{E}[\mathbf{z}\mathbf{z}^H]
R=E[zzH]
=
A
E
[
s
s
H
]
A
H
+
σ
2
I
= \mathbf A\mathbb{E}[\mathbf{s}\mathbf{s}^H]\mathbf A^H + \sigma^2\mathbf{I}
=AE[ssH]AH+σ2I
=
A
R
S
A
H
+
σ
2
I
=\mathbf A \mathbf R_S\mathbf A^H + \sigma^2\mathbf{I}
=ARSAH+σ2I 其中
R
S
\mathbf{R}_S
RS表示入射信号的协方差矩阵,
σ
2
I
\sigma^2\mathbf{I}
σ2I表示功率为
σ
2
\sigma^2
σ2的高斯白噪声的协方差矩阵。
\;\;\;\;\;
实际应用中天线阵列获取的信息是有限次的快拍,因此只能得到协方差矩阵的估计值
R
^
\hat{\mathbf{R}}
R^,其计算公式如下:
R
^
=
1
J
∑
j
=
1
J
z
(
j
)
z
H
(
j
)
\hat{\mathbf{R}} = \frac{1}{J}\sum_{j=1}^{J}\mathbf{z}(j)\mathbf{z}^H(j)
R^=J1j=1∑Jz(j)zH(j)
\;\;\;\;\;
由于接收信号的协方差矩阵
R
\mathbf{R}
R是对称矩阵,因此可以对其进行特征值分解,可以得到:
R
=
U
Λ
U
T
\mathbf{R} = \mathbf{U}\boldsymbol{\Lambda}\mathbf{U}^T
R=UΛUT 其中
U
\mathbf{U}
U为
R
\mathbf{R}
R的特征向量构成的矩阵,
Λ
\boldsymbol{\Lambda}
Λ是一个由特征值构成的对角矩阵。
Λ
=
d
i
a
g
{
λ
1
,
λ
2
,
.
.
.
,
λ
M
N
}
\boldsymbol{\Lambda} = diag\{ \lambda_1,\lambda_2,...,\lambda_{MN} \}
Λ=diag{λ1,λ2,...,λMN}
\;\;\;\;\;
假设对角矩阵中的特征值降序排列,满足如下关系:
λ
1
≥
λ
2
≥
⋯
≥
λ
K
>
λ
K
+
1
=
⋯
=
λ
M
N
=
σ
2
\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_K > \lambda_K + 1 = \cdots = \lambda_{MN} = \sigma^2
λ1≥λ2≥⋯≥λK>λK+1=⋯=λMN=σ2 由前
K
K
K个较大的特征值构成的对角矩阵
Λ
S
\boldsymbol{\Lambda}_S
ΛS,其对应的特征向量构成的矩阵
U
S
\mathbf U_S
US为信号子空间。由后
M
−
K
M-K
M−K个较小的特征值构成的对角矩阵
A
N
\mathbf A_N
AN,其对应的特征向量构成的矩阵
U
N
\mathbf U_N
UN为噪声子空间。
\;\;\;\;\;
根据前文假设,信号与噪声相互独立,因此信号子空间与噪声子空间是相互正交的,故信号阵列流矢量与噪声子空间也具有正交性。同一维MUSIC算法一样,可构造二维空间谱函数:
P
2
D
−
M
U
S
I
C
(
θ
,
ϕ
)
=
1
a
H
(
θ
,
ϕ
)
U
N
U
N
H
a
(
θ
,
ϕ
)
P_{2D-MUSIC}(\theta, \phi) = \frac{1}{\mathbf a^{H}(\theta, \phi) \mathbf U_N \mathbf U_N^{H} \mathbf a(\theta, \phi)}
P2D−MUSIC(θ,ϕ)=aH(θ,ϕ)UNUNHa(θ,ϕ)1
\;\;\;\;\;
当天线阵列的方向矢量与噪声子空间近似正交时,上式分母部分取极小值,空间谱函数在此时取得极大值,得到空间谱的谱峰。对空间谱进行谱峰搜索,就能够得到入射信号的方位角与俯仰角的角度,至此完成了对于信源的二维 DOA估计。
二、二维MUSIC算法MATLAB仿真
\;\;\;\;\; 参数设置如下:改变任何一个参数,仿真结果都会跟着改变,可以通过修改参数观察不同条件对估计结果的影响。
M=3; % x轴阵元个数
N=2; % y轴阵元个数
K=1024; % 快拍数
fc=100e+6; % 载波
fs=300e+6; % 采样频率
Pn=1; % 噪声功率
fines=[45 180 250 300]; % 信号入射方位角
thetas=[5 30 55 75]; % 信号入射俯仰角
signal_f=[15e6 30e6 45e6 60e6]; % 信号频率
signal_SNR=[30 30 30 30]; % 信噪比
m=(0:M-1)'; % x轴坐标
n=(0:N-1)'; % y轴坐标
c=3e+8; % 光速
lamda=c/fc; % 波长
dx=1/2*lamda; % x轴阵元间距
dy=1/2*lamda; % y轴阵元间距
\;\;\;\;\;
通过观察参数,可以得出以下结论,可以自己通过改变参数来验证,这里就不贴图了。
1、随着阵元数目的增大,MUSIC 算法的分辨率逐渐增强。
2、随着信号信噪比的增大,MUSIC 算法的分辨率逐渐增强。
3、当阵元间距与波长的比值为二分之一时,MUSIC算法能够有效进行 DOA 估计;当阵元间距小于波长的二分之一时,MUSIC 算法的分辨率会降低;当阵元间距大于波长的二分之一时,由于采样严重不足,MUSIC算法可能会丧失分辨能力。
三、MATLAB源代码
均匀面阵MUSIC算法DOA估计MATLAB仿真源代码
总结
\;\;\;\;\; 以上就是今天记录的所有内容,分享了均匀面阵MUSIC算法DOA估计的原理及其在MATLAB软件上仿真的结果。