经典传统music算法

目录

1 传统music算法

2 协方差矩阵

3 协方差

4 期望值

5 方差

6 两个子空间之间的正交性

7 奇异值分解(SVD)

8 奇异值

9 主成分分析(PCA)

10 流型矩阵

11 MUSIC算法的基本步骤

12 在MATLAB中实现MUSIC算法步骤

13 信号模型

14 假设不同信号源的信号之间是相互正交的,噪声与信号之间是正交的,则阵列输出信号 𝑥(𝑡) 的协方差矩阵为

15 MUSIC 算法的谱估计公式

16 输出信号协方差矩阵的近似值

17 输出信号协方差矩阵R的特征分解

18 正规阵

19 正定矩阵或者半正定矩阵

20 厄米特矩阵(Hermite矩阵)

21 MATLAB代码示例


1 传统music算法

MUSIC算法(Multiple Signal Classification)是一种用于估计信号来波方向(DOA估计)的多信号分类算法。它的基本思想是将任意阵列输出数据的协方差矩阵进行特征值分解,对应不同特征值的特征向量构成相互正交的信号子空间和噪声子空间。大特征值对应的特征向量构成信号子空间,小特征值对应的特征向量构成噪声子空间。通过两个子空间之间的正交性来估计信号的方向。
在MUSIC算法中,我们首先构建流型矩阵,然后对协方差矩阵进行特征分解,得到信号子空间和噪声子空间。这使我们能够准确估计信号源的方向。

传统的MUSIC(Multiple Signal Classification)算法是一种用于方向估计(DOA)的空间谱估计方法。它基于信号的协方差矩阵进行特征值分解,将观测空间分解为信号子空间和噪声子空间。这两个子空间是正交的,信号子空间由与信号相关的特征向量组成,而噪声子空间则由与噪声功率相对应的最小特征值的特征向量组成。

当声源发声时,阵列会接收到来自目标方向的信号,但是也会接受到不同方向的反射信号。MUSIC算法可以剔除掉其余的反射信号,选出来自目标方向的那个信号,从而得到目标的方向。

MUSIC算法适用于来波为平行波,即目标与阵列的距离L远大于阵元之间间距d。此时来自目标的信号相对于每个阵元的方位角基本可视为相同。

天线阵元数,阵元间距,快拍数和信噪比都对MUSIC算法的性能影响很大。

信号到达阵元2多走了一段距离,那么信号到达阵元2的时间总会相比阵元1延迟\Delta t(我们通常称之为时延,实际上相位差就是时延导致的)

MUSIC算法的最终目的:计算\theta

可以知道\theta和两个阵元接收信号的相位差\phi紧密相关。能求到\phi,就能求得\theta

那么在理想条件,也就是没有任何的反射和折射,且只有一个声源,这时直接用两阵元接收信号相除y_{2}/y_{1}就能得到相位差\phi,从而得到目标方位角\theta

假设信号x_{1}和信号x_{2}是不相关(MUSIC算法的假设条件1)的(当信号x_{1}和信号x_{2}线性相关时,可以找到一个非零复数c,使得x_{1}=c*x_{2}

阵元个数 > 声源信号的个数(MUSIC算法的假设条件2)

MUSIC算法是通过求1/\begin{vmatrix} \vec c*\vec a \end{vmatirx}1/\begin{vmatrix} \vec c*\vec a \end{vmatirx}的最大值(谱峰搜索)来找相应的解\phi,对应的也就是相应的目标方位角。

2 协方差矩阵

协方差矩阵(covariance matrix)在统计学和概率论中是一个方阵,用于描述两列随机变量之间的协方差。它是协方差的直接推广。

定义:
假设有两列实数随机变量序列,分别记为 X 和 Y。
若二者对应的期望值分别为 μ_X 和 μ_Y,则这两列随机变量的协方差矩阵为: [ \Sigma = \begin{bmatrix} \sigma(X, X) & \sigma(X, Y) \ \sigma(Y, X) & \sigma(Y, Y) \end{bmatrix} ],Sigma(大写Σ,小写σ)。
其中,对角线上的元素是各个变量自身的方差,非对角线上的元素是两两随机变量之间的协方差。

这个形式的协方差矩阵是一个 2x2 的矩阵,其中包含了两个随机变量 X 和 Y 之间的协方差信息。让我们来看看每个元素的含义:

  1. [ \sigma(X, X) ]:这是 X 变量自身的方差,表示 X 的变化幅度。
  2. [ \sigma(Y, Y) ]:这是 Y 变量自身的方差,表示 Y 的变化幅度。
  3. [ \sigma(X, Y) ]:这是 X 和 Y 之间的协方差,表示二者的联合变化情况。

如果协方差矩阵是对称的,那么它的特征向量将构成信号子空间和噪声子空间,用于估计信号源的方向。

协方差矩阵反映了两个随机变量变化时是同向还是反向的(相关性):
若协方差 > 0,则说明这两个随机变量同向变化。
若协方差 < 0,则说明是反向变化。
若协方差 = 0,则说明是两个随机变量的变化方向没有任何相似度.
协方差矩阵在多元正态分布、线性变换等领域具有重要应用。

3 协方差

协方差(Covariance)用于衡量两个变量的总体误差。如果两个变量的变化趋势一致,即其中一个大于自身的期望值,另外一个也大于自身的期望值,那么两个变量之间的协方差就是正值。如果两个变量的变化趋势相反,即其中一个大于自身的期望值,另外一个却小于自身的期望值,那么两个变量之间的协方差就是负值。
协方差的绝对值越大,两个变量相互影响越大。总之,协方差是用来度量两个变量之间的“协同变异”大小的总体参数。

4 期望值

在概率论和统计学中,期望值(或数学期望、或均值,亦简称期望,物理学中称为期待值)是指在一个离散性随机变量试验中每次可能结果的概率乘以其结果的总和。换句话说,期望值是随机试验在同样的机会下重复多次的结果计算出的等同“期望”的平均值。
对于离散随机变量,期望值是所有可能结果的概率加权和。它是每个结果乘以其发生概率的总和。
例如,掷一枚公平的六面骰子,其每次「点数」的期望值是3.5。
期望值帮助我们理解和预测随机变量的平均行为,是概率分布中的一个重要特征。

5 方差

在概率论和统计学中,方差是与平均值的平方距离的平均值。它表示随机变量如何在均值附近分布。小方差表明随机变量分布在平均值附近,而大方差则表明随机变量的分布远离平均值。例如,对于正态分布,窄钟形曲线将具有较小的方差,而宽钟形曲线将具有较大的方差。

6 两个子空间之间的正交性

正交子空间:
如果两个子空间 V 和 W 中的每个向量都垂直于另一个子空间中的每个向量,那么这两个子空间是正交的。
换句话说,V 中的每个向量与 W 中的每个向量都正交,或者说它们的内积为零。
你可以想象在一个房间里,地面是一个子空间 V,两面墙的交线是另一个子空间 W,这两个子空间是正交的。
性质:
正交子空间具有以下性质:
行空间和零空间是正交的:如果矩阵 A 的每行与零空间向量 x 的内积都等于零,那么行空间与零空间正交。
列空间和左零空间是正交的:如果矩阵 A 的每列与左零空间向量 y 的内积等于零,那么列空间与左零空间正交。
行空间与列空间的秩相等,都等于矩阵的秩。
每个向量都可以由行空间和零空间的组分构成。
正交子空间在线性代数和数据科学中具有重要应用,例如奇异值分解(SVD)主成分分析(PCA)

7 奇异值分解(SVD)

奇异值分解(Singular Value Decomposition,缩写为SVD)是线性代数中一种重要的矩阵分解方法。给定一个矩阵 A,SVD将其分解为三个矩阵的乘积,即:
[ A = U \Sigma V^T ]
其中:
U 是一个 m×m 的矩阵。
Σ 是一个 m×n 的矩阵,除了主对角线上的元素以外全为零,主对角线上的每个元素都称为奇异值
V 是一个 n×n 的矩阵。
U 和 V 都是酉矩阵,即满足 U^TU = I 和 V^TV = I。奇异值分解在降维算法、推荐系统、自然语言处理等领域有广泛应用,是很多机器学习算法的基石。

8 奇异值

奇异值(Singular Value)是矩阵中的一个概念,通常通过奇异值分解来求得。设矩阵 A 为 m×n 阶矩阵,其中 q=min(m,n),则 A*A 的 q 个非负特征值的算术平方根即为 A 的奇异值。奇异值分解在信号处理、统计学等领域有广泛应用。

9 主成分分析(PCA)

主成分分析(PCA,Principal Component Analysis)是一种数学降维方法,用于将高维度的数据映射到较低维度的空间。让我来简单解释一下它的原理和作用:
原理:
PCA的目标是找到一组新的正交特征,称为主成分,以最大程度地保留原始数据的方差。
第一个主成分是原始数据中方差最大的方向,第二个主成分是与第一个正交的方向,依此类推。
通过这种方式,我们可以用较少的主成分来表示原始数据,实现降维。
作用:
数据降维:将高维数据转换为低维,减少计算复杂度。
去除噪音:滤去变化幅度较小的噪音,提高信噪比。
可视化:利用主成分展示多维数据的关系。
发现隐性相关变量:识别存在协同或拮抗关系的变量。
实现方法:
通过计算数据矩阵的协方差矩阵,然后得到协方差矩阵的特征值和特征向量,选择特征值最大的前几个特征向量作为主成分。
也可以使用奇异值分解(SVD)来实现PCA。

10 流型矩阵

流型矩阵(Manifold Matrix)是一个与李群和李代数相关的概念。
流型:
流型是一种数学结构,用于描述在局部上类似于欧几里德空间的空间。
李群的流型是一种特殊类型的流型,它在局部上具有李群的性质。
流型矩阵:
在李群的研究中,我们通常将李群的元素表示为矩阵。
流型矩阵是一种描述李群元素在流型上的变化的工具。
它可以用来映射李群的切空间到流型上的局部坐标系。
应用:
流型矩阵在李群表示、李代数、机器人学、SLAM(同时定位与地图构建)等领域有广泛应用。
通过流型矩阵,我们可以理解李群元素之间的关系以及它们在流型上的变化。

11 MUSIC算法的基本步骤

数据收集:收集由天线阵列接收到的信号数据。
协方差矩阵计算:根据收集到的数据计算信号的协方差矩阵。
特征值分解:对协方差矩阵进行特征值分解,得到信号子空间和噪声子空间。
构建谱函数:利用噪声子空间构建空间谱函数,该函数将用于确定信号源的方向。
谱峰搜索:通过搜索空间谱函数的峰值来估计信号源的方向。
结果分析:分析得到的方向估计结果,并进行必要的后处理。

12 在MATLAB中实现MUSIC算法步骤

定义信号源和阵列参数。
生成或收集信号数据。
计算信号的协方差矩阵。
进行特征值分解以分离信号子空间和噪声子空间。
构建并搜索空间谱函数以估计信号源的方向。
MUSIC算法的性能受多个因素影响,包括信号源数目、阵列的几何形状、信号与噪声的比例等。在实际应用中,可能需要根据具体情况调整算法的参数和实现细节。

13 信号模型

假设空间中存在 𝑀 个不同方向的信号,入射到由 𝑁 个天线单元构成的均匀直线阵上。第 𝑖 个信号源的方向为 𝜙𝑖(𝑖=1,…,𝑀),第 𝑖 个信号源的信号为 𝛼𝑖(𝑡)。假设 𝑀<𝑁 。

令第 𝑛 个天线单元的噪声为 𝑛𝑛(𝑡)。在窄带远场条件下,第 𝑛 个天线单元的输出信号 𝑥𝑛(𝑡) 可表示为:

\begin{gathered} x_{n}(t) =\sum_{i=1}^M\alpha_i(t)e^{jkz_n\sin\phi_i}+n_n(t) \\ =\sum_{i=1}^M\alpha_i(t)s_n(\phi_i)+n_n(t) \end{gathered}

将 𝑁 个天线的输出信号表示成向量形式 𝑥(𝑡),则上式可归纳为: 

\begin{gathered} \mathbf{x}(t)=\mathbf{S}\alpha(t)+\mathbf{n}(t) \\ \begin{bmatrix}x_1\left(t\right)\\x_2\left(t\right)\\\cdots\\x_N\left(t\right)\end{bmatrix}=\begin{bmatrix}&e^{jkz_1\sin\phi_1}&e^{jkz_1\sin\phi_2}&\ldots&e^{jkz_1\sin\phi_M}\\&e^{jkz_2\sin\phi_1}&e^{jkz_2\sin\phi_2}&\ldots&e^{jkz_2\sin\phi_M}\\&\cdots&\cdots&\ldots&\ldots\\&e^{jkz_N\sin\phi_1}&e^{jkz_N\sin\phi_2}&\ldots&e^{jkz_N\sin\phi_M}\end{bmatrix}\begin{bmatrix}\alpha_1\left(t\right)\\\alpha_2\left(t\right)\\\cdots\\\cdots\\\alpha_M(t)\end{bmatrix} \\ +\begin{bmatrix}n_1(t)\\n_2(t)\\\ldots\\n_N(t)\end{bmatrix} \end{gathered}

其中,𝑆 为阵列的流型矩阵,矩阵规模为 𝑁×𝑀,具体可表示为 𝑀 个不同方向对应的阵列导向矢量:

\mathbf{S}=[\mathbf{s}(\phi_1),\mathbf{s}(\phi_2),\ldots,\mathbf{s}(\phi_M)]

\varphi_\mathrm{k}=\frac{2\pi\mathrm{d}}\lambda\sin\theta_\mathrm{k}

由于 𝑀<𝑁,流型矩阵 𝑆 为列满秩矩阵,Rank(𝑆)=𝑀。

X_{N1}=A_{NM}S_{M1}+N_{N1}

其中X为阵列接收到的信号矩阵,两个维度分别代表:阵元个数(number of array elements)、采样点数(snapshots);A为阵列方向矩阵,两个维度分别代表:阵元个数、信号方向的方向向量;s为信号源发射信号矩阵,两个维度分别代表:信号源个数、采样点数;N为噪声矩阵,两个维度分别为阵元个数、采样点数。

14 假设不同信号源的信号之间是相互正交的,噪声与信号之间是正交的,则阵列输出信号 𝑥(𝑡) 的协方差矩阵为

\begin{aligned} \text{R}& =\operatorname{E}[\mathbf{x}(t)\mathbf{x}(t)^H] \\ &=\mathrm{E}[\mathbf{S}\alpha(t)\alpha(t)^H\mathbf{S}^H+\mathbf{n}(t)\mathbf{n}(t)^H] \\ &=\mathbf{SAS}^H+\sigma^2\mathbf{I} \\ &=\mathbf{R}_S+\sigma^2\mathbf{I} \end{aligned}

噪声是高斯白噪声(服从均值为0,方差为 σ^2) ,其中H表示矩阵的共轭转置。

σ ^2是噪声功率、I是M × M阶的单位矩阵。

其中,𝐴 为不同信号源之间的协方差矩阵,由于不同信号源之间是相互正交的,𝐴 为正定对角矩阵

\mathbf{A}=\begin{bmatrix}&\operatorname{E}[|\alpha_1(t)|^2]&\cdots&\cdots&\cdots\\&\cdots&\operatorname{E}[|\alpha_2(t)|^2]&\cdots&\cdots\\&\cdots&\cdots&\cdots&\cdots\\&\cdots&\cdots&\cdots&\operatorname{E}[|\alpha_M(t)|^2]\end{bmatrix}

15 MUSIC 算法的谱估计公式

\mathrm{A^Tq_n=0,n=M+1,\ldots,N}

信号来波方向的相位矢量A与噪声特征值对应的特征向量正交。对协方差矩阵R进行特征值分解,并将特征值从到到小排列,其中第M+1到第N个特征值对应的特征向量就是噪声子空间的基。定义一个噪声矩阵:

\mathrm E_n=[\mathrm q_{\mathrm M+1},\ldots,\mathrm q\mathrm N]

\mathrm{A^H(\theta_i)E_n^*E_n^TA(\theta_i)=0}

\mathrm{P_{music}~=1/(A^H(\theta_i)P_NA(\theta_i))}

\mathrm{P_{N}~=E_{n}^{*}E_{n}^{\mathrm{T}}}

P_\text{MUsIC}(\phi)=\frac{1}{\left|\left|\mathbf{Q}_n^H\mathbf{s}(\phi)\right|\right|^2}

 𝜙 与信号源方向相同时,分母为零,此时 MUSIC 谱估计为无穷大。因此,MUSIC 谱估计的尖峰数目与信源数目相同,尖峰对应的方向即为信号源的方向

16 输出信号协方差矩阵的近似值

\mathbf{R}=\frac{1}{K}\sum_{k=1}^K\mathbf{x}_k\mathbf{x}_k^H

17 输出信号协方差矩阵R的特征分解

定理1:方阵A酉相似于对角矩阵的充要条件是:A为正规阵(实或复)

定理2:对于实矩阵A,若满足\mathrm{A^{T}A=AA^{T}} ,则A为正规阵

对于复矩阵A,若满足\mathrm{A^{H}A}=\mathrm{AA^{H}},则A为正规阵

由于协方差矩阵为厄米矩阵(\mathrm{R_s^H~=~R_s~}),厄米矩阵是正规矩阵,所以协方差矩阵Rs酉相似于对角矩阵\mathrm{R_s=Q\Lambda Q^H}

\mathbb{Q}^\mathrm{H}\mathbb{Q}=\mathbb{I}\Lambda=\operatorname{diag}(\lambda_1,\lambda_2,\ldots,\lambda_\mathrm{L})

定理3:R为正定或者半正定矩阵,所有特征值为大于等于零的正实数

定理4:R不同特征值对应的特征向量相互正交,将一组完备的正交基分为信号子空间和噪声子空间

\mathbf{R}=\mathbf{Q}(\mathbf{\Lambda}+\sigma^2\mathbf{I})\mathbf{Q}^H

=\mathbf{Q}\begin{bmatrix}\lambda_1+\sigma^2&0&\ldots&0&0&\ldots&0\\0&\lambda_2+\sigma^2&\ldots&0&0&\ldots&0\\\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots\\0&0&\ldots&\lambda_M+\sigma^2&0&\ldots&0\\0&0&\ldots&0&\sigma^2&\ldots&0\\\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots\\0&0&\ldots&0&0&\ldots&\sigma^2\end{bmatrix}\mathbf{Q}^H

将输出信号矩阵 𝑅 进行特征分解,得到的 𝑁−𝑀 个较小且相等的特征值对应的特征向量即可构成 𝑄𝑛。其中,Λ的对角线上M个元素的正数,其余元素均为零。R的M个大的特征值对应的特征向量张成的子空间为信号子空间,其余的N+1-M个小特征值对应的特征向量张成的子空间为噪声子空间。

a^H(\theta)U_N=0

Un是由R的所有特征值中娇小的(阵元个数-信号源个数)个特征向量组成的子空间,称为噪声子空间。

要求理想情况下信号子空间和噪声子空间正交,也就是信号子空间中的方向向量a(theta)和噪声子空间正交:由于噪声的存在,实际上a(theta)和Un并不能完全正交。因此实际上是通过进行最小优化搜索来实现的:

\theta_{MUSIC}=argmin_{\theta}a^{H}(\theta)U_NU_N^{H}a(\theta)

由于实际中阵列接受数据是有限的,所以通常由协方差矩阵的最大似然估计来代替协方差矩阵:

R_x=\frac{1}{N}\sum_{i=1}^{N}X(i)X^H(i)

R x 是R 的最大似然估计。
\mathbb{N}\to\infty,他们是一致的,但实际情况将由于样本数有限而造成误差。

18 正规阵


正规阵,无论是实数矩阵还是复数矩阵,指的是与自己的共轭转置矩阵(对于实数矩阵则是转置矩阵)相乘的结果与其共轭转置矩阵相乘的结果相同的矩阵。数学上,如果矩阵 ( A ) 满足 ( AA^* = A^A ),其中 ( A^ ) 表示 ( A ) 的共轭转置,那么 ( A ) 就是一个正规矩阵。正规矩阵的一个重要性质是它们可以通过酉变换(对于实数矩阵是正交变换)对角化。
正规矩阵包括几个特殊类别的矩阵:
对角矩阵:所有非对角线元素都是零。
实对称矩阵:等于其转置的实数矩阵。
厄米特矩阵:等于其共轭转置的复数矩阵。
反厄米特矩阵:其共轭转置等于其本身取负的复数矩阵。
正交矩阵:矩阵的转置是其逆矩阵的实数矩阵。
酉矩阵:矩阵的共轭转置是其逆矩阵的复数矩阵。
这些矩阵都有一个共同的性质,即它们与自己的共轭转置(或转置)可交换

19 正定矩阵或者半正定矩阵

正定矩阵和半正定矩阵是线性代数中的重要概念,它们在多个领域如优化问题、统计学和机器学习中都有广泛应用。
正定矩阵的定义是:对于一个实对称矩阵 ( A ),如果对于所有非零实向量 ( x ),都有 ( x^T A x > 0 ),则称矩阵 ( A ) 为正定矩阵。这意味着正定矩阵的所有特征值都是正的。
半正定矩阵的定义是:对于一个实对称矩阵 ( A ),如果对于所有向量 ( x ),都有 ( x^T A x \geq 0 ),则称矩阵 ( A ) 为半正定矩阵。半正定矩阵包括了正定矩阵,因为它允许特征值为零。

20 厄米特矩阵(Hermite矩阵)

厄米特矩阵(Hermite矩阵)是一种特殊的复方阵,满足矩阵等于其自身的共轭转置。数学上,如果一个矩阵 ( A ) 满足 ( A = A^H ),其中 ( A^H ) 表示 ( A ) 的共轭转置(即矩阵的转置并取每个元素的共轭),那么 ( A ) 就是一个厄米特矩阵。厄米特矩阵的对角元素必须是实数,而且它的特征值都是实数。此外,厄米特矩阵的任意两个不同特征值对应的特征向量是正交的。
厄米特矩阵在物理学中尤为重要,因为它们经常出现在量子力学中,描述可观测量的算符。在数学和工程学中,它们也有广泛的应用,例如在信号处理中的协方差矩阵通常是厄米特的。

21 MATLAB代码示例

clc; clear; close all;
%MUSIC算法(Multiple Signal Classification)用于方向估计(Direction of Arrival, DoA)

%% 参数设置
%%% 工作频率
c = 3e8; %光速
freq = 10e9; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以调参 设置了工作频率为 10 GHz。在这个上下文中,它用于计算波长和波数,以便在MUSIC算法中进行方向估计。
lambda = c/freq;    % 波长
k = 2*pi/lambda;    % 波数 波数是波动的一种性质,表示每 2π 长度的波长数量(即每单位长度的波长数量乘以 2π)。更明确地说,波数是每 2π 长度内,波动重复的次数(一个波动取同样相位的次数)。

%%% 阵列参数
N = 12;                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以调参 阵元数量 阵列数量不能小于5
d = 0.5*lambda;        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以调参 阵元间隔 
z = (0:d:(N-1)*d)';     % 阵元坐标分布 在阵列信号处理中,我们通常需要知道每个阵元(传感器)的位置。z 是一个列向量,其中包含了阵列中每个阵元的坐标。0:d:(N-1)*d 表示从0开始,以步长 d 递增,直到 (N-1)*d。这样就生成了一系列坐标值,对应于阵列中不同阵元的位置。

%%% 信号源参数
phi = [-75,-40,35,90,-20, 50, -10, -30, 60]'*pi/180;   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以调参 来波方向  信号数量不能多于4个
M = length(phi);                % 信号源数目

%%% 仿真参数
SNR = 10;             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以调参% 信噪比(dB)
K = 1024;     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以调参% 采样点数

%% 阵列接收信号仿真模拟
S = exp(1j*k*z*sin(phi'));          % 流型矩阵 计算了给定阵列的流型矩阵。exp(1j*k*z*sin(phi')):这个表达式涉及几个组成部分:sin(phi'):计算到达角度(phi)的正弦值(以弧度为单位)。z:表示阵列元素的坐标。k:波数(与波长有关)。
%1j:表示虚数单位(即 √(-1))。
%exp(...):计算复指数。
%得到的矩阵 S 包含了与不同到达角度对应的复数相位。

Alpha = randn(M, K);         % 输入信号
%生成了一个随机输入信号矩阵 Alpha
%randn(M, K):这个函数生成一个大小为 M 行、K 列的随机矩阵,其中的元素是从标准正态分布(均值为0,方差为1)中随机抽取的。
%Alpha:这是一个包含随机输入信号的矩阵,用于模拟阵列接收信号。

X = S*Alpha;                        % 阵列接收信号、输出信号
X1 = awgn(X, SNR, 'measured');      % 加载高斯白噪声
%加载了高斯白噪声到阵列接收信号 X 
%awgn(X, SNR, 'measured'):这个函数将高斯白噪声添加到信号 X 中,以达到指定的信噪比(SNR)。
%SNR:信噪比(以分贝为单位)。
%X1:这是添加了高斯白噪声的阵列接收信号。

%% MUSIC 算法
%%% 阵列接收信号的协方差矩阵的特征分解
R = X1*X1'/K;    % 阵列接收信号的协方差矩阵 在阵列信号处理中,协方差矩阵描述了不同阵元之间的相关性。R = X1*X1'/K 计算了阵列接收信号的协方差矩阵,其中 X1 是加入了高斯白噪声的阵列接收信号。
[EV, D] = eig(R);       % 特征值分解 eig(R) 对协方差矩阵进行特征值分解,得到特征值和特征向量。EV 是特征向量构成的矩阵,按降序排列。
EVA = diag(D);          % 提取特征值 EVA 是提取的特征值。
[EVA, I] = sort(EVA, 'descend');   % 降序排序 对特征值进行降序排序。
Q = EV(:, I);           % 特征向量构成的矩阵EV EV(:, I):这一步从特征向量矩阵 EV 中选择了特定的列,这些列对应于排序后的特征值。Q:这是一个包含特征向量的矩阵,用于后续的MUSIC算法。
Q_n = Q(:, M+1:N);      % 噪声子空间 Q_n 是由特征向量构成的矩阵,其中包含了噪声子空间的信息。噪声子空间用于估计信号源方向时,将噪声成分从信号中分离出来。
%提取了噪声子空间 Q_n。
%Q(:, M+1:N):这一步从特征向量矩阵 Q 中选择了一部分列,对应于排序后的噪声子空间。
%Q_n:这是一个包含噪声子空间信息的矩阵,用于在估计信号源方向时将噪声成分与信号分离开来。

%%% 计算MUSIC谱估计函数
phi_list = linspace(-pi/2, pi/2, 200)'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以调参 % 创建了一个从 (-\frac{\pi}{2}) 到 (\frac{\pi}{2}) 的等间隔的向量,包含了200个元素。这些元素表示不同的方向角度,用于计算MUSIC谱估计函数。在MUSIC算法中,我们需要在一定范围内对信号源的方向进行估计。这里的 phi_list 就是用来表示不同方向的候选角度,以便计算MUSIC谱估计函数。
S1 = exp(1j*k*z*sin(phi_list'));    % 不同方向对应的流型矢量构成矩阵 建了一个流型矢量矩阵 S1,其中包含了不同候选角度对应的复数相位。
P_MUSIC = 1./sum(abs(Q_n'*S1).^2);     % MUSIC 谱估计公式计算了 MUSIC 谱估计。
%Q_n'*S1:这一步计算了噪声子空间和流型矢量构成矩阵之间的内积。
%abs(...).^2:计算了内积的模的平方。
%sum(...):对所有内积的模的平方进行求和。
%1./sum(...):取倒数,得到 MUSIC 谱估计。
%P_MUSIC:这是包含 MUSIC 谱估计值的向量。

%%% 转换为dB
P_MUSIC = abs(P_MUSIC);
P_MUSIC_max = max(P_MUSIC);
P_MUSIC_dB = 10*log10(P_MUSIC/P_MUSIC_max);

%%% 提取信号源方向
[P_peaks, P_peaks_idx] = findpeaks(P_MUSIC_dB);     % 提取峰值 MUSIC算法通过计算不同方向的伪谱估计来估计信号源的方向P_MUSIC_dB 是一个包含不同方向伪谱估计值的向量。findpeaks(P_MUSIC_dB) 用于从伪谱估计函数中找到峰值。P_peaks 是提取的峰值值。P_peaks_idx 是对应于这些峰值的索引。
[P_peaks, I] = sort(P_peaks, 'descend');    % 峰值降序排序 
P_peaks_idx = P_peaks_idx(I); %将 P_peaks 中的值按照降序排列。排序后,I存储了排序后的索引,以便在后续步骤中使用。
P_peaks = P_peaks(1:M);             % 提取前M个
P_peaks_idx = P_peaks_idx(1:M);
phi_e = phi_list(P_peaks_idx)*180/pi;   % 估计方向
disp('信号源估计方向为:'); %disp 是MATLAB中的一个函数,用于在命令行中显示文本消息。
disp(phi_e); %phi_e 是一个包含估计的信号源方向的向量。这行代码将 phi_e 的值显示在命令行中。

%%% 绘图
figure;
plot(phi_list*180/pi, P_MUSIC_dB, 'k', 'Linewidth', 1); %绘制了一个图形,其中包含了不同方向的伪谱估计值。phi_list*180/pi 是x轴上的角度值,以度为单位。P_MUSIC_dB 是y轴上的伪谱估计值,以分贝(dB)为单位。
xlabel('\phi方位角(deg)');
ylabel('Pseudo-spectrum幅值(dB)');
axis([-90, 90, -40, 40]); % 设置了x轴和y轴的范围。具体来说:x轴范围:从 -90 度到 90 度。y轴范围:从 -40 到 0db。这样设置坐标轴范围可以确保绘制的MUSIC谱估计函数在合适的区域内显示,使得图形更易于理解和分析。
grid on;
hold on;
plot(phi_e, P_peaks, 'r.', 'MarkerSize', 25);%在图形上绘制了一系列红色的点。phi_e 是估计的信号源方向的角度值。P_peaks 是对应于这些方向的伪谱估计峰值。
hold on;
for idx = 1:M %在图形上添加了文本标签。phi_e(idx) 是估计的信号源方向的角度值。P_peaks(idx) 是对应于这些方向的伪谱估计峰值。sprintf('%0.1f°', phi_e(idx)) 用于格式化文本,将角度值显示为带有一位小数的度数。
    text(phi_e(idx)+3, P_peaks(idx), sprintf('%0.3f°', phi_e(idx)));
end

参考文章:DoA 估计:多重信号分类 MUSIC 算法(附 MATLAB 代码)

参考文章:MUSIC算法原理以及详细推导

参考文章:关于MUSIC算法的知识点和MATLAB程序详细注释——学习笔记

参考文章:较为详细的MUSIC算法原理及MATLAB实现

  • 24
    点赞
  • 38
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值