MUSIC 与 ESPRIT 算法对比：阵列信号DOA估计的两种经典方法（含MATLAB代码）

MUSIC和ESPRIT算法是阵列信号DOA估计的两种经典高分辨率方法。MUSIC算法基于信号与噪声子空间正交性，通过谱峰搜索实现高精度角度估计；ESPRIT则利用阵列旋转不变特性，通过特征值分解直接计算角度，计算效率更高。文章首先建立了均匀线阵接收模型，详细推导了两种算法的数学原理和实现步骤，并提供了MATLAB仿真代码。结果表明，两种算法都能有效估计信号来向角，MUSIC具有更高的分辨率，而E

weixin_63763813

1253人浏览 · 2025-09-11 21:07:49

weixin_63763813 · 2025-09-11 21:07:49 发布

MUSIC 与 ESPRIT 算法对比：阵列信号DOA估计的两种经典方法（含MATLAB代码）

一、引言

在现代雷达、无线通信、声学定位等领域中，阵列信号处理扮演着至关重要的角色。其中，来波方向估计（Direction of Arrival, DOA） 问题是核心任务之一。通过估计信号源的入射角度，系统可以实现目标定位、波束形成、干扰抑制等关键功能。

传统的波束形成方法（如 Capon、CBF）在分辨率和稳健性上存在一定局限。为突破这一瓶颈，研究者提出了基于子空间分解的高分辨率 DOA 估计算法，其中最具代表性的就是 MUSIC（Multiple Signal Classification） 和 ESPRIT（Estimation of Signal Parameters via Rotational Invariance Techniques）。

MUSIC 算法 通过利用信号子空间与噪声子空间的正交性，构造空间谱函数，从而在真实来向角处获得清晰的谱峰，具有极高的分辨率。
ESPRIT 算法 则依赖于阵列的旋转不变特性，通过特征值分解直接估计来向角，避免了谱峰搜索，计算效率更高。

二、阵列信号接收模型

在讨论 MUSIC 与 ESPRIT 之前，我们首先需要明确 阵列接收模型。本文以均匀线阵（Uniform Linear Array, ULA）为例，建立其数学描述。

2.1 阵列接收信号模型

假设阵列共有 $M$ 个阵元，阵元间距为 $d$ ，有 $P$ 个窄带平面波从方向 ${θ1,θ2,…,θP}\{\theta_1, \theta_2, \dots, \theta_P\}$ 入射，接收信号可表示为：

$\mathbf{x}(t) = \mathbf{A}(\theta)\mathbf{s}(t) + \mathbf{n}(t)$

其中：

$x(t)∈CM×1\mathbf{x}(t) \in \mathbb{C}^{M\times 1}$ ：阵列接收信号向量
$A(θ)=[a(θ1),a(θ2),…,a(θP)]\mathbf{A}(\theta) = [\mathbf{a}(\theta_1), \mathbf{a}(\theta_2), \dots, \mathbf{a}(\theta_P)]$ ：阵列流形矩阵
$s(t)∈CP×1\mathbf{s}(t) \in \mathbb{C}^{P\times 1}$ ：入射信号向量
$n(t)∈CM×1\mathbf{n}(t) \in \mathbb{C}^{M\times 1}$ ：加性白噪声向量

2.2 阵列流形向量

对于均匀线阵（ULA），阵列流形向量为：

$\mathbf{a}(\theta) = \begin{bmatrix} 1 \\ e^{-j \frac{2\pi d}{\lambda}\sin\theta} \\ e^{-j \frac{2\pi d}{\lambda}2\sin\theta} \\ \vdots \\ e^{-j \frac{2\pi d}{\lambda}(M-1)\sin\theta} \end{bmatrix}$

其中 $λ=c/f0\lambda = c/f_0$ 为载波波长。

2.3 协方差矩阵

在快拍数为 $L$ 的情况下，阵列接收信号的协方差矩阵可表示为：

$\mathbf{R}_x = E[\mathbf{x}(t)\mathbf{x}^H(t)] = \mathbf{A}(\theta)\mathbf{R}_s \mathbf{A}^H(\theta) + \sigma_n^2\mathbf{I}_M$

其中：

$Rs=E[s(t)sH(t)]\mathbf{R}_s = E[\mathbf{s}(t)\mathbf{s}^H(t)]$ 为信号协方差矩阵
$σn2\sigma_n^2$ 为噪声方差
$IM\mathbf{I}_M$ 为 $M$ 阶单位阵

在实际仿真中，通常通过样本协方差来估计：

$R^x=1L∑t=1Lx(t)xH(t) \hat{\mathbf{R}}_x = \frac{1}{L}\sum_{t=1}^L \mathbf{x}(t)\mathbf{x}^H(t)$

三、MUSIC 算法原理

3.1 基本思想

MUSIC（Multiple Signal Classification，多重信号分类）是基于子空间分解的高分辨率 DOA 估计算法。其核心思想是：

信号协方差矩阵 $Rx\mathbf{R}_x$ 可通过特征分解分为 信号子空间 与 噪声子空间。
当扫描方向角 $θ\theta$ 等于真实来向角时，阵列流形向量 $a(θ)\mathbf{a}(\theta)$ 落在信号子空间内，因此与噪声子空间正交：
$\mathbf{a}^H(\theta) \mathbf{U}_n \approx 0$
基于这一正交性，可以构造空间谱函数，利用谱峰位置估计 DOA。

3.2 数学公式

设协方差矩阵特征分解为：

$\mathbf{R}_x = \mathbf{U}_s \mathbf{\Lambda}_s \mathbf{U}_s^H \;+\; \mathbf{U}_n \mathbf{\Lambda}_n \mathbf{U}_n^H$

其中：

$Us∈CM×P\mathbf{U}_s \in \mathbb{C}^{M \times P}$ ：信号子空间
$Un∈CM×(M−P)\mathbf{U}_n \in \mathbb{C}^{M \times (M-P)}$ ：噪声子空间

MUSIC 空间谱定义为：

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

在真实来向角 $θk\theta_k$ 处，分母最小，从而在谱函数中出现峰值。

3.3 算法流程

计算样本协方差矩阵 $R^x\hat{\mathbf{R}}_x$ 。
对 $R^x\hat{\mathbf{R}}_x$ 进行特征分解，得到信号子空间与噪声子空间。
构造 MUSIC 空间谱 $PMUSIC(θ)P_{MUSIC}(\theta)$ 。
在扫描角度网格中搜索谱峰，获得 DOA 估计。

四、ESPRIT 算法原理

4.1 基本思想

ESPRIT（Estimation of Signal Parameters via Rotational Invariance Techniques，旋转不变子空间法）是一种经典的高分辨率 DOA 估计算法。与 MUSIC 依赖谱峰搜索不同，ESPRIT 通过 阵列结构的旋转不变性（Shift-Invariance） 直接估计信号来向角，因而具有计算效率高、无需扫描的优势。

对于均匀线阵（ULA），如果将其划分为两个高度相关的子阵列，则信号子空间在这两个子阵之间存在如下关系：

$\mathbf{U}_{s2} = \mathbf{U}_{s1} \mathbf{\Phi}$

其中：

$Us1,Us2\mathbf{U}_{s1}, \mathbf{U}_{s2}$ 分别为信号子空间在两个子阵中的投影
$Φ\mathbf{\Phi}$ 为一个对角矩阵，其特征值与来向角直接相关

4.2 数学公式推导

对 $Φ\mathbf{\Phi}$ 做特征分解：

$\psi_k = e^{-j \frac{2\pi d}{\lambda} \sin\theta_k}$

由此可得 DOA 估计公式：

$\theta_k = \arcsin\!\left(-\frac{\lambda}{2\pi d}\arg(\psi_k)\right)$

ESPRIT 主要有两种实现方式：

LS-ESPRIT（最小二乘）
$\mathbf{\Phi}_{LS} = (\mathbf{U}_{s1}^\dagger)\mathbf{U}_{s2}$
TLS-ESPRIT（总最小二乘）
通过对 $[Us1,Us2][\mathbf{U}_{s1}, \mathbf{U}_{s2}]$ 做奇异值分解 (SVD)，得到更加稳健的估计。

4.3 算法流程

构造两个相邻子阵的接收数据 $X_1, X_2$
通过 SVD 提取信号子空间 $Us\mathbf{U}_s$
拆分为 $Us1,Us2\mathbf{U}_{s1}, \mathbf{U}_{s2}$
构造旋转矩阵 $Φ\mathbf{\Phi}$ （LS 或 TLS）
由 $ψk\psi_k$ 的相位直接计算 $θk\theta_k$

五、代码

在这里插入图片描述

close all; clear; clc;
rng(2);
%% ---------------- 参数定义 ----------------
c = physconst('LightSpeed'); % 光速
f0 = 10e9; % 载波频率 Hz
fg = 50e9; % 采样率（未显式用到）
BW = 1e8; % 带宽（未显式用到）
src_angle = deg2rad([-10, -20]); % 信源DOA (rad)
M = 32; % 阵元数
L = 128; % 快拍数
snr = 20; % 信噪比 dB
scan_angle = deg2rad(-60:0.1:60); % 扫描角度

%% ---------------- 参数计算 ----------------
lambda = c / f0; % 波长
d = lambda/2; % 阵元间距
P = numel(src_angle); % 信源数

%% ---------------- 回波生成 ----------------
% 阵列流形
i = (0:M-1).';
A = zeros(M,P);

for k = 1:P
A(:,k) = exp(-1j*2*pi*d/lambda*i*sin(src_angle(k)));
end

% 信号矩阵
S = (randn(P,L)+1j*randn(P,L))/sqrt(2);

% 噪声
sig_pow = mean(abs(S(:)).^2);
snr_lin = 10^(snr/10);
noise_pow = sig_pow/snr_lin;
N = sqrt(noise_pow/2)*(randn(M,L)+1j*randn(M,L));

% 接收信号
X = A*S + N;

% 协方差矩阵
Rxx = (X*X')/L;

%% ---------------- MUSIC 算法 ----------------
% 特征分解
[V,D] = eig((Rxx+Rxx')/2);
[~,idx] = sort(diag(D),'descend');
V = V(:,idx);
Un = V(:,P+1:end); % 噪声子空间

% 计算谱函数
P_music = zeros(1,numel(scan_angle));
for g = 1:numel(scan_angle)
a = exp(-1j*2*pi*d/lambda*i*sin(scan_angle(g)));
P_music(g) = 1/(a'*(Un*Un')*a);
end

P_music = abs(P_music)/max(abs(P_music));

%% ---------------- ESPRIT 算法 ----------------
% 子阵划分
X1 = X(1:M-1,:);
X2 = X(2:M,:);
[U,Sv,~] = svd([X1;X2]/sqrt(L),'econ');
Us = U(:,1:P);
Us1 = Us(1:M-1,:);
Us2 = Us(M:end,:);

% LS-ESPRIT
Phi_ls = pinv(Us1)*Us2;
psi_ls = eig(Phi_ls);
doa_ls = asin(-real(angle(psi_ls))*lambda/(2*pi*d)); 
doa_ls_deg = sort(rad2deg(doa_ls));

% TLS-ESPRIT
[~,~,Vt] = svd([Us1,Us2],'econ');
V12 = Vt(1:P,P+1:end);
V22 = Vt(P+1:end,P+1:end);
Phi_tls = -V12/V22;
psi_tls = eig(Phi_tls);
doa_tls = asin(-real(angle(psi_tls))*lambda/(2*pi*d)); 
doa_tls_deg = sort(rad2deg(doa_tls));

%% ---------------- 绘图：MUSIC谱 ----------------
figure('Color','w'); hold on; grid on; box on;
plot(rad2deg(scan_angle), P_music,'b','LineWidth',1.5);
xlabel('扫描角度 (°)'); ylabel('归一化MUSIC空间谱');
title(sprintf('MUSIC vs ESPRIT (M=%d, L=%d, SNR=%ddB)',M,L,snr));

% 标注真实DOA (竖虚线)
for k = 1:P
xline(rad2deg(src_angle(k)),'--k',...
sprintf('真实 %d°',round(rad2deg(src_angle(k)))),...
'LineWidth',1.2);
end

% 标注 ESPRIT-LS 结果 (红色星号)
plot(doa_ls_deg, ones(size(doa_ls_deg)), 'r*', 'MarkerSize',10,...
'DisplayName','ESPRIT-LS');

% 标注 ESPRIT-TLS 结果 (绿色圆点)
plot(doa_tls_deg, ones(size(doa_tls_deg))*0.95, 'go', 'MarkerSize',8,...
'LineWidth',1.5,'DisplayName','ESPRIT-TLS');
legend('MUSIC谱','真实DOA','ESPRIT-LS','ESPRIT-TLS','Location','best');

%% ---------------- 结果打印 ----------------
fprintf('\n=== DOA估计结果 (单位：度) ===\n');
fprintf('真实 DOA : %s\n',num2str(rad2deg(src_angle)));
fprintf('ESPRIT-LS估计 : %s\n',num2str(doa_ls_deg(:).'));
fprintf('ESPRIT-TLS估计: %s\n',num2str(doa_tls_deg(:).'));

六、总结

本文围绕 MUSIC 与 ESPRIT 两种经典 DOA 估计算法展开了系统介绍与对比。首先，我们从阵列信号接收模型出发，推导了子空间分解与旋转不变性两种基本原理。随后，分别阐述了 MUSIC 的空间谱构造方法以及 ESPRIT 的直接特征值估计方法，并通过 MATLAB 仿真进行了对比实验。