MUSIC 算法由 Schmidt 于 1979 年提出，是子空间类算法中最经典和广泛应用的一种。它通过空间谱估计实现高分辨率的 DOA（Direction of Arrival，波达方向）估计，尤其在多径效应等复杂场景中表现出色。

一、原理

MUSIC 的核心思想是利用接收信号的协方差矩阵分解出信号子空间和噪声子空间，并基于方向向量与噪声子空间的正交性构造空间谱函数，通过谱峰搜索估计 DOA。

1.信号模型

1.1接收信号的基本模型为：

$$
\mathbf{X}(t) = \mathbf{A} \mathbf{S}(t) + \mathbf{N}(t)
$$

- \(\mathbf{X}(t)\)：\(M \times 1\) 的接收信号向量，\(M\) 是阵元数。
- \(\mathbf{A}\)：\(M \times D\) 的阵列流形矩阵，\(D\) 是信号源数，每列 \(\mathbf{a}(\theta_i)\) 是方向向量。
- \(\mathbf{S}(t)\)：\(D \times 1\) 的信号向量，表示 \(D\) 个信号源的波形。
- \(\mathbf{N}(t)\)：\(M \times 1\) 的噪声向量，假设为零均值白噪声，方差为 \(\sigma^2\)。  

假设噪声与信号不相关，且信号之间在理想情况下不相关（多径环境下可能相干）。

1.2阵列流形矩阵 \(\mathbf{A}\)

阵列流形矩阵 \(\mathbf{A}\) 定义为：

$$
\mathbf{A} = [\mathbf{a}(\theta_1), \mathbf{a}(\theta_2), \dots, \mathbf{a}(\theta_D)]
$$

对于均匀线阵，方向向量 \(\mathbf{a}(\theta)\) 表示为：

$$
\mathbf{a}(\theta) = [1, e^{-j 2\pi \frac{d}{\lambda} \sin\theta}, e^{-j 4\pi \frac{d}{\lambda} \sin\theta}, \dots, e^{-j 2\pi (M-1) \frac{d}{\lambda} \sin\theta}]^T
$$

- \(d\)：阵元间距（通常 \(d = \lambda/2\) 以避免空间混叠）。  
- \(\lambda\)：信号波长。  
- \(\theta\)：信号的 DOA。  
- \(e^{-j 2\pi k \frac{d}{\lambda} \sin\theta}\)：第 \(k\) 个传感器的相位延迟，来源于波程差 \(k d \sin\theta\)。

注：相位延迟：

对于第 \(k\) 个传感器（位置 \(x_k = (k-1)d\)），波程差是 \((k-1) d \sin\theta\)。

相位差为：

$$
\phi_k = 2\pi \frac{(k-1) d \sin\theta}{\lambda}
$$

信号是复数形式，延迟表现为相位乘以 \(e^{-j\phi_k}\)。

2.协方差矩阵与特征分解

理论协方差矩阵定义为：

$$
\mathbf{R}_x = E[\mathbf{X}(t) \mathbf{X}^H(t)] = \mathbf{A} \mathbf{R}_s \mathbf{A}^H + \sigma^2 \mathbf{I}
$$

- \(\mathbf{R}_s = E[\mathbf{S}(t) \mathbf{S}^H(t)]\)：信号协方差矩阵。
- \(\sigma^2 \mathbf{I}\)：噪声协方差矩阵，\(\mathbf{I}\) 是单位矩阵。

实际中，使用有限快拍数 \(N\) 估计样本协方差矩阵：

$$
\hat{\mathbf{R}}_x = \frac{1}{N} \mathbf{X} \mathbf{X}^H
$$

对 \(\hat{\mathbf{R}}_x\) 进行特征值分解：

$$
\hat{\mathbf{R}}_x = \mathbf{U} \mathbf{\Lambda} \mathbf{U}^H
$$

- \(\mathbf{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_M)\)：特征值矩阵，\(\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_M\)。  
- \(\mathbf{U} = [\mathbf{u}_1, \mathbf{u}_2, \dots, \mathbf{u}_M]\)：特征向量矩阵。

分解结果：

- 信号子空间：由前 \(D\) 个较大特征值对应的特征向量组成，\(\mathbf{U}_s = [\mathbf{u}_1, \mathbf{u}_2, \dots, \mathbf{u}_D]\)，维度 \(M \times D\)。  
- 噪声子空间：由后 \(M - D\) 个较小特征值对应的特征向量组成，\(\mathbf{U}_n = [\mathbf{u}_{D+1}, \mathbf{u}_{D+2}, \dots, \mathbf{u}_M]\)，维度 \(M \times (M - D)\)。

注：确定信号源数 \(D\)

- 理论：\(D\) 已知。  
- 实际：通过特征值分布估计：
  - 观察 \(\lambda_i\) 的跳变点。  
  - 或使用信息准则（如 AIC）：

$$
\text{AIC}(k) = -2N(M - k) \log\left(\frac{\prod_{i=k+1}^M \lambda_i^{1/(M-k)}}{\frac{1}{M-k} \sum_{i=k+1}^M \lambda_i}\right) + 2k(2M - k)
$$

\(k\) 取最小化 AIC 的值作为 \(D\)。

3.空间谱函数

MUSIC 利用方向向量 \(\mathbf{a}(\theta)\) 与噪声子空间 \(\mathbf{U}_n\) 的正交性，构造空间谱函数：

$$
P_{\text{MUSIC}}(\theta) = \frac{1}{\mathbf{a}^H(\theta) \mathbf{U}_n \mathbf{U}_n^H \mathbf{a}(\theta)}
$$

- 当 \(\theta = \theta_i\)（真实信号方向）时，\(\mathbf{a}(\theta_i)\) 与 \(\mathbf{U}_n\) 正交，分母接近 0，\(P_{\text{MUSIC}}(\theta)\) 出现谱峰。  
- \(\mathbf{a}^H(\theta)\) 是 \(\mathbf{a}(\theta)\) 的共轭转置。

4.谱峰搜索

通过扫描所有可能的 \(\theta\)（通常为 \([-90^\circ, 90^\circ]\)），找到 \(P_{\text{MUSIC}}(\theta)\) 的峰值位置，即为 DOA 估计。

二、仿真代码

import numpy as np
import matplotlib.pyplot as plt
# 参数设置
M = 8          # 阵元数
D = 2          # 信号源数
N = 100        # 快拍数
d = 0.5        # 阵元间距（归一化波长）
theta = [20, 30]  # 真实DOA（度）
snr = 10       # 信噪比 (dB)

# 生成阵列流形向量
def steering_vector(theta, M, d):
    theta_rad = np.deg2rad(theta)
    k = np.arange(M)
    return np.exp(-1j * 2 * np.pi * d * k * np.sin(theta_rad))

# 生成信号
A = np.array([steering_vector(t, M, d) for t in theta]).T  # M x D
S = np.random.randn(D, N) + 1j * np.random.randn(D, N)     # 信号
noise = (10**(-snr/20)) * (np.random.randn(M, N) + 1j * np.random.randn(M, N))  # 噪声
X = A @ S + noise  # 接收数据

# 1. 计算协方差矩阵
R_x = (1/N) * (X @ X.conj().T)

# 2. 特征值分解
eigenvalues, U = np.linalg.eigh(R_x)
idx = np.argsort(eigenvalues)[::-1]  # 降序排列
eigenvalues = eigenvalues[idx]
U = U[:, idx]

# 3. 分离子空间（假设D已知）
U_s = U[:, :D]      # 信号子空间
U_n = U[:, D:]      # 噪声子空间

# 4. MUSIC谱函数
theta_range = np.linspace(-90, 90, 1000)  # 扫描角度
P_music = []
for th in theta_range:
    a = steering_vector(th, M, d)
    P_music.append(1 / np.abs(a.conj().T @ U_n @ U_n.conj().T @ a))
P_music = np.array(P_music)

# 5. 寻找谱峰
peaks = np.where((P_music[:-2] < P_music[1:-1]) & (P_music[1:-1] > P_music[2:]))[0] + 1
estimated_doa = theta_range[peaks]
print("Estimated DOA:", estimated_doa)

# 可视化
plt.plot(theta_range, 10*np.log10(P_music))
plt.xlabel("Angle (degrees)")
plt.ylabel("MUSIC Spectrum (dB)")
plt.title("MUSIC Algorithm DOA Estimation")
plt.grid()
plt.show()