当前位置: 首页 > article >正文

基于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,M1(θ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,N1(θ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),,sK1(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=1Jz(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 MK个较小的特征值构成的对角矩阵 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)} P2DMUSIC(θ,ϕ)=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软件上仿真的结果。


http://www.kler.cn/a/553624.html

相关文章:

  • 基于SpringBoot+vue粮油商城小程序系统
  • 汇能感知的光谱相机/模块产品有哪些?
  • 【机器学习】K折交叉验证(K-Fold Cross-Validation)
  • 网工项目实践2.4 北京公司安全加固、服务需求分析及方案制定
  • linux内核数据结构之哈希表
  • 【Cesium学习(十二)】Cesium常见问题整理总结
  • dockerfile2-15
  • Vue 3最新组件解析与实践指南:提升开发效率的利器
  • Excel如何给单元格填色,以及如何用Python 3实现单元格填色
  • deepseek帮我设计物理量采集单片机口保护电路方案
  • 市场波动中的数据分析与策略优化
  • Spring源码分析のBean创建流程(下)
  • 科普:Docker run的相关事项
  • Vue前端开发-Vant之Image组件
  • 【基于SprintBoot+Mybatis+Mysql】电脑商城项目之设置默认收货地址及删除收货地址
  • React 渲染 Flash 接口数据
  • AI技术生成【参考】
  • K8s学习总结
  • 智慧园区安全调度的重要性
  • Oppo手机投屏到Windows的两个方法,Windows系统自带投屏的替补!