DOA估计算法解析：MUSIC、ESPRIT与压缩感知的MATLAB实现

简介：一套面向信号处理与阵列测向学习者的MATLAB源码包，聚焦DOA（到达方向）估计中的三类典型方法——MUSIC、ESPRIT与压缩感知。MUSIC基于子空间正交性搜索谱峰，ESPRIT利用阵列的旋转不变性，而压缩感知借助角度域稀疏性重构来波方向，三者分别展示了不同技术路线的代码实现，适用于无线通信、雷达及声学场景中的测向实验。压缩包采用zip格式，共包含3个.m文件，分别对应上述三种算法的核心脚本，整体仅4KB，轻量易读，便于在MATLAB中直接运行或修改参数对比效果。目前已有1474人学习下载，适合需要快速上手DOA算法并验证原理的初学者，也可作为课程设计或课题预研的参考基底。通过这三个脚本，可直观观察子空间类与稀疏重构类方法的输出差异，节省从零搭建实验环境的时间，为后续算法选型提供实际依据。

1. DOA估计到底是什么：从阵列信号到角度谱

DOA（Direction of Arrival）估计在雷达、5G定位和声学监测里解决同一个问题：接收阵列只有几个阵元的采样电压，怎么把来波方向算出来。传统波束形成把能量扫一圈也能出峰，但分辨率受限于阵列孔径；而MUSIC、ESPRIT这类子空间算法基于一个反直觉的事实——真正有用的不是信号本身，而是“没有信号的那个子空间”。这份MATLAB源码集（MUSIC.m、esprit.m、CS.m）恰好把三种典型路线放在一起，适合刚接触DOA的通信、雷达方向学生，也适合要把算法移植到自研仿真平台的工程师。先搞清楚原理，再对照代码调参，最后处理相干源和网格失配，这可能比单纯跑通脚本省下更多时间。

2. MUSIC算法：子空间正交与谱峰搜索，附MATLAB实现

2.1 信号模型与协方差矩阵：特征分解为何能拆出噪声子空间

一个M元均匀线阵（ULA），来波方向为θk，导向矢量写成a(θk)。M个阵元的一次快拍可以写成：

X(t) = A S(t) + N(t)

其中A是M×K的导向矢量矩阵，第k列由阵元间相位差exp(-j2πd sinθk / λ)累积而成。阵元间距d与波长λ满足d≤λ/2时，相位差和角度才是一一对应的。接收数据的协方差矩阵Rxx = E[XX^H]，如果信号彼此独立，信号部分A Rs A^H是秩K的Hermitian矩阵，噪声部分是σ²I，所以Rxx的特征值分布会呈现K个大特征值和M-K个接近σ²的小特征值。与小特征值对应的特征向量张成噪声子空间Un，因为信号子空间和噪声子空间在理想条件下正交，任何真实导向矢量a(θk)投影到Un的范数都接近零，而非真实方向时投影不为零。MUSIC谱就是构造这个投影的倒数，在真实角度处形成尖峰。

实际工程中拿不到统计协方差，只能用N个快拍估计样本协方差Rxx = X X^H / N。快拍数不足时，特征值会发生泄漏，小特征值不再严格等于σ²，这是后面避坑章节要处理的第一类问题。

2.2 MUSIC谱计算与MATLAB代码

下面这段脚本是可运行的MUSIC基线版本，处理8元均匀线阵、两个非相干信号源：

% music_doa_demo.m 针对8元均匀线阵、两个独立源 M = 8; % 阵元数 d_lambda = 0.5; % 阵元间距/波长 theta_true = [-10 20]; % 真实角度（度） K = numel(theta_true); % 信号个数，先手工给 N = 500; % 快拍数 SNR_dB = 10; % 信噪比（dB） angs = theta_true*pi/180; A = exp(-1j*2*pi*d_lambda*(0:M-1)'*sin(angs)); % 正交复高斯源，功率归一化 S = (randn(K,N) + 1j*randn(K,N))/sqrt(2); sigma_n = sqrt(10^(-SNR_dB/10)); X = A*S + sigma_n*(randn(M,N) + 1j*randn(M,N))/sqrt(2); % 样本协方差矩阵 Rxx = (X*X')/N; % 特征分解，按特征值降序排列 [V,D] = eig(Rxx); [ev, idx] = sort(diag(D),'descend'); V = V(:,idx); Un = V(:, K+1:M); % 噪声子空间 % 谱搜索：-90度到90度，步长0.1度 theta_grid = -90:0.1:90; Pm = zeros(size(theta_grid)); for ii = 1:numel(theta_grid) a = exp(-1j*2*pi*d_lambda*(0:M-1)'*sin(theta_grid(ii)*pi/180)); Pm(ii) = 1/real(a'*(Un*Un')*a); % 取实部避免微小虚数杂散 end Pm = Pm/max(Pm); % 峰值提取与显示 [pks, locs] = findpeaks(10*log10(Pm), theta_grid, ... 'MinPeakHeight', -10, 'Annotate','extents', 'MinPeakDistance', 1); disp(locs); plot(theta_grid, 10*log10(Pm));

代码的逻辑拆开看：前三步在构造数据，其中K必须提前知道，这是子空间类算法的前提。第四步的谱公式是1/a^H Un Un^H a，Un的列数等于M-K，所以Un*Un'把特征向量合并成一个投影矩阵。最后一步findpeaks返回峰值角度，MinPeakDistance设为1度能防止同一个谱峰被切出多个局部极大值。这里用了0.1度的扫描步长，实际误差还包括网格量化误差，峰值定位精度通常优于步长的一半，但不会好于克拉美罗界（CRB）。

2.3 信号源个数估计与谱搜索代价

上面代码直接写了K=2，真实系统里信号个数不会提前知道。常见的做法是特征值拉平检验或信息论准则，AIC和MDL都用协方差矩阵的特征值构造代价函数，在M-K比较大的时候，MDL的过估惩罚更重，因此工程上更常用MDL。本质是看特征值序列有没有明显台阶：如果小特征值聚类在同一水平，台阶前的个数就是K。谱搜索前把这个数算错，噪声子空间就会丢列或者混入信号子空间，谱峰会直接消失或多出假峰。

谱搜索的计算量也值得注意。90度范围、步长0.1度一共901个网格点，每个点算一次导向矢量投影，MATLAB里是微秒级，但放进实时系统就变成负担。减少扫描网格到0.5度，再在峰值附近做抛物线插值，网格点降到原来的五分之一，峰值损失很小。正是这个搜索开销，让ESPRIT有了存在的理由：它完全不需要扫描。

3. ESPRIT算法：旋转不变性省掉一次谱搜索

3.1 为什么需要ESPRIT：MUSIC的搜索开销在哪里

MUSIC的问题有两处，一是每个网格点都要计算一次导向矢量投影，二是在低信噪比下谱峰可能被噪声底抬高，导致findpeaks提取的峰值角度抖动。ESPRIT走的是另一条路线：利用均匀线阵相邻阵元的固定相位差，把角度估计从“搜索谱峰”变成“特征值分解”。整个流程只有一次协方差矩阵构造、一次特征分解、一次小矩阵求逆，没有网格循环。所以对实时性敏感的项目，比如车载雷达目标角度估计，ESPRIT的优先级通常比MUSIC高。同一组数据对比，900点的MUSIC扫描大约需要几毫秒，ESPRIT在数百微秒量级，相差约一个数量级。

3.2 子阵划分与旋转方程推导

把M元阵列按“去掉首尾各一个”的方式拆成两个子阵：子阵1取第1到M-1个阵元，子阵2取第2到M个阵元。同一个信号到达两个子阵之间存在固定相位差exp(-j2πd sinθ/λ)，因为两子阵完全平行、间距正好是d。在无噪声条件下，两次快拍矩阵满足X2 = X1 Φ，Φ是K×K对角阵，对角元正好是各源的旋转因子。实际数据有噪声，所以要先投影到信号子空间再做处理。

先对全阵列协方差做特征分解，取前K个特征向量组成Es。由于Es的列是信号子空间的一组基，两个子阵的信号子空间分别是Es的第一部分和第二部分，记E1 = Es(1:M-1,:)，E2 = Es(2:M,:)。理想情况下二者满足E2 = E1 Ψ，Ψ是K×K矩阵，其特征值就是旋转因子。这就是LS-ESPRIT的核心式子，Ψ的最小二乘解写作E1 \ E2。得到Ψ后对它做特征分解，得到K个复特征值φ，每个φ的相位角包含角度信息：

θ = -asin( angle(φ) / (2πd/λ) )

负号的来历和构造导向矢量时的相位符号约定有关，实际使用时以你写A的符号为准。有一个细节值得展开：LS-ESPRIT假设E1没有误差只把误差放在E2，这是普通最小二乘的标准假设；如果M较大但快拍数少，E1和E2都有扰动，用TLS-ESPRIT会更稳。TLS的做法是对扩展矩阵[E1 E2]做奇异值分解，再从右奇异向量的后K列里提取旋转算子，MATLAB实现多几行，但低信噪比下方差能压下来一点。

另一个常见预处理是协方差矩阵的Toeplitz化。有限快拍下样本协方差并不严格满足Toeplitz结构，把Rxx的对角线按同一方向取平均再重排成实对称Toeplitz矩阵，等效于强制恢复平稳性假设，能改善小快拍下的ESPRIT估计稳定性。这个技巧在MUSIC里也可以直接用，代价是若信号存在多径时平滑作用会改变秩结构，需要配合空间平滑一起判断。

3.3 ESPRIT的MATLAB代码与角度换算

% esprit_doa_demo.m 与MUSIC共用同一组数据模型 M = 8; d_lambda = 0.5; K = 2; N = 500; SNR_dB = 10; theta_true = [-10 20]; angs = theta_true*pi/180; A = exp(-1j*2*pi*d_lambda*(0:M-1)'*sin(angs)); S = (randn(K,N) + 1j*randn(K,N))/sqrt(2); sigma_n = sqrt(10^(-SNR_dB/10)); X = A*S + sigma_n*(randn(M,N) + 1j*randn(M,N))/sqrt(2); Rxx = (X*X')/N; [V,D] = eig(Rxx); [ev, idx] = sort(diag(D),'descend'); Es = V(:, idx(1:K)); % 信号子空间 E1 = Es(1:M-1, :); E2 = Es(2:M, :); Psi = E1 \ E2; % 最小二乘旋转算子 [Vec, rPhi] = eig(Psi); phi = diag(rPhi); theta_est = -asin( angle(phi) / (2*pi*d_lambda) ) * 180/pi; theta_est = sort(theta_est.'); disp('ESPRIT估计角度:'); disp(theta_est);

几个值得说清的点。第一，Es到E1、E2的划分必须严格对齐，如果你先对Es做了某种线性变换再去划分子阵，旋转不变性就不成立了。第二，Psi = E1 \ E2本质是解超定最小二乘，因为E1是(M-1)×K矩阵，行数大于列数，用反斜杠前最好加一句if rank(E1)<K, return; end检查秩。第三，angle(phi)得到的是(-π,π]之间的相位，当d/λ大于0.5或者真实角度接近±90度时，相位可能越过周期边界，出现“角度模糊”。所以ESPRIT和MUSIC一样默认d≤λ/2。

ESPRIT也有自己的短板：它对相干源更敏感。原因在于相干信号使得信号子空间维度不足，前K个特征向量不再能正确张成旋转不变子空间。另一个坑是特征值的输出顺序不固定，好在DOA场景下每个φ自带角度信息，按角度排序即可；如果你之后要做二维DOA，还需要专门的配对步骤。

4. 压缩感知DOA：稀疏重构与网格失配

4.1 稀疏基与观测模型：为什么角度域能写成稀疏问题

子空间算法依赖信号源个数K的特征值分布，而压缩感知DOA换了一个视角：把整个角度范围离散成G个网格点，每个网格对应一个导向矢量，把所有导向矢量排成一个M×G的过完备字典A_dict。真实信号只占据其中K个网格，所以字典域上的信号系数向量s是稀疏的。一次快拍模型变成x = A_dict s + n，DOA估计变成从恢复出的s里找非零位置。这个模型的好处是不需要预先知道K，稀疏重构会自动决定支撑集个数。

但字典的过完备性是一把双刃剑。网格越多，原子之间越接近线性相关，矩阵的相干性上升，精确恢复条件更难满足。理论上要求A_dict满足限制等距性质（RIP），实际阵列流形矩阵很难严格满足，所以CS-DOA在低信噪比下的表现和噪声上界的取值关系很大。这是用CS做DOA时最容易踩的坑：不是算法本身不好，而是参数场景没对上。

4.2 用CVX求解L1最小化：代码与参数说明

% cs_doa_demo.m 单快拍CS-DOA重构 M = 8; d_lambda = 0.5; N = 200; K = 2; SNR_dB = 15; theta_true = [-10 20]; theta_grid = -90:0.5:90; G = numel(theta_grid); % 过完备字典 A_dict = exp(-1j*2*pi*d_lambda*(0:M-1)'*sin(theta_grid*pi/180)); % 生成数据 angs = theta_true*pi/180; A_true = exp(-1j*2*pi*d_lambda*(0:M-1)'*sin(angs)); S = (randn(K,N) + 1j*randn(K,N))/sqrt(2); sigma_n = sqrt(10^(-SNR_dB/10)); X = A_true*S + sigma_n*(randn(M,N) + 1j*randn(M,N))/sqrt(2); % 取第100个快拍做单快拍重构 x = X(:, 100); epsilon = sqrt(M)*sigma_n*1.2; % 噪声上界，带安全系数 cvx_begin quiet variable s_hat(G) complex; minimize( norm(s_hat, 1) ); subject to norm(x - A_dict*s_hat, 2) <= epsilon; cvx_end mag = abs(s_hat); [pks, locs] = findpeaks(mag, theta_grid, 'MinPeakHeight', 0.3*max(mag)); disp('CS估计角度: '); disp(locs);

CVX代码的逻辑分三段：构造字典和数据是前置条件；cvx_begin块里，目标函数是s_hat的L1范数，约束是残差的L2范数小于等于epsilon，这是标准的基追踪去噪BPDN形式；恢复结果画成幅度谱后，非零峰值所在网格就是DOA。这里的epsilon是唯一的自由参数，取sqrt(M)*sigma_n是噪声理论上界，乘以1.2的安全系数是为了防止噪声实际超过理论值。epsilon放得太大，残差约束过于宽松，s_hat会被压成全零或只有一个主峰；epsilon太小，则会把噪声也解释成信号，出现密集的伪峰。

多快拍场景下，逐快拍做L1重构再平均峰值位置，效果往往不如联合稀疏重构。联合稀疏模型把s_hat扩展成G×N矩阵，目标函数用L2,1范数（每一行的L2范数之和），多个快拍共享同一个支撑集。CVX里用sum(norms(s_hat,2,2))可以表达，但求解规模变大，N个快拍要存出一个G×N复数矩阵，内存和耗时成倍增加。因此实际项目中，我一般先试单快拍，确定支撑集稳定了，再决定要不要上联合稀疏。

4.3 网格失配与离网处理：最难绕开的坑

经典网格化CS-DOA假设真实角度正好落在网格上。用0.5度网格时，如果真实角度是20.2度，恢复结果会在20度和20.5度各分一部分能量，选主峰要么整体偏移，要么在两次仿真之间来回跳。处理网格失配有三种常见做法。

第一种是把网格加密到0.1度，代价是字典相关性变高，相邻原子几乎相同，L1解的支撑会在两个原子间抖动。第二种是两阶段法，先用0.5度网格粗定位，然后在粗估计附近取±3度的小区间，再以0.05度细分做第二次重构，既控制字典规模又降低量化误差。第三种是离网模型，把角度偏差直接放进优化变量联合估计，比如原子范数最小化，实现复杂度明显上升。我一般先观察粗网格恢复出的支撑是否落在网格点上，如果峰呈现“跨格分布”而不是集中单点，再决定是否上两阶段。不要一上来就把网格做到0.05度，字典规模和内存占用很快会让你后悔。

5. 常见问题与参数避坑：阵列、相干源与快拍数

5.1 阵元间距超过半波长导致角度模糊

现象：把d_lambda设为1.0，MUSIC谱在正确角度附近出现多个等幅峰，ESPRIT算出的角度也明显超出真实值。

原因：相邻阵元相位差是2πd sinθ/λ，当d/λ>0.5时，sinθ在[-1,1]内变化会让相位差跨过2π周期，不同角度可能产生相同相位差，这就是空间混叠。

解决：仿真中固定d=0.5λ；如果实际阵列间距受物理布置限制无法缩小，就只做限定扇区的DOA，利用先验信息排除模糊角度。

提示：做二维DOA时，x轴和y轴的阵元间距都要满足半波长条件，只约束一个方向没有意义。

5.2 相干源导致秩亏、谱峰消失

现象：两个信号来自同一发射机，经过多径到达阵列，二者强相干，MUSIC谱只有一个峰，幅度还很低。

原因：相干使得信号协方差矩阵Rs奇异，Rxx中A Rs A^H的秩从K降到1，特征分解只能分出一个大特征值，Un里混入了第二个信号的分量，正交性被破坏。

解决：做前向空间平滑，把M元阵划分成L个重叠子阵，每个子阵长度M_sub = M-L+1，对各子阵协方差求平均。平滑次数L要满足L≥K，子阵长度M_sub≥K+1，平滑后等效阵元数下降但秩恢复。代码可以插在MUSIC特征分解前：

% 前向空间平滑示意：在MUSIC特征分解之前插入 L = 2; % 平滑次数（子阵个数） M_sub = M - L + 1; % 子阵长度 R_sm = zeros(M_sub, M_sub); for ii = 1:L X_sub = X(ii:ii+M_sub-1, :); R_sm = R_sm + (X_sub*X_sub')/N; end Rxx = R_sm / L; % 平滑后的协方差

平滑后Rxx变成M_sub×M_sub，后续特征分解和谱搜索的维数全部要换成M_sub。前向平滑会牺牲阵列孔径，L越大，有效阵元数越少，角度分辨率越低，所以L只要取刚好满足秩恢复的值，不要贪多。前后向平滑在独立源情况下能把有效快拍翻倍，但在相干场景里要配合前向平滑一起用，单独使用对秩亏的帮助有限。我自己的血泪经验是：先判断信号是否相干，再决定这行代码加不加，否则可能把本来正常的角度谱弄变形。

5.3 快拍数太少导致谱峰分裂

现象：N从500降到20，MUSIC谱主峰变平、旁边多出小峰，峰值角度偏移3度以上。

原因：样本协方差和真实协方差偏差变大，小特征值不再集中，噪声子空间被扰动，谱底不再平坦。

解决：优先增加快拍数到100以上；如果不允许，采用前后向平均把有效快拍翻倍；还不行就在Rxx对角线上加对角加载量βI，β取0.1到0.5倍小特征值均值，压制噪声子空间抖动。注意对角加载会抬高谱底，峰高相对下降，所以加载量不应超过最小信号特征值。调试时先画一次谱看谱底高度，再决定加载量，不要盲目取固定值。

5.4 三种算法的适用边界对比

这张表是选型参考，不是死规矩。MUSIC在低信噪比下的表现其实和ESPRIT接近，二者都受快拍数限制，差别在于MUSIC的谱峰形态更直观，便于人工观察；CS在低快拍数下有理论优势，但前提是网格失配可控。如果你的场景是阵元数只有4、信噪比低于5dB，先做平滑和累积快拍，再谈算法选择，这一条经验比任何算法改进都管用。

6. 三种算法同台验证：仿真中值得盯住的指标

6.1 用同一组数据算RMSE，评价调参是否有效

把三种算法放进同一个蒙特卡洛循环，跑200次，每次重新生成噪声但保持真实角度不变，计算估计角度的均方根误差：

% compare_doa.m 三种算法RMSE对比骨架 nMC = 200; theta_true = [-10, 20]; err_music = zeros(nMC,2); err_esprit = zeros(nMC,2); err_cs = zeros(nMC,2); for mc = 1:nMC X = genData(M, theta_true, N, SNR_dB); % 同前文的数据生成函数 est_mu = sort(music_doa(X, M, K)); est_es = sort(esprit_doa(X, M, K)); est_cs = sort(cs_doa(X(:,1), A_dict, sigma_n)); err_music(mc,:) = abs(est_mu - theta_true); err_esprit(mc,:) = abs(est_es - theta_true); err_cs(mc,:) = abs(est_cs - theta_true); end rmse_music = sqrt(mean(err_music.^2)); rmse_esprit = sqrt(mean(err_esprit.^2)); rmse_cs = sqrt(mean(err_cs.^2));

注意排序要对齐，因为ESPRIT和CS返回的角度顺序不固定。RMSE之外还要统计“失败率”，比如定义偏航大于3度为一次失败，这个指标比平均误差更能反映工程可靠性。调参顺序我固定为：先确认K，再调d_lambda和去相干，最后优化网格步长或epsilon。K错了，后面全是白做。

6.2 网格步长、快拍数、加载系数：哪些参数值得最后打磨

仿真阶段最值得打磨的三个参数依次是信号源数K的估计（用MDL做参考）、平滑子阵数L（和阵元数成反比）、CS的epsilon（先按理论值，再用小范围扫描）。这三个变量直接决定最终是收敛到CRB还是在3度附近游走。经历过一次ESPRIT把20度估计成17度的事故后，我养成了一个习惯：每次改完参数，先把估计结果和真实角度的差值打印成数组扫一眼，不要只看RMSE的平均数，因为平均数的光滑会掩盖个别仿真的大偏差。从那以后我每次跑DOA仿真都强制走一遍“先盯发散点、再调参数”的流程。这套源码包里的MUSIC.m、esprit.m、CS.m可以作为这个流程的起点，把阵元数、快拍数和角度改成你自己的场景，三种算法的差异很快就能摸透。希望帮到你。

本文还有配套的精品资源，点击获取