MUSIC算法相关原理知识（物理解读+数学推导+Matlab代码实现）

部分来自于网络教程，如有侵权请联系本人删除

教程链接：MUSIC算法的直观解释：1，MUSIC算法的背景和基础知识_哔哩哔哩_bilibili

MUSIC算法的直观解释：2，我对于MUSIC算法的理解_哔哩哔哩_bilibili

https://blog.csdn.net/zhangziju/article/details/100730081

一、MUSIC算法作用

MUSIC（Multiple Signal Classification），多重信号分类，是一类空间谱估计算法。其思想是利用接收数据的协方差矩阵(Rx)进行特征分解，分离出信号子空间和噪声子空间，利用信号方向向量与噪声子空间的正交性来构成空间扫描谱，进行全域搜索谱峰，从而实现信号的参数估计。

MUSIC算法通常被用来使用麦克风阵列进行声源定位。

例如当麦克风阵列放在一个房间中，房间中存在一个声源。当声源发声时，阵列会接收到来自目标方向的信号，但是也会接受到不同方向的反射信号。MUSIC算法可以剔除掉其余的反射信号，选出来自目标方向的那个信号，从而得到目标的方向。

声波是机械波，通常利用麦克风阵列接受后转换为电信号进行处理。当信号是电磁波时，例如wifi信号，我们拿天线阵列进行接收，此时仍然可以利用MUSIC算法把不同方向信号的角度计算出来。

二、MUSIC算法原理

MUSIC算法适用于来波为平行波，即目标与麦克风阵列的距离L远大于阵元之间间距d。此时来自目标的信号相对于每个阵元的方位角基本可视为相同。具体如下图所示：

图1

1.时延、相位差和目标方位角关系

假设信号源发射信号为 $x=e^{jf}$

当信号从声源目标S传播到阵元1时，信号传播了 $L_{1}$ 距离，假设声速为c，则耗时 $L_{1}/c$

这样便会导致阵元1接受信号 $y_{1}$ 的相位和发射信号 $x$ 不一致，会延迟 $e^{j 2\pi f(L_{1}/c)}$

那么阵元1最终接收信号为: $y_{1}=e^{jf}*e^{j2\pi f(L_{1}/c)}=x*e^{j2\pi f(L_{1}/c)}=x*\beta(\beta=e^{j2\pi f(L_{1}/c)})$

阵元2相对于阵元1多走了 $dcos\theta$ ，那么阵元2接受信号则为：

$y_{2}=e^{jf}*e^{j2\pi f((L_{1}+dcos\theta)/c)}=x*e^{j2\pi f((L_{1}+dcos\theta)/c)}=y_{1}*\phi_{1}=x*\beta*\phi_{1}$

其中 $\phi_{1}=e^{j2\pi f(dcos\theta/c)}$ ，代表阵元2接收信号和阵元1接收信号之间的相位差

阵元3相对于阵元1多走了 $2dcos\theta$ ，那么相应的相位差为： $\phi_{2}=e^{j2\pi f(2dcos\theta/c)}=\phi_{1}^{2}$

则阵元3接收信号为： $y_{3}=y_{1}*\phi_{1}^{2}=x*\beta*\phi_{1}^{2}$

PS：如果还不好理解因为多走了一段距离导致的相位差怎么计算，可以这样理解(以阵元1和2为例)：

假设阵元1接收信号为 $y_{1}=sin(wt)$

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

那么阵元2接收信号则为 $y_{2}=sin(w(t+\Delta t))$

很明显， $y_{1}$ 和 $y_{2}$ 的相位差为 $\phi=w(t+\Delta t)-wt=w\Delta t$

由数字信号处理知识可以知道 $w=2\pi f$ ，同时在图1中阵元1和阵元2的时延 $\Delta t=(dcos\theta)/c$

那么就可以得到相位差 $\phi=w\Delta t=2\pi f*(dcos\theta/c)$

2.MUSIC算法核心原理（思路来源）

MUSIC算法的最终目的：计算 $\theta$

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

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

但是实际上会有很多反射折射信号被麦克风阵列所接收，而且声源可能不止一个，此时该怎么办呢？这就是MUSIC算法需要解决的问题。

好，那么我们假设一共有两个声源A,B，发射信号分别为 $x_{1}$ 和 $x_{2}$ （暂时不考虑反射和折射）

那么在某时刻t，三个阵元接收信号分别为:

$y_{1}[t]=\beta _{1}x_{1}[t]+ \beta _{2}x_{2}[t]$

$y_{2}[t]=\beta _{1}x_{1}[t]\phi_{1}+ \beta _{2}x_{2}[t]\phi_{2}$

$y_{3}[t]=\beta _{1}x_{1}[t]\phi_{1}^{2}+ \beta _{2}x_{2}[t]\phi_{2}^{2}$

那么某段时间内，麦克风阵列所接收信号为:

$\begin{bmatrix} y_{1}[1]&y_{1}[2]&...&y_{1}[n]\\ y_{2}[1]&y_{2}[2]&...&y_{2}[n]\\ y_{3}[1]&y_{3}[2]&...&y_{3}[n] \end{bmatrix} = \begin{bmatrix} 1 &1\\ \phi_{1} & \phi_{2}\\ \phi_{1}^{2} & \phi_{2}^{2} \end{bmatrix}* \begin{bmatrix} \beta_{1}x_{1}[1]&\beta_{1}x_{1}[2]&...&\beta_{1}x_{1}[n]\\ \beta_{2}x_{2}[1]&\beta_{2}x_{2}[2]&...&\beta_{2}x_{2}[n] \end{bmatrix}$

对应可以写成: $Y=\Phi*X$

其中 $Y$ 我们已知， $\Phi$ 是需要求得的， $X$ 可能已知，可能不知（当 $X$ 已知且可逆时，可以直接用逆矩阵求 $\Phi=Y*X^{-1}$ ，但是这样的情况较少）

能不能通过处理 $Y$ ，将等式右边的 $X$ 消除掉？（MUSIC算法的核心）

怎么处理呢？？？

如果能找到三个复数 $c_{1}$ 、 $c_{2}$ 和 $c_{3}$ ，分别对3个阵元接收信号进行幅度和相位变换（用1个复数去乘以一个信号，则是对这个信号进行了幅度和相位变换），且变换后接收信号之间完全抵消了，即：

$y_{1}c_{1}+y_{2}c_{2}+y_{3}c_{3}=0$

或用矩阵表示：

$\begin{bmatrix} y_{1}[1]&y_{2}[1]&y_{3}[1]\\ y_{1}[2]&y_{2}[2]&y_{3}[2]\\ ...&...&... \\y_{1}[n]&y_{2}[n]&y_{3}[n] \end{bmatrix} * \begin{bmatrix} c_{1}\\c_{2}\\c_{3} \end{bmatrix}=\vec0$

把 $y_{1}$ 、 $y_{2}$ 和 $y_{3}$ 分别代入上式可以得到：

$(\beta _{1}x_{1}+ \beta _{2}x_{2})c_{1}+(\beta _{1}x_{1}\phi_{1}+ \beta _{2}x_{2}\phi_{2})c_{2}+(\beta _{1}x_{1}\phi_{1}^{2}+ \beta _{2}x_{2}\phi_{2}^{2})c_{3}=0$

即：

$(\beta_{1}c_{1}+\beta_{1}\phi_{1}c_{2}+\beta_{1}\phi_{1}^{2}c_{3})x_{1}+(\beta_{2}c_{1}+\beta_{2}\phi_{2}c_{2}+\beta_{2}\phi_{2}^{2}c_{3})x_{2}=0$

MUSIC算法在此时进行了一个假设，即假设信号 $x_{1}$ 和信号 $x_{2}$ 是不相关（MUSIC算法的假设条件1）的（当信号 $x_{1}$ 和信号 $x_{2}$ 线性相关时，可以找到一个非零复数 $c$ ，使得 $x_{1}=c*x_{2}$ ）

那么此时上面公式中 $x_{1}$ 和 $x_{2}$ 的系数必须都为0，即：

$\begin{Bmatrix} \beta_{1}c_{1}+\beta_{1}\phi_{1}c_{2}+\beta_{1}\phi_{1}^{2}c_{3}=0\\ \beta_{2}c_{1}+\beta_{2}\phi_{2}c_{2}+\beta_{2}\phi_{2}^{2}c_{3}=0 \end{matrix}$

上式中系数 $\beta_{1}$ 和 $\beta_{2}$ 可以直接消除掉，那么可以看出，只要找到 $c_{1}$ 、 $c_{2}$ 和 $c_{3}$ 便可求出 $\phi_{1}$ 和 $\phi_{2}$

那么现在问题就转换为了，如何找到这一组复数 $c_{1}$ 、 $c_{2}$ 和 $c_{3}$ ？？？

要能找到这一组复数，必须满足：阵元个数 > 声源信号的个数（MUSIC算法的假设条件2）

其实最终就是解：

$(c_{1},c_{2},c_{3})*(1,\phi_{1},\phi_{1}^{2})^{T}=0$ 和 $(c_{1},c_{2},c_{3})*(1,\phi_{2},\phi_{2}^{2})^{T}=0$

即 $\begin{vmatrix} \vec c*\vec a \end{vmatirx}=0$ ，那么MUSIC算法是通过求 $1/\begin{vmatrix} \vec c*\vec a \end{vmatirx}$ 的最大值（谱峰搜索）来找相应的解 $\phi$ ，对应的也就是相应的目标方位角 $\theta$ 。

3.MUSIC算法步骤总结

窄带远场信号的DOA数学模型为：

$X(t) = A(\theta)s(t)+N(t)$

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

那么阵列接收数据的协方差矩阵为：

$R = E[XX^H]=AE[SS^H]A^H+\sigma^2I=AR_{s}A^H+\sigma^2I$

由于信号和噪声互相独立，数据协方差矩阵可以分解为信号、噪声相关的两部分，其中Rs是信号的协方差矩阵，ARsA^H是信号部分。

对R进行特征分解有：

$R = U_{s}\Lambda_{s}U_{s}^{H}+U_{N}\Lambda_{N}U_{N}^{H}$

式中，Us是由R的所有特征值中较大的（信号源个数）个特征向量组成的子空间，称为信号子空间；Un是由R的所有特征值中娇小的（阵元个数-信号源个数）个特征向量组成的子空间，称为噪声子空间。

根据之前我们所推导的MUSIC算法的条件，要求理想情况下信号子空间和噪声子空间正交，也就是信号子空间中的方向向量a(theta)和噪声子空间正交：

$a^{H}(\theta)U_{N}=0$

由于噪声的存在，实际上a(theta)和Un并不能完全正交。因此实际上是通过进行最小优化搜索来实现的：

$\theta_{MUSIC} = argmin_{\theta}\ a^{H}(\theta)U_{N}U_{N}^{H}a(\theta)$

和我们上文所说一样，MUSIC实际上是通过谱峰搜索来求最优解theta：

$P_{MUSIC}=\frac{1}{a^{H}(\theta)U_{N}U_{N}^{H}a(\theta)}$

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

$\hat{R} = \frac{1}{L}\sum_{i=1}^{L}XX^H$

总结以上算法原理，MUSIC算法的步骤为：

1.根据N个接收信号矢量得到下面协方差矩阵的估计值：

$R_{x}=\frac{1}{N}\sum_{i=1}^{N}X(i)X^{H}(i)$

2.对第1步得到的协方差矩阵进行特征分解

$R_{x}=AR_{s}A^{H}+\sigma^{2}I$

3.矩阵 $R_{x}$ 会有M个特征值。将其从大到小进行排列： $\lambda_{1}>\lambda_{2}>...>\lambda_{M}>0$

其中D(D=信号源个数)个较大的特征值对应信号，将其对应的特征向量看做信号部分空间。

M-D(M=阵元个数)个较小的特征值对应噪声，将其对应的特征向量看做信号部分空间，得到噪声矩阵 $E_{n}$

4.使 $\theta$ 不断变化，计算谱函数：

$P(\theta)=\frac{1}{a^{H}(\theta)E_{n}E_{n}^{H}a(\theta)}$ 通过寻找谱峰来计算波达方向的估计值。此处的 $a(\theta)$ 为阵元的方向相应向量。

$a(\theta_{k})=[1,e^{-j\phi_{k}},...,e^{-j(M-1)\phi_{k}}]$ ， $\phi_{k}=\frac{2\pi d}{\lambda}sin\theta_{k}$

4.Matlab代码实现

clear all
close all
clc
%----------------均匀线列阵实现MUSIC算法------------------%
ang2rad = pi/180;                   % 角度转弧度系数
N = 10;                             % 阵元个数
M = 3;                              % 信源个数
theta = [-65,0,45];                 % 来波方向(角度)
snr = 10;                           % 信号信噪比dB
K = 512;                            % 总采样点
delta_d = 0.05;                     % 阵元间距
f = 2400;                           % 信号源频率
c = 340;                            % 声速d = 0:delta_d:(N-1)*delta_d;
A = exp(-1i*2*pi*(f/c)*d.'*sin(theta*ang2rad));   % 接收信号方向向量
S = randn(M,K);                     % 阵列接收到来自声源的信号
X = A*S;                            % 最终接收信号，是带有方向向量的信号
X1 = awgn(X,snr,'measured');        % 在信号中添加高斯噪声
Rx = X1*X1'/K;                      % 协方差矩阵
[Ev,D] = eig(Rx);                   % 特征值分解
% [V,D] = eig(A) 返回特征值的对角矩阵 D 和矩阵 V
% 其列是对应的右特征向量，使得 AV = VD
EVA = diag(D)';                     % 将特征值提取为1行
[EVA,I] = sort(EVA);                % 对特征值排序,从小到大。其中I为index：1,2，...，10
EV = fliplr(Ev(:,I));               % 对应特征矢量排序En = EV(:,M+1:N);                  % 取特征向量矩阵的第M+1到N列特征向量组成噪声子空间% 遍历所有角度，计算空间谱
for i = 1:361angle(i) = (i-181)/2;           % 映射到-90度到90度theta_m = angle(i)*ang2rad;a = exp(-1i*2*pi*(f/c)*d*sin(theta_m)).';p_music(i) = abs(1/(a'*En*En'*a));
end
p_max = max(p_music);
p_music = 10*log10(p_music/p_max);  % 归一化处理
figure()
plot(angle,p_music,'b-')
grid on
xlabel('入射角/度')
ylabel('空间谱/dB')

MUSIC算法相关原理知识（物理解读+数学推导+Matlab代码实现）相关推荐

前端面试：你应该了解的JS算法相关的知识
前端面试:你应该了解的JS算法相关的知识时间复杂度通常使用最差的时间复杂度来衡量一个算法的好坏. 常数时间 O(1) 代表这个操作和数据量没关系,是一个固定时间的操作,比如说四则运算. 对于一个算 ...
RFID相关原理知识及RFID应用领域详细介绍
RFID工作原理和相关原理知识其中,电子标签又称为射频标签.应答器.数据载体:阅读器又称为读出装置,扫描器.通讯器.读写器(取决于电子标签是否可以无线改写数据).电子标签与阅读器之间通过耦合组件实现 ...
原创 | 机器学习数学推导与代码实现30讲.pdf
机器学习 Author:louwill Machine Learning Lab 机器学习数学推导与代码实现30讲已完成,主要包括监督学习模型.无监督学习模型.集成学习模型和概率模型四个大类29个模型 ...
ACO蚁群算法优化BP神经网络（ACO-BPNN）回归预测MATLAB代码（有优化前后的对比）
ACO蚁群算法优化BP神经网络(ACO-BPNN)回归预测MATLAB代码(有优化前后的对比) 代码注释清楚. main为运行主程序,可以读取本地EXCEL数据. 很方便,容易上手. (以电厂运行数据 ...
PSO粒子群算法优化BP神经网络（PSO-BPNN）回归预测MATLAB代码
PSO粒子群算法优化BP神经网络(PSO-BPNN)回归预测MATLAB代码(有优化前后的对比) 代码注释清楚. main为运行主程序,可以读取本地EXCEL数据,也可以加载本地数据.mat,使用换自 ...
【MIMO】两种空间相关信道生成方式的记录（公式+MATLAB代码）
文章目录前言一. Kronecker相关信道模型二.生成方式1 1.公式 2.MATLAB代码三.生成方式2-complex correlation 1.公式 2.MATLAB代码四.生成方 ...
SSA麻雀算法优化BP神经网络（SSA-BPNN）回归预测MATLAB代码（有优化前后的对比
SSA麻雀算法优化BP神经网络(SSA-BPNN)回归预测MATLAB代码(有优化前后的对比) 代码注释清楚. main为运行主程序,可以读取本地EXCEL数据. 很方便,容易上手. (以电厂运行数据 ...
MATLAB应用实战系列NSGA-II多目标优化算法原理及应用实例（附MATLAB代码）
前言 NSGA-Ⅱ是最流行的多目标遗传算法之一,它降低了非劣排序遗传算法的复杂性,具有运行速度快,解集的收敛性好的优点,成为其他多目标优化算法性能的基准. NSGA-Ⅱ算法是 Srinivas 和 D ...
【滤波跟踪】基于Huber函数和最大相关熵的抗差滤波算法实现GNSS导航定位粗差处理附matlab代码
1 内容介绍雷达系统中跟踪滤波器的设计通常依赖于线性高斯系统.一旦系统为非线性且受到非高斯噪声干扰时,雷达跟踪性能便出现严重恶化.为了提高目标在非线性非高斯环境下跟踪的精度,将最大相关熵扩展卡尔曼滤 ...
滤波系列（一）卡尔曼滤波算法（KF）：详细数学推导
滤波系列(一)卡尔曼滤波算法(KF) 在本文,将给出卡尔曼滤波算法的详细数学推导过程,如果想直接了解卡尔曼滤波算法的应用,请看博客:卡尔曼滤波算法的应用(python代码)或者直接可以调用Python ...

MUSIC算法相关原理知识（物理解读+数学推导+Matlab代码实现）

MUSIC算法相关原理知识（物理解读+数学推导+Matlab代码实现）相关推荐

最新文章

热门文章