RELAX算法强弱信号DOA估计的MATLAB实现与参数详解

简介：RELAX算法在雷达DOA估计中的应用实践资源，面向信号处理、阵列测向与雷达目标识别方向的研究者和工程师，重点解决强信号与弱信号并存时弱信号定位困难的问题。该算法通过迭代松弛过程逐步提高弱信号估计精度，相比传统FFT方法在低信噪比下具有更好的分辨能力，适用于超分辨测向、一维信号重构、谐波检测以及基于HRRP的雷达目标识别等场景。资源共17个文件，总大小195KB，包含6个MATLAB源文件、5个文本说明、4个RAR压缩包和1个PDF文档；其中M脚本为核心算法与仿真代码，TXT描述实验原理和操作步骤，RAR包内为完整工程或对比数据，PDF供算法理论参考。内容涵盖RELAX基本实现、CAvsRELAX对比分析以及针对强弱信号DOA的实验代码，读者可直接运行查看效果，也可修改信号参数、阵元配置进行二次实验。已有859人学习下载，适合需要掌握迭代类DOA估计方法并希望获得可直接运行MATLAB代码的初学者和进阶用户。

1. RELAX算法在强弱信号DOA估计里的定位：先看清它能补谁的短板

阵列信号处理里做DOA估计，最头疼的不是两个等强信号靠得太近，而是一个强信号旁边趴着一个弱信号。MUSIC、ESPRIT这类子空间方法在这种场景下常常翻车——强信号的主导特征向量把噪声子空间带偏，弱信号对应的谱峰直接沉到噪声底下面，你怎么调快拍数和窗函数都救不回来。RELAX算法走的是另一条路：它不构造协方差矩阵，而是把每个信源的幅度和角度当成一组未知参数，用“估计最强分量→从数据里减掉→再估计下一个”的迭代方式逐个抠信号。这个思路对强弱信号共存、相干信号、低快拍数场景特别实用，而且不需要训练数据，给一组阵列快照就能跑。这篇笔记直接给你一个能复现的MATLAB实现，把搜索网格、迭代次数、收敛门限这些参数怎么设、坑在哪一次说清楚。

2. RELAX算法原理拆解：迭代相消与幅度-角度联合估计

2.1 从信号模型看RELAX为什么不需要协方差矩阵求逆

考虑M元均匀线阵，接收N个快拍的数据矩阵X，尺寸M×N。假设有K个远场窄带信号从角度θ1…θK入射，阵列流型矩阵A = [a(θ1), …, a(θK)]，其中方向向量a(θ) = [1, e^{j2πd sinθ/λ}, …, e^{j2π(M-1)d sinθ/λ}]^T。数据模型写成：

X = A S + N

S是K×N的复幅度矩阵，N是噪声。MUSIC家族的方法要从X估计协方差矩阵R = XX^H / N，再做特征分解，把空间分成信号子空间和噪声子空间。这个流程有两个先天毛病：第一，R的特征分解需要快拍数N至少大于阵元数M，否则协方差估计秩亏；第二，如果两个信号完全相干（比如多径），R的信号子空间会退化，MUSIC的谱峰直接消失。ESPRIT稍好但也离不开协方差矩阵的构造。

RELAX完全不碰协方差矩阵。它的代价函数是非线性最小二乘形式：

C(θ, s) = ||X - Σ_{k=1}^{K} a(θ_k) s_k||_F^2

其中s_k是第k个信号的复幅度（1×N的行向量）。RELAX用坐标下降法最小化这个代价函数：固定其他K-1个源的参数，只优化第k个源的θ_k和s_k；算完一个马上更新，再算下一个，多轮迭代直到收敛。因为始终面对的是“原始数据减掉其他源贡献后的残差”，强信号残留对弱信号的遮蔽效应会在迭代中被逐步削弱。这跟subspacenet这类端到端深度学习方法也不一样——深度学习DOA需要大量带标注训练数据和场景适配，RELAX只要阵列标定没问题，换频段换阵型改个方向向量就能跑。

2.2 迭代主循环：估计-相消-再估计的三步闭环

RELAX的核心执行顺序可以拆成两步：初始化，然后进入主循环。初始化阶段最关键，我一般用“顺序相消法”：先在原始数据上扫一遍全角度谱，找到最强峰归为第1个源，估计它的幅度并从数据里减掉；再在残差上找下一个峰，重复K次。这样初始角度不会扎堆落在同一个峰上。

% 初始化：顺序相消，逐个找初始角度 theta_init = zeros(1, K); X_res = X; for k = 1:K P = abs(A_grid' * X_res).^2 / N; % A_grid 是预计算的搜索网格方向矩阵 [~, idx] = max(P); theta_init(k) = grid(idx); a = A_grid(:, idx); s_hat = a' * X_res; % 最小二乘投影得到复幅度 X_res = X_res - a * s_hat; % 从数据中相消 end

这里A_grid = exp(1j * 2π * d_lambda * (0:M-1)' * sin(grid * π / 180))，每一列对应一个搜索角度。A_grid需要按列归一化，使a' * a = 1，这样幅度估计s_hat = a' * X_res就是最小二乘解，不需要额外除以M。为什么初始化要顺序相消而不是直接取前K个谱峰？因为强信号的谱峰旁瓣可能比弱信号的主瓣还高，直接取前K个峰大概率会把旁瓣当信号，后面迭代就很难纠正回来了。顺序相消虽然会让初始化偏保守，但胜在稳定。

主循环的每一轮里，对每个k依次执行：构造去掉第k个源后的残差X_others = X - Σ_{j≠k} a(θ_j) s_j，在残差上重新扫描角度谱找最大值位置更新θ_k，再用新角度估计s_k。所有K个源更新完算一次结束。角度搜索和幅度估计是交替进行的，这保证了代价函数单调不增，也就是残差能量只会下降不会上升。

2.3 收敛判据与信源数K怎么处理

收敛判据我习惯用残差能量比而不是角度变化量：每一轮迭代结束后算hist_norm = ||X - Σ a_k s_k||_F^2 / ||X||_F^2，连续两轮的残差比变化小于tol（比如1e-3），就认为收敛。单纯看角度变化不够稳，因为角度在小数点后变动时残差可能还在明显下降；反之，角度已经不变但幅度还在微调，残差比也在降。所以残差比是更可靠的收敛信号。

% 用残差比判断信源数是否足够 K_max = 5; ratio = zeros(1, K_max); for K = 1:K_max [~, ~, hist_norm] = relax_doa(X, d_lambda, K, opts); ratio(K) = hist_norm(end); end % 看残差比随K的下降曲线，找拐点 d_ratio = abs(diff(ratio)); K_est = find(d_ratio < 0.05 * d_ratio(1), 1);

这段代码的思路很简单：K给多了，多出来的源会去拟合噪声，残差比下降幅度明显变小；K给少了，残差比会处在一个高位下不去。实际使用中，如果已知信号数上限，可以用MDL/AIC先粗定阶，再用这个残差比曲线做二次确认。我通常在仿真里直接把K当已知参数传进去，但在实测数据里一定跑一遍这个拐点诊断——实测环境的信源数往往比你以为的多。

3. 用MATLAB实现RELAX-DOA估计：核心函数与可运行脚本

3.1 relax_doa主函数：参数、循环与输出

下面这个函数是我常用的RELAX实现，核心逻辑是初始化顺序相消加主循环交替投影。输入参数全部收敛到一个opts结构体里，避免函数签名过长。

function [theta_est, s_est, hist_norm] = relax_doa(X, d_lambda, K, opts) % RELAX_DOA 基于RELAX迭代相消的波达方向估计 % 输入: % X - M x N 复数矩阵，M个阵元N个快拍 % d_lambda - 阵元间距与波长之比，通常取0.5 % K - 信源数目 % opts - 结构体，可选字段: % .grid_res 搜索步长(度)，默认0.1 % .max_iter 最大迭代轮数，默认10 % .tol 残差比收敛门限，默认1e-3 % .verbose 是否打印迭代信息，默认false % 输出: % theta_est - 1 x K 角度估计(度) % s_est - 1 x K 复幅度估计 % hist_norm - 1 x iter 每轮残差能量比 if nargin < 4, opts = struct(); end grid_res = 0.1; max_iter = 10; tol = 1e-3; verbose = false; if isfield(opts, 'grid_res'), grid_res = opts.grid_res; end if isfield(opts, 'max_iter'), max_iter = opts.max_iter; end if isfield(opts, 'tol'), tol = opts.tol; end if isfield(opts, 'verbose'), verbose = opts.verbose; end [M, N] = size(X); grid = -90:grid_res:90; A_grid = exp(1j * 2 * pi * d_lambda * (0:M-1)' * sin(grid * pi / 180)); A_grid = A_grid / sqrt(M); % 归一化使 a'*a = 1 % ---------- 初始化：顺序相消 ---------- theta = zeros(1, K); X_res = X; for k = 1:K P = abs(A_grid' * X_res).^2 / N; [~, idx] = max(P); theta(k) = grid(idx); a = A_grid(:, idx); s = a' * X_res; X_res = X_res - a * s; end % ---------- 主循环：交替投影 ---------- s = zeros(K, N); hist_norm = zeros(1, max_iter); for it = 1:max_iter for k = 1:K X_others = X; for j = 1:K if j ~= k idx_j = round((theta(j) + 90) / grid_res) + 1; idx_j = max(1, min(length(grid), idx_j)); a_j = A_grid(:, idx_j); X_others = X_others - a_j * s(j, :); end end P = abs(A_grid' * X_others).^2 / N; [~, idx] = max(P); theta(k) = grid(idx); s(k, :) = A_grid(:, idx)' * X_others; end X_res = X; for j = 1:K idx_j = round((theta(j) + 90) / grid_res) + 1; idx_j = max(1, min(length(grid), idx_j)); a_j = A_grid(:, idx_j); X_res = X_res - a_j * s(j, :); end hist_norm(it) = norm(X_res, 'fro')^2 / norm(X, 'fro')^2; if verbose fprintf('iter %2d, 残差比 %.4e, 角度=[%s]\n', ... it, hist_norm(it), num2str(theta, '%.3f ')); end if it > 1 && abs(hist_norm(it) - hist_norm(it-1)) < tol break; end end theta_est = theta; s_est = s(1, :); % 这里只返回行向量示意，多源场景建议返回完整s矩阵 end

几个参数设计的逻辑说明。grid_res默认0.1度，对应-90到90度共1801个网格点，A_grid矩阵是M×1801，内存占用不大。如果阵元数M超过32，建议把A_grid改成在每次搜索时逐列计算，否则矩阵乘法的开销会拖慢速度。max_iter设10通常够用，因为RELAX是坐标下降，前三轮下降最快，后面主要是微调。tol设1e-3是基于残差能量比的相对变化，比固定绝对阈值更适应不同信噪比场景。收敛后返回的hist_norm(end)就是最终的残差能量比，可以用来判断信号是否被提取干净。

3.2 两源强弱信号仿真脚本：从数据生成到画图

拿到函数之后，先用一个典型场景做冒烟测试：10阵元、200快拍、两个信号一个强一个弱，信噪比差20dB。真实角度设在-10度和12度，间隔22度。

% 强-弱信号DOA仿真脚本 rng(2024); M = 10; % 阵元数 N = 200; % 快拍数 d_lambda = 0.5; % 半波长布阵 theta_true = [-10, 12]; % 真实角度 amp = [1, 0.1]; % 强信号幅度1，弱信号幅度0.1 noise_power = 0.01; % 噪声功率，弱信号SNR约等于0dB % 生成阵列数据 A = exp(1j * 2 * pi * d_lambda * (0:M-1)' * sin(theta_true * pi / 180)); S = amp(:) .* (randn(length(theta_true), N) + 1j * randn(length(theta_true), N)) / sqrt(2); X = A * S + sqrt(noise_power / 2) * (randn(M, N) + 1j * randn(M, N)); % 调用RELAX opts = struct('grid_res', 0.1, 'max_iter', 10, 'tol', 1e-4, 'verbose', true); [theta_est, ~, hist_norm] = relax_doa(X, d_lambda, 2, opts); fprintf('估计角度: %.2f°, %.2f°, 残差比: %.4e\n', ... theta_est(1), theta_est(2), hist_norm(end));

代码里S的构造要让两个信号都保持随机相位，幅度分别1和0.1，这样功率比是20dB。噪声功率取0.01后，弱信号信噪比约0dB，强信号约20dB，正好还原“强弱悬殊”场景。跑完你会看到第一次迭代时弱信号角度可能偏到10度附近，但迭代到第4轮左右会拉回到12度附近，残差比稳定在1e-3量级。如果这一步不收敛，优先检查A_grid的归一化——a'*a必须是1，否则幅度估计s会带一个固定增益偏置。

画出扫描谱用下面这段，注意迭代结束后再扫一次残差谱，能直观看到弱信号峰是否冒出来：

% 画迭代结束后的残差扫描谱 X_res = X; for j = 1:2 idx_j = round((theta_est(j) + 90) / 0.1) + 1; a_j = exp(1j * 2 * pi * d_lambda * (0:M-1)' * sin(theta_est(j) * pi / 180)); X_res = X_res - a_j * (a_j' * X_res) / M; end P_scan = abs(A_grid' * X_res).^2 / N; figure; plot(-90:0.1:90, 10 * log10(P_scan / max(P_scan)), 'LineWidth', 1.5); xlabel('角度 (°)'); ylabel('归一化功率 (dB)'); title('RELAX残差扫描谱'); grid on;

3.3 参数表：阵元数、快拍数、搜索步长、迭代上限怎么定

下面这张表是多次实验后总结的经验范围，不同场景直接按表里选起始值再微调。

参数	推荐范围	调整建议	踩坑提示
阵元数M	8～32	源数K不超过M-2，否则自由度不足	M太小会让弱信号的谱峰被旁瓣淹没
快拍数N	5～500	实测数据少于5快拍时先做空间平滑	N=1时RELAX方差很大，不建议直接用
搜索步长grid_res	0.05～0.5度	先粗扫到1度找峰再局部加密到0.05度	步长过粗导致角度量化误差，过细拖慢收敛
最大迭代max_iter	5～15	残差比连续两轮变化小于tol即可提前停	超过15轮还在动，说明K给多了或初始化失败
收敛门限tol	1e-4～1e-2	高信噪比用1e-4，低信噪比放松到1e-2	门限太严会进入噪声拟合阶段
阵元间距d_lambda	0.5	大于0.5会引入栅瓣，小于0.5损失分辨率	非均匀阵需要改方向向量表达式

这张表的调整逻辑是：先定M和N，再定K，最后调搜索步长和迭代上限。很多人在仿真里一上来就设0.01度步长，结果A_grid变成1万多列，每次迭代多耗几十倍时间，精度却只提了零点零几度。正确姿势是先1度粗扫确定峰的大致位置，再在峰附近局部用细网格。

4. 强弱信号场景的仿真实验：验证RELAX的边界

4.1 场景一：20dB信噪比差下的强-弱源分离

把3.2的脚本跑起来，你会看到两个关键现在。第一个是迭代过程中的谱变化：第1轮扫描时，-10度强信号峰非常尖锐，12度弱信号位置只有一个微微凸起的包，肉眼几乎分不清是信号还是噪声旁瓣。顺序相消初始化会先把强信号角度定为-10度附近，幅度约等于1，相消后残差里强信号能量被干掉30dB以上，弱信号峰才开始明显。第二个现像是弱信号的收敛轨迹：初始角度可能落在10.8度，偏离真实值1度多，之后每轮迭代向12度靠拢，到第5轮基本稳住。

这个实验里MUSIC作为对照组很有意思。同样数据跑MUSIC，在12度方向要么没有峰，要么有个比强信号低18dB的小峰，而且这个小峰的宽度特别大，定位精度远不如RELAX相消后出来的峰。原因就是前面说的：强弱信号共存时，噪声子空间被强信号特征向量污染，弱信号在正交谱里的投影被压掉了。RELAX没有这个机制，它做的是时域/空域相消而不是子空间投影，弱信号能量没有被“正交感”吃掉。

需要提醒的是，RELAX在这个场景能成功有个前提——两个信号角度间隔要大于阵列的瑞利分辨极限。如果间隔缩小到5度以内，强信号的相消残留会直接盖住弱信号峰，这时候RELAX也会失效，单纯靠迭代救不回来。想分辨更近的源，需要更大的阵列孔径或者更小的阵元间距生理极限。

4.2 场景二：相干信号与低快拍数下的对比

相干信号场景是MUSIC的“职业杀手”。多径反射会让同一个信号从两个方向到达阵列，两个源完全相干，协方差矩阵的信号子空间退化成一个维度，MUSIC只能看到一个峰。RELAX的数据模型里A S N这三项都是直接面对原始快照，没有任何协方差构造，所以相干信号对它来说只是两个非线性参数，迭代照样逐个抠。

% 相干信号场景：两个入射角，同一相位中心 phi = exp(1j * 0.3); % 相对相位差 A_co = exp(1j * 2 * pi * d_lambda * (0:M-1)' * sin(theta_true * pi / 180)); s_co = [1; 0.8 * phi]; % 相干源，幅度可比 X_co = A_co * s_co * (randn(1, N) + 1j * randn(1, N)) ... + sqrt(noise_power / 2) * (randn(M, N) + 1j * randn(M, N));

这段代码里s_co是K×1的向量乘以同一个随机行向量，两个信号完全相干。跑完relax_doa，两个角度都能找回来，残差比也能收到1e-3附近。低快拍场景类似：把N改成5，MUSIC的协方差矩阵直接秩亏没法用，RELAX虽然方差变大，但角度估计的均值还在真实值附近。

相干和低快拍这两个场景合起来，就是RELAX在实测环境里最值钱的地方。实测数据往往快拍少、多径重，MUSIC在这些条件下需要额外做前向/后向空间平滑、去相干处理，而RELAX只需要常规迭代就能用。

4.3 结果怎么看：谱图、RMSE与收敛曲线

跑仿真不能只看一次估计结果，要看统计量。我一般做200次蒙特卡洛，每次重新生成随机幅度和噪声，统计RMSE。

n_trial = 200; err = zeros(n_trial, 2); for trial = 1:n_trial S = amp(:) .* (randn(2, N) + 1j * randn(2, N)) / sqrt(2); X = A * S + sqrt(noise_power / 2) * (randn(M, N) + 1j * randn(M, N)); [th, ~] = relax_doa(X, d_lambda, 2, opts); err(trial, :) = th - theta_true; end rmse = sqrt(mean(err.^2, 1)); fprintf('RMSE: 强信号 %.3f°, 弱信号 %.3f°\n', rmse(1), rmse(2));

判断结果时三条经验。第一，强信号RMSE通常在0.1度以内，弱信号RMSE在0.3到0.6度之间都算正常，如果弱信号RMSE超过1度，优先怀疑搜索步长或K给错。第二，看收敛曲线hist_norm，前2轮应快速下降，之后进入平缓期；如果曲线在第5轮突然反弹，说明某个源的角度在来回跳，这是初始值选错或网格太粗的信号。第三，对比克拉美罗界（CRB），RELAX的RMSE通常在CRB的1.5到3倍之间，如果超过5倍说明有系统性偏差，比如搜索网格量化误差没被插值去掉。

5. RELAX算法避坑指南：5个常见问题与排查

5.1 弱信号角度总往强信号方向偏

现象：弱信号真实在12度，估计出来老在10.5到11度之间，迭代好几轮都拉不回来。强信号角度估得很准，弱信号却带一个固定方向的偏置。

原因：强信号相消不彻底。网格搜索是离散的，强信号真实角度落在两个网格点之间时，估计角度会有最多半个步长的量化误差，这个误差对应的方向向量和真实方向向量不完全正交，相消后残留的强信号能量在弱信号方向形成一个小斜坡，把弱信号峰往强信号方向拽。

解决：先对已估计的强信号角度做抛物线插值，把量化误差压到0.01度量级，再做相消。插值公式在6.1节给出。另外可以把弱信号所在局部区域的搜索步长加密到0.01度，单独扫一遍，避免全局细网格带来的计算量爆炸。

提示：强信号相消后的残留功率如果能保持在噪声功率以下，弱信号偏置问题基本消失。

5.2 迭代发散：谱峰越找越歪

现象：第一轮迭代角度还算正常，第二轮开始某个源的角度跳到完全错误的位置，残差比不降反升，甚至输出全是NaN。

原因：RELAX是坐标下降，代价函数非凸，对初始值敏感。初始化如果用了“直接取前K个谱峰”的偷懒做法，两个角度可能同时落在同一个峰的相邻网格点上；主循环里更新了角度但没有立即更新幅度，相消时用的还是旧幅度，会把本来正确的残差结构破坏掉。

解决：初始化一定用顺序相消，并且主循环里每更新一个θ_k就立刻用新角度重新估计s_k。代码里我用s(k, :) = A_grid(:, idx)' * X_others，就是这个目的。另外给角度搜索加一个先验扇区，比如已知信号不会出现在负仰角区域，就把grid限制在0到90度，能避免伪峰干扰。加了扇区限制后迭代速度也会变快，因为搜索网格变短了。

5.3 信源数K给错导致谱峰乱跳

现象：K设1时，估计角度稳定在两个真实角度之间，偏向能量更强那个；K设3时，第三源的估计角度每次蒙特卡洛都不一样，残差比却降得很低。

原因：K设少了，算法被迫把两个真实源的能量合并成一个源估计，等价于拟合模型的偏差太大；K设多了，多余的源会去拟合噪声，由于每次噪声实现不同，多余源的参数自然到处乱跳。

解决：先用MDL/AIC对协方差矩阵特征值做一次粗定阶，再用2.3节的残差比拐点法二次确认。拐点法里，如果K从2增加到3时残差比下降小于20%，基本可以判定K=3是过估计。实测数据里我还会额外看一个指标：多余源估计出来的幅值是否比噪声底还低，如果低，直接把它从结果里去掉重启一轮RELAX。

5.4 单快拍下RELAX失效

现象：N=1时，强信号角度估计还行，弱信号角度方差变得极大，残差比卡在0.01附近下不去。

原因：RELAX的扫描谱P = |a'X_others|^2 / N，N=1时这个量只有一个快拍的投影能量，噪声的随机起伏没有被平均掉，弱信号峰和噪声尖峰在幅度上分不开。多快拍时N次平均能把噪声投影压到原来的1/N，单快拍没有这个平均增益。

解决：至少用3到5个快拍做累积。如果系统架构限制了每次只能出一个快拍，可以对连续几个快拍做滑动窗累积：X_acc = [X(t), X(t-1), X(t-2)]，把三个快拍拼成M×3矩阵再进RELAX。注意快拍之间的信号幅度相位要一致，否则直接拼会破坏相干性，这种情况需要先做相位校正。

5.5 网格搜索精度卡死：0.1度步长就是精度上限

现象：无论怎么加密网格，角度估计误差始终在grid_res/2附近徘徊，把步长从0.1度细化到0.01度，误差只下降了不到20%，计算时间却涨了10倍。

原因：网格搜索的本质是量化估计，输出角度只能是网格点上的值。虽然加密网格能降低量化误差，但RELAX的角度谱在峰附近不是尖锐的δ函数，而是有一定宽度的光滑峰，网格点上的最大值位置受谱形状影响，不是简单加密就能逼近真实值。

解决：粗网格全局搜索后，对峰值点做抛物线插值或高斯插值。抛物线插值只需要峰和左右两个邻点的功率值，公式在6.1节。实测中这一步能把角度误差从半个步长压到0.01度量级，代价是额外三次乘法。注意插值只在谱峰形状近似抛物线时有效，如果峰被旁瓣污染严重，插值结果不可信，这时应该先加窗减小旁瓣。

注意：RELAX本身不是超分辨算法，它的精度上限由阵列孔径和信噪比决定，网格加密和插值只是去除量化误差，不能突破物理分辨极限。

6. 把RELAX用到极致：两阶段搜索、FFT加速与联合MUSIC

6.1 两阶段搜索：粗网格兜底、抛物线插值细化

全局细网格是RELAX性能的最大杀手。我现在的标准流程是：先用1度网格扫描全角度域，找到K个粗峰；然后对每个峰在±1度局部范围内用0.01度网格扫描，最后做抛物线插值。局部细扫的角度范围小，A_grid只覆盖几十个角度，计算量比全局细扫小两个数量级。

% 抛物线插值细化角度估计 [~, idx] = max(P_local); % P_local 是局部细扫的谱 if 1 < idx && idx < length(P_local) p1 = P_local(idx-1); p2 = P_local(idx); p3 = P_local(idx+1); denom = p1 - 2 * p2 + p3; if denom ~= 0 delta = 0.5 * (p1 - p3) / denom * grid_res_local; theta_fine = theta_coarse + delta; else theta_fine = theta_coarse; end end

delta的符号决定了插值方向：分母为负时抛物线开口向下，峰值在p2左侧还是右侧由分子决定。这个公式默认谱峰是对称抛物线，实际中强信号谱峰可能略带偏斜，但作为亚网格修正已经够用。如果角度谱在峰附近被旁瓣污染，denom可能接近零或变正，这时插值结果不可信，退回到粗网格估计值。

6.2 均匀线阵下的FFT加速实现

对于均匀线阵，方向向量a(θ) = [1, e^{j2πd sinθ/λ}, …, e^{j(M-1)2πd sinθ/λ}]^T正好是空间频率域的傅里叶基。对单快拍数据做FFT，得到的空间谱可以直接对应角度谱，把网格搜索变成一次FFT。

% 空间FFT粗扫（零填充提高插值密度） X_fft = fft(X, 2^nextpow2(M * 16), 1); % 按阵元维度做FFT P_fft = sum(abs(X_fft).^2, 2); % 累积所有快拍的能量 u = (0:size(X_fft,1)-1) / size(X_fft,1) * 2 - 1; % 归一化空间频率 theta_tmp = asin(u / (2 * d_lambda)) * 180 / pi; % 映射到角度

这段代码有两点要注意。第一，FFT输出的是空间频率u从-1到1的范围，对应角度从-90度到90度，但asin的自变量要在[-1,1]内，超出部分对应栅瓣区域，要剔除。第二，FFT的空间分辨率只有2π/M，零填充只能插值不能提高真实分辨力，所以FFT只适合做初始化粗扫。我会在初始化阶段用FFT代替顺序相消里的完整网格扫描，把初始化时间从秒级压到毫秒级，主循环里再用局部网格精细搜索。

6.3 用蒙特卡洛验证：RMSE对CRB的收敛性

判断RELAX实现得好不好，终极标准是让RMSE逼近CRB。确定性信号模型下，角度估计的CRB是：

CRB(θ) = σ²/(2N) * [Re{D^H Π_A^⊥ D ⊙ (A^H A)^T}]^(-1)

其中D是方向向量对θ的偏导数矩阵，Π_A^⊥ = I - A(A^H A)^(-1)A^H是A的正交补投影矩阵。实际验证时不用手工算这个矩阵求逆，直接对比不同信噪比下RELAX的RMSE曲线：高信噪比（20dB以上）RMSE和CRB的差距应该小于3倍；差距过大，问题出在网格量化或相消残留；差距过小反而要怀疑是不是偶然事件太少，至少跑500次再下结论。

我习惯把RMSE/CRB倍率作为RELAX实现的验收指标。调试新阵列配置时，先跑这个倍率；如果倍率超5，我优先查初始化顺序相消是否把旁瓣误当信号，其次查迭代是否在角度和幅度更新之间失配。这样一步一步查，基本能收敛到正确实现。

做DOA估计这几年，我最大的教训是RELAX不是超分辨算法，它是一个把能量分干净的算法。别拿它跟MUSIC比分辨率，要比鲁棒性——强弱信号、相干源、少快拍，这些才是它的主场。把搜索步长、迭代上限和信源数三个参数当作一组整体来调，先粗后细，大概率不会翻车。希望帮到你。

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