简介：针对阵列信号处理中的来波方向估计问题，这份MATLAB源码包为学习者和工程师提供了三条经典技术路线的可执行参考。资源共3个文件，均为.m脚本，压缩包仅4KB，分别对应多重信号分类算法、旋转不变参数估计算法与压缩感知方法。包内代码紧凑，每行均有清晰结构，便于结合协方差矩阵分解、旋转不变性与稀疏重构等原理逐步理解，压缩感知脚本还演示了低采样率下利用稀疏性恢复角度。用户可以修改阵元数、快拍数或信噪比进行对比实验，观察分辨率、稳健性及计算量的差别。已有1474人学习下载，适合用于课程设计、科研入门或工程预研。运行这些脚本，不仅能巩固理论认识，还能掌握在MATLAB中搭建来波方向估计实验的基本流程，为进一步开发改进算法提供起点。

1. DOA算法是什么：它不是测距，是测方向

DOA（Direction of Arrival）算法在阵列信号处理里的地位，和FFT在频谱分析里的地位差不多。你有一排天线或者麦克风，收到的不只是一路信号，而是多条带有时间延迟的副本——DOA算法要回答的问题只有一个：信号从哪个角度打过来的。它不关心目标离你多远，不关心信号强到什么程度，只关心方向。这个能力在雷达、声呐、被动侦测、麦克风阵列语音增强里都是刚需。

很多人第一次接触DOA时会被“算法”两个字吓住，觉得要先把什么子空间、特征分解啃明白才能动代码。实际上用MATLAB跑通一个最基础的MUSIC算法，核心代码不超过30行。真正的难点不在代码本身，而在三件事：信号模型怎么建、协方差矩阵怎么算对、参数怎么调得让谱峰干净。这篇文章就按这个顺序，把DOA算法的原理和MATLAB实现串成一条能直接复现的路径。

2. 均匀线阵与DOA信号模型：先有阵列，才有“方向”这一说

2.1 为什么阵元间距要卡死在半波长

DOA估计的前提是“同一个信号被多个阵元接收时，阵元之间存在稳定的相位差”。这个相位差来自信号到达不同阵元的时间差。要想让时间差和入射角度一一对应而不产生模糊，阵元间距d必须满足d ≤ λ/2，λ是信号波长。工程上最常见的选择就是d = λ/2，既满足不模糊条件，又能最大化阵列孔径，从而获得最好的角度分辨率。

注意这个半波长约束是硬性的。如果你设d = 0.8λ，那么某些入射角会出现“栅瓣”，你会在谱里看到一个假峰，真目标旁边多出一个鬼影。我见过不少人在仿真里为了省事随便定d，结果谱峰乱成一团，还以为是算法写错了——其实从头到尾就是阵元间距不符合采样定理。阵列采样在空间上同样遵循Nyquist定理，这算是最容易被忽略的一层物理含义。

2.2 窄带远场模型的数学表达：导向矢量

经典的DOA分析都建立在窄带远场模型上，也就是：信号频率单一，且波前面在阵列孔径范围内可以看成平面波。这个模型下，第m个阵元相对第1个阵元的相位延迟写成：

[ a(\theta) = \exp(-j \cdot 2\pi \cdot d \cdot (m-1) \cdot \sin(\theta)/\lambda) ]

把M个阵元的响应叠成一个列向量，就是导向矢量 (a(\theta))。整个阵列接收到K个来自不同方向 (\theta_1, \theta_2, ..., \theta_K) 的信号时，单次快拍的阵列输出可以写成：

[ \mathbf{x} = \mathbf{A}(\theta) \cdot \mathbf{s} + \mathbf{n} ]

其中 (\mathbf{A}) 是M×K的导向矢量矩阵，(\mathbf{s}) 是K×1的信号矢量，(\mathbf{n}) 是噪声。这个公式是后面所有DOA算法的出发点。MUSIC和ESPRIT都是在这条基本等式上做文章，区别在于谁在利用噪声子空间，谁在利用信号子空间的旋转不变性。

在MATLAB里构建导向矢量时，最容易出错的地方是exp里面的符号。不同教材定义阵列方向的方式不同，有的从左边入射为正，有的从右边为正，但只要你保持一致，最后谱峰位置就不会受符号影响——除非你在构建导向矢量和在谱峰搜索时用了两套不同的符号约定，那就一定会翻车。

2.3 协方差矩阵其实是“能量布局图”

阵列输出的协方差矩阵定义为：

[ \mathbf{R} = E[\mathbf{x} \mathbf{x}^H] ]

它是阵列信号处理的“黑匣子里最核心的黑匣子”。理论上，(\mathbf{R}) 包含信号方向的全部信息；实践上，我们只有N次快拍，只能用样本协方差矩阵 (\hat{\mathbf{R}} = (1/N) \sum \mathbf{x}_n \mathbf{x}_n^H) 去估计它。这个估计的质量直接决定DOA算法的角度估计精度。

一个关键性质是：(\mathbf{R}) 的特征分解会自然分出两个子空间。大特征值对应的特征向量张成信号子空间，小特征值对应的特征向量张成噪声子空间。信号子空间和噪声子空间在理想条件下是完全正交的——这正是MUSIC算法能分辨任意角度间隔的信号的根本原因。你在MATLAB里用eig(R)也好，用svd(R)也好，都是在拿这套正交性去做文章。

写代码时要注意R = (X * X') / N和R = (X * X') / (N-1)的区别：仿真里差一个常数因子无所谓，因为后面谱峰搜索都是比值运算；但如果你同时用R做信噪比估计，就应该用/N，因为( \mathbf{R} )的定义里就是期望，不是统计学的样本方差。这个细节虽然不影响谱峰位置，但会影响你对噪声功率的读数。

3. 三种DOA算法的原理与MATLAB实现：从功率谱到子空间

3.1 Capon波束形成：有约束的最小方差估计

Capon算法，也叫MVDR（Minimum Variance Distortionless Response），思路很直观：我想听某个方向 (\theta) 的信号，那就设计一组权矢量 (\mathbf{w})，让输出功率 (P = \mathbf{w}^H \mathbf{R} \mathbf{w}) 最小，同时约束 (\mathbf{w}^H \mathbf{a}(\theta) = 1)，保证想听的方向信号不衰减。

这其实是一个带等式约束的二次优化问题，用拉格朗日乘子法可以直接解出：

[ \mathbf{w}_{capon} = \frac{\mathbf{R}^{-1} \mathbf{a}(\theta)}{\mathbf{a}^H(\theta) \mathbf{R}^{-1} \mathbf{a}(\theta)} ]

对应的输出功率谱为：

[ P_{capon}(\theta) = \frac{1}{\mathbf{a}^H(\theta) \mathbf{R}^{-1} \mathbf{a}(\theta)} ]

Capon的意义在于它自动在干扰方向形成了零陷，比传统波束形成的角度分辨率高。但是它对阵列误差和协方差矩阵估计误差比较敏感，快拍数少时谱峰容易变钝。

3.2 MUSIC算法：信号子空间与噪声子空间的正交性来寻峰

MUSIC（Multiple Signal Classification）是DOA领域最出名的子空间类方法。它的核心洞察是：导向矢量 (\mathbf{a}(\theta)) 如果正好对应某个真实信号方向，那么它与噪声子空间正交，即 (\mathbf{a}^H(\theta) \mathbf{U}_n = \mathbf{0})。

在实现时，我们没法验证严格为0，只能找让 (\mathbf{a}^H(\theta) \mathbf{U}_n) 的范数最小的角度。所以MUSIC谱定义为：

[ P_{music}(\theta) = \frac{1}{\mathbf{a}^H(\theta) \mathbf{U}_n \mathbf{U}_n^H \mathbf{a}(\theta)} ]

分母接近0时，谱值趋向无穷大，谱峰就出来了。由于MUSIC用上了噪声子空间的正交性质，它的分辨率比Capon高得多，同样是两个间隔5度的信号，MUSIC能分开而Capon很可能只有一个馒头状谱峰。

3.3 MATLAB仿真代码：一个脚本同时跑通Capon和MUSIC

%% DOA仿真：均匀线阵，两个信号源，对比Capon与MUSIC谱 clear; clc; close all; %% 参数设置 M = 8; % 阵元数 d_lambda = 0.5; % 阵元间距 / 波长 N = 200; % 快拍数 K = 2; % 信源数 theta = [-10 25]; % 真实入射角（度） SNR_dB = 10; % 信噪比 %% 生成窄带远场信号 lambda = 1; % 归一化波长 d = d_lambda * lambda; % 实际阵元间距 angles = theta * pi/180; % 角度转弧度 % 构建导向矢量矩阵 A：每列对应一个来波方向 A = exp(-1j * 2 * pi * d * (0:M-1)' * sin(angles) / lambda); % 注意共轭转置 % 生成信号与噪声 S = (randn(K, N) + 1j * randn(K, N)) / sqrt(2); % K路零均值复信号 X = A * S; % 阵列输出 sigma2 = 10^(-SNR_dB/10); % 噪声功率 X = X + sqrt(sigma2/2) * (randn(M, N) + 1j * randn(M, N)); %% 样本协方差矩阵 R = (X * X') / N; % 注意是共轭转置，不是转置 %% 角度搜索 theta_scan = linspace(-90, 90, 181); % 从-90到90度，每度1个点 P_capon = zeros(size(theta_scan)); P_music = zeros(size(theta_scan)); % 对协方差矩阵做特征分解供MUSIC使用 [U, D] = eig(R); [eigval, idx] = sort(diag(D), 'descend'); U_s = U(:, idx(1:K)); % 特征值最大的K个向量 => 信号子空间 U_n = U(:, idx(K+1:end)); % 其余M-K个向量 => 噪声子空间 for i = 1:length(theta_scan) a = exp(-1j * 2 * pi * d * (0:M-1)' * sin(theta_scan(i)*pi/180) / lambda); P_capon(i) = 1 / abs(a' * (R \ a)); % Capon谱 P_music(i) = 1 / abs(a' * (U_n * U_n') * a); % MUSIC谱 end %% 归一化画图 P_capon_dB = 10*log10(abs(P_capon) / max(abs(P_capon))); P_music_dB = 10*log10(abs(P_music) / max(abs(P_music))); plot(theta_scan, P_capon_dB, 'b', 'LineWidth', 1.5); hold on; plot(theta_scan, P_music_dB, 'r', 'LineWidth', 1.5); xline(theta(1), '--k'); xline(theta(2), '--k'); xlabel('入射角 (deg)'); ylabel('归一化谱 (dB)'); legend('Capon','MUSIC','真实方向'); grid on;

这段代码的核心逻辑是三步：先生成阵列输出数据，再算协方差矩阵，最后在角度网格上遍历，用各自的谱公式算功率。MUSIC部分注意特征分解后的排序：eig函数返回的特征值不一定是降序排列，必须先排序再切分信号子空间和噪声子空间，否则你切出来的子空间是乱的，谱峰位置会完全错掉。

参数说明里最值得留意的是SNR_dB：这里的信号功率设为1，噪声功率按公式10^(-SNR_dB/10)换算。如果你想看高信噪比下的谱峰锐度，把SNR拉到20；想看低信噪比下限，SNR=0也能出峰，但MUSIC会出现明显的旁瓣起伏。

3.4 ESPRIT算法：免栅格搜索的旋转不变法

MUSIC要在每个角度上算一次谱，计算量大；ESPRIT则是绕开了搜索的另一种思路。它把阵列分成两个完全一样的子阵，利用两个子阵之间的旋转不变性来直接以闭式解求角度，不需要遍历角度网格。

ESPRIT的求解过程是：取前M-1行作为子阵1，后M-1行作为子阵2。对信号子空间做两个子阵的映射 (\mathbf{U}{s1}, \mathbf{U}{s2})，然后求解 (\mathbf{\Psi} = \mathbf{U}{s1}^+ \mathbf{U}{s2})，它的特征值包含所有到达角度信息。入射方向为：

[ \theta_k = \arcsin\left( \frac{\lambda}{2\pi d} \cdot \arg(\lambda_k) \right) ]

其中 (\lambda_k) 是 (\mathbf{\Psi}) 的第k个特征值。MATLAB实现时把特征值取相位再除以(2\pi d/\lambda)，用asin转回角度即可。由于不需要网格搜索，ESPRIT的运算速度远快于MUSIC，特别适合实时处理；代价是它对子阵划分方式有要求，且两个子阵之间必须完全一致，工程上硬件不匹配时会引入系统性偏差。

4. DOA算法的三个必调参数：阵元数、快拍数、信噪比

4.1 阵元数M限制最大可分辨信源数

均匀线阵条件下，M阵元的阵列最多能分辨M-1个独立信号源。这是MUSIC算法的硬限制——因为你需要至少1个噪声子空间向量，如果K = M，信号子空间占了全部维度，噪声子空间就变成零维了。在这个代码框架里，要注意K值必须小于M，否则后面按U(:, K+1:end)切片会出错。

有些教材里说MUSIC能分辨M-1个源，也有说到M个的，区别在于你是否把全零特征值对应的向量也算进噪声子空间。实际工程中，通常建议K不超过M-2，毕竟协方差矩阵估计有误差，留点冗余子空间能让谱峰更稳。

4.2 快拍数N带来的方差变化

快拍数决定了协方差矩阵 (\hat{\mathbf{R}}) 对真实R的逼近程度。快拍数太少，(\hat{\mathbf{R}}) 的噪声特征值发散得厉害，MUSIC谱峰旁边会出现许多毛刺，甚至出现伪峰。我一般建议：仿真验证至少取N = 100，工程实时处理至少N = 20～50，但N再低就不要再往下压了。

如果你想做一个实验观察快拍数的实际影响，可以让N从10变化到500，记录峰位置与真实角度的偏差。你会发现N < 30时方差飙升，N > 200以后曲线的斜率明显变缓——这说明快拍数不是越多越好，存在一个边际收益递减的拐点。做法是把N写成一个循环，每次重新生成信号并算一次DOA估计，统计标准差曲线。

4.3 信噪比与角度分辨率的经验关系

角度分辨率并不只由阵元数决定，信噪比同样在起作用。MUSIC理论优势在于它能把两个紧邻的信号分开，但前提是它们的空间特征在噪声背景下仍然可辨识。经验法则是：两个信号的角度间隔至少大于(\sim (2\pi/M))乘以由信噪比决定的一个系数。画成曲线的话，你会看到低信噪比时分辨率急剧恶化，高信噪比时逼近阵列孔径的物理极限。

想在仿真里量化这个现象，固定M = 8，让两个信号角间隔从10度逐步降到1度，在SNR = 5dB和SNR = 20dB下分别跑MUSIC，看哪个间隔下谱峰开始合并。你会发现SNR=20时1度还能勉强分开，SNR=5时3度以内就已经分辨不出来了。

5. DOA仿真避坑指南：五个常把结果带偏的工程细节

5.1 协方差矩阵用错转置：找到谱峰的“玄学”偏差

现象：MUSIC谱峰位置与真实入射角始终差一个固定量，且信号入射角越接近零度越明显。

原因：MATLAB里X'是共轭转置，X.'才是普通转置。如果信号是复数，用X.'算协方差矩阵会丢掉共轭信息，(\hat{\mathbf{R}}) 不再是Hermitian矩阵，特征分解产生虚部误差，谱峰自然偏移。

解决：一律用R = (X * X') / N，写代码时给这行加注释“共轭转置勿写成X.'”。这个坑在实数域仿真里不出现，所以很难通过简化模型来发现，只有上复数信号时才暴露。

5.2 角度搜索步长太粗：谱峰直接“丢”在网格缝隙里

现象：真实角度是-3.7度，搜索网格粒度是1度，MUSIC谱峰出现在-4度和-3度之间的某个位置，但幅值远低于旁边的伪峰。

原因：MUSIC谱的峰很尖锐，尤其是高信噪比时，峰值宽度远小于1度。网格太粗可能根本没有采样到真正的峰值点，最大输出指向了离网格最近的旁瓣。

解决：第一轮用linspace(-90, 90, 181)做粗扫描找出峰的位置，然后在峰附近再做一轮细扫描，linspace(peak-2, peak+2, 501)，步长精确到0.01度量级。这样既控制总计算量，又能拿到精确的峰值位置。

5.3 信号源数K设定错误：噪声子空间被污染

现象：谱图全高，几乎没有明显的峰值；或者峰的数量明显多于真实信号数。

原因：K设大了，把信号子空间的维度撑大，真正的信号特征向量被挤进噪声子空间，正交性被破坏，谱自然乱。K设小了，信号子空间容量不足，剩余信号泄漏到噪声子空间，会出现伪峰。

解决：仿真中K已知，直接赋值没问题；真实数据不知道K，用MDL或AIC准则估计。MDL一般比AIC更保守，过拟合概率更低。简单做法是把[U, D] = eig(R)的特征值画出来，找明显的“拐点”，拐点之后的特征值数量就是M-K，对应噪声子空间维度。

5.4 相干信号（多径）让MUSIC失效：预处理不做等于白跑

现象：同一个信号源经过反射形成多条路径，入射角不同但信号完全相同，MUSIC谱只出一个宽峰，第二个方向丢失。

原因：相干信号使协方差矩阵的信号子空间塌缩，其秩降到1而不是信号路径数，子空间分解自然失效。

解决：做空间平滑预处理。把M元阵列划分成L个重叠子阵，对每个子阵分别算协方差矩阵，然后平均，恢复协方差矩阵的秩。代价是有效阵元数从M降到M-L+1，MUSIC可分辨的信号数也因此减少。这是多径环境下必须付出的代价。

5.5 特征分解后的特征向量排序：看似不起眼，错则全盘皆错

现象：MUSIC代码逻辑看起来完全正确，但谱峰始终出现在错误角度或完全没有峰。

原因：MATLAB的eig返回的特征值和特征向量，顺序是升序的，这是常见诱发点。但你没排序就取后K列当信号子空间，等于把最小特征值对应的特征向量当成了信号——子空间完全反了。

解决：加一行[~, idx] = sort(diag(D), 'descend')做显式降序排序。这个坑在经典教科书代码里常被忽略，因为很多教材用伪代码描述，而伪代码里不会告诉你eig的输出顺序。踩过一回之后，我写任何特征分解代码都会把这行排序放在紧挨着solve上面的位置。

6. 将DOA算法送进工程处理：相干源、宽带信号与未知信源数的三个常用扩展技巧

6.1 相干信号的解算：前向空间平滑怎么嵌进现有MUSIC流程

当信号来自多径时，按常规方式算协方差矩阵不可用，需要先做前向空间平滑。把M个阵元划分成L个长度为M-L+1的重叠子阵，第l个子阵数据是X的第l行到第l+M-L行，对每个子阵算协方差矩阵再求平均：

%% 前向空间平滑预处理 L = 4; % 子阵数，可根据M和信号数调整 sub_len = M - L + 1; % 每个子阵的长度 R_fb = zeros(sub_len, sub_len); for l = 1:L X_sub = X(l:l+sub_len-1, :); R_fb = R_fb + (X_sub * X_sub') / N; end R_fb = R_fb / L;

这段代码的关键是子阵数L的选取。L越大，去相干能力越强，但有效阵元数sub_len越小，能分辨的信号数也跟着下降。经验值是让sub_len至少大于两倍信号数，不然MUSIC的子空间维度不够用。做完空间平滑之后，后面走的MUSIC流程和普通版本一模一样，只需要把谱估计里的(\mathbf{R})换成本文里的R_fb。

工程上，相干源场景里还可以做前后向空间平滑，即把前向平滑矩阵和它的共轭翻转矩阵也做一次平均，额外获得约2倍的平滑效果。代价是计算量翻倍，但对静态场景来说几乎可以忽略。

6.2 宽带信号与窄带模型的对接：频域子带法是底线

DOA模型建立在窄带上，实际声呐和语音信号往往是宽带的。常见做法是把宽带信号分段FFT后，在每个频点单独跑MUSIC，然后将多个频点的空间谱做平均。具体步骤是：数据分帧、FFT到频域、对每个频点构建协方差矩阵、用同一个导向矢量公式在该频点计算MUSIC谱、最后把所有频点的谱加权平均。

这个过程中需要特别注意的是，导向矢量里的波长是当前频点对应的波长，不同频点下的阵元间距与波长的比值d/λ并不一样——你用300Hz的波长构建的导向矢量，不能直接套用900Hz频点的数据。正确写法是对每个频点单独计算d/λ，再带入a = exp(-1j * 2 * pi * d_lambda * (0:M-1)' * sin(theta)*pi/180)的公式里。

频带平均的好处是抗噪，坏处是分辨率被低频段拉低。如果你不追求极致精度，只取3～5个信噪比较高的频点做简单平均即可，不必全频段覆盖。

6.3 信号源数未知：MDL准则自动推断K

真实场景里你事先不知道有几个来波方向，而MUSIC要求输入K。我习惯用MDL（Minimum Description Length）准则做自动推断，代码很短：

%% MDL准则估计信号源数 N_MDL = size(X, 2); % 快拍数 [~, D_mdl] = eig(R); eigvals = diag(D_mdl); [eigvals, ~] = sort(eigvals, 'descend'); mdl_val = zeros(1, M-1); for k = 0:M-2 Lk = M - k; % 当前假设的自由度 avg_log = mean(log(eigvals(k+1:end))); mdl_val(k+1) = -N_MDL * sum(log(eigvals(k+1:end)) - avg_log) ... + 0.5 * k * (2*M - k) * log(N_MDL); end [~, K_est] = min(mdl_val); K_est = K_est - 1; % 下标从1开始对应k从0开始

MDL的核心思想是对模型复杂度和拟合优度做权衡。信源数增加时，拟合优度改善但模型惩罚项变大，MDL曲线的极小值点就是最优点。这段代码在快拍数大于100时表现稳定；快拍数太少时MDL容易低估信号数，这时我宁可在MDL估计值上加1，也不愿意漏掉一个信号。

我自己的习惯是，每做一个DOA项目，都会先写一个类似第3章的通用仿真脚本，把M、N、SNR、角度、K都参数化，然后才在它的基础上加平滑、MDL、宽带扩展这些模块。这样每次拿到新数据，先跑一遍基线流程，确认信号模型匹配，再决定要不要加复杂度。这次分享的这些方法我自己踩过不少坑，特别是特征值排序和相干信号预处理那两个，写在这里希望能帮到你。

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