简介面向MATLAB信号处理学习者的最大似然波达角估计实用脚本聚焦于利用最大似然估计MLE算法从多天线接收信号中解析波达角适用于无线通信、雷达探测与声纳系统等领域的学习与研究。脚本涵盖数据预处理、模型设定、似然函数构建、优化求解及误差计算等完整流程通过相位差模型将波达角与观测数据关联并借助fminunc等优化工具进行迭代估计帮助读者掌握MLE在阵列信号处理中的具体实现与性能评估方法。最大似然估计作为经典参数估计方法其在波达角问题中的代码实现具有较强的工程参考价值。包内仅有1个m文件文件大小仅2KB无需额外数据即可运行适合作为入门或课程设计的参考模板。已有184人学习对于希望快速理解最大似然估计原理及波达角估计代码结构的开发者是一份轻量但功能完整的示例。1. 当MUSIC失效时最大似然估计为什么值得自己写一遍在很多阵列信号处理项目里MUSIC、ESPRIT这类子空间算法是默认首选因为它们算得快、不用迭代。但一旦遇到低信噪比、少快拍或者相干信源子空间分解会直接丢掉信号子空间维度谱峰分裂估计偏差非常大。这时候最大似然估计MLE的价值就体现出来了它在模型正确的前提下是渐进最优的方差能逼近克拉美罗界。Maximum_likelyhood_estimation.m就是这样一个用MATLAB从零实现最大似然波达角估计的完整脚本覆盖数据生成、似然函数构造、fminunc优化和MSE评估。适合做雷达测向、声纳定位和5G阵列通信的工程师也适合写课程设计时想弄懂MLE而不是只调工具箱的人。这篇笔记我会从阵列模型讲到集中似然再给出可直接跑的MATLAB代码和性能验证方法最后把我调试时踩过的优化器坑一并说清。2. 均匀线阵模型与集中似然函数MLE为什么能逼近克拉美罗界最大似然估计的第一步不是写优化代码而是把接收信号模型写到能求导的程度。这里考虑最常见的均匀线阵ULA假设有M个阵元阵元间距d信号波长为lambda有一个远场窄带信号从角度theta入射。那么第t个快拍的接收向量可以写成x(t) a(theta) * s(t) n(t)其中a(theta)是导向矢量s(t)是信号复幅度n(t)是零均值复高斯白噪声。这个模型看起来简单但它决定了后面所有推导导向矢量里只有角度这一个未知参数而幅度、噪声功率都算扰动项。MLE的思路就是找到一个theta让观测数据x的联合概率密度最大。2.1 导向矢量与阵列流型矩阵的实现对于ULA第m个阵元相对参考点的相位延迟是2*pi*d*(m-1)*sin(theta)/lambda所以导向矢量写出来是一个复指数向量。工程实现时要注意角度要转成弧度阵元间隔以波长为单位时代码里可以直接用sin(theta)乘以系数不需要真的去计算每条传播路径的延迟。下面这段MATLAB代码生成单信号源的导向矢量function a steering_vector(theta_deg, M, d_over_lambda) % theta_deg: 入射角单位度 % M: 阵元数量 % d_over_lambda: 阵元间距除以波长通常取0.5 theta deg2rad(theta_deg); m (0:M-1).; a exp(1j * 2 * pi * d_over_lambda * m * sin(theta)); end这段函数里m是列向量sin(theta)是标量两者相乘得到M维相位差向量。exp(1j * ...)生成复指数导向矢量它是后续所有似然函数计算的基础。如果信源不止一个就把多个导向矢量横向拼成矩阵A [a(theta1), a(theta2), ...]这个矩阵叫阵列流型矩阵。在最大似然估计中我们最终要优化的是这个矩阵里包含的角度参数。2.2 集中似然把多维优化降成一维直接对theta、信号幅度S和噪声功率sigma^2求联合最大似然参数太多优化不现实。标准做法是先用最小二乘解把线性参数消掉。对固定theta信号幅度的最大似然解是s_hat (a^H * a)^(-1) * a^H * x对于单信源a^H * a就是个复数值等于M所以a^H * x就是匹配滤波输出。把S的解析解代回原似然函数就得到只关于theta的集中对数似然函数。忽略常数项后最大化似然等价于最大化J(theta) x^H * a * (a^H * a)^(-1) * a^H * x这个表达式在DOA文献里常写成tr(P_A * R_xx)其中P_A A*(A^H*A)^(-1)*A^H是导向矢量张成空间的投影矩阵R_xx是样本协方差矩阵。为什么要用多个快拍呢单快拍噪声影响太大用N个快拍平均得到协方差矩阵后投影能量能被平滑出来。下面的代码计算负对数似然代价函数方便后面用fminunc求最小值function val neg_log_likelihood(theta_deg, X, M, d_over_lambda) % X: M x N 复数接收矩阵N为快拍数 % 返回标量代价优化时需要最小化该值 a steering_vector(theta_deg, M, d_over_lambda); P_A a * (a * a)^(-1) * a; % 投影矩阵 R_xx X * X / size(X, 2); % 样本协方差矩阵 val -real(trace(P_A * R_xx)); % 取负号是因为fminunc找最小值 end这里P_A是M维复方阵trace(P_A * R_xx)计算信号在导向矢量方向上的能量。由于我们要求最小值所以把要最大化的目标函数取负。real()是为了把浮点运算中微小的虚部误差去掉不影响梯度方向。注意a是共轭转置而不是普通转置这是复信号处理最容易出错的地方。下表列出了集中似然推导中几个关键参数的角色弄混了会在优化时得到奇怪的结果参数含义维度/取值影响M阵元数量正整数导向矢量长度决定阵列孔径N快拍数正整数协方差矩阵估计质量d_over_lambda阵元间距与波长比通常0.5过大产生栅瓣过小降低分辨率P_A导向投影矩阵M x M集中似然的核心消去信号幅度R_xx样本协方差矩阵M x M由快拍数据统计得到在低快拍场景下R_xx的估计误差会直接传导到代价函数里导致MLE出现多个局部极值。这也是为什么后面必须做多起点搜索而不是只调用一次fminunc。3. 从仿真数据到fminunc求解最大似然角估计的完整MATLAB实现这一章直接给出能跑的脚本。设计思路是先用已知角度生成仿真数据再用最大似然估计把它还原出来最后和真实值对比。这样你能清楚看到每一步数据在哪、参数怎么改。3.1 生成带噪接收数据仿真数据质量直接影响优化收敛。这里设置M8个阵元目标信号来自theta_true10度快拍数N200信噪比设为10dB。信号幅度用复高斯随机变量噪声是独立复高斯白噪声。代码如下rng(0); % 固定随机种子方便复现 M 8; % 阵元数量 N 200; % 快拍数 theta_true 10; % 真实波达角度 d_over_lambda 0.5; % 半波长间距 snr 10; % 信噪比dB a steering_vector(theta_true, M, d_over_lambda); s sqrt(10^(snr/10)) * (randn(1, N) 1j*randn(1, N)) / sqrt(2); n (randn(M, N) 1j*randn(M, N)) / sqrt(2); X a * s n;这里信号幅度sqrt(10^(snr/10))是把dB信噪比转成线性幅度比。噪声功率归一化为1所以信号功率直接由这个幅度系数控制。randn(1,N)1j*randn(1,N)得到复高斯信号除以sqrt(2)保证实部和虚部总功率为1。X的每一列是一个快拍每一行是一个阵元的采样。3.2 一维网格扫描加fminunc精估计一维角度搜索的初值比二维、多维场景容易处理。我的做法是先做粗网格扫描用扫描结果作为fminunc的初始点。这样比随机初值稳定得多也不会陷到远离真实角度的局部极值。下面这段代码完成网格扫描和精估计theta_grid -90:0.5:90; % 粗网格 cost_grid zeros(size(theta_grid)); for k 1:length(theta_grid) cost_grid(k) neg_log_likelihood(theta_grid(k), X, M, d_over_lambda); end [~, idx] min(cost_grid); theta_init theta_grid(idx); % 网格最优值作为初值 options optimoptions(fminunc, ... Algorithm, quasi-newton, ... Display, off, ... MaxIterations, 200, ... OptimalityTolerance, 1e-8, ... StepTolerance, 1e-8); theta_hat fminunc((t) neg_log_likelihood(t, X, M, d_over_lambda), ... theta_init, options); fprintf(估计角度: %.4f deg, 真实角度: %.4f deg\n, theta_hat, theta_true);这段代码里theta_grid步长0.5度网格扫描的协方差矩阵重复计算了181次对于单信源MLE完全够快。fminunc的Algorithm选quasi-newton因为代价函数光滑不需要提供解析梯度。MaxIterations设200次对一维问题足够了。如果优化结果停在网格初值附近说明网格分辨率不够需要加密。3.3 优化器参数怎么选下表是我在调试这个脚本时常用的参数组合不同值会影响收敛速度和精度但不是越大越好参数推荐值作用调试要点Algorithmquasi-newton用BFGS更新Hessian近似对一维平滑问题稳定MaxIterations200限制迭代次数太小会提前停OptimalityTolerance1e-8梯度范数阈值太严会拖慢结束StepTolerance1e-8参数变化阈值配合上面的值保持一致Displayoff不打印每次迭代批量跑时避免刷屏如果你的MATLAB没有优化工具箱也可以换成一维黄金分割搜索或fminbnd。fminbnd不需要初值只需要角度区间但收敛速度比fminunc慢一点。实际测试中fminunc在半波长间距、单信源情况下基本无偏偏差主要来自有限快拍和噪声实现。提示不要用fminunc默认的trust-region算法因为它要求目标函数返回梯度值这里我们只给了函数值会直接报错。4. 用MSE和克拉美罗界评估最大似然估计的性能写完了优化流程下一步是量化估计算法到底准不准。只看一次运行结果没有说服力必须做蒙特卡洛统计。我会在这个章节给出MSE计算代码和克拉美罗界CRB的数值实现方便你画性能曲线。4.1 蒙特卡洛MSE统计固定theta_true10度信噪比从0dB到20dB每个信噪比下重复L200次独立实验记录每次估计结果最后计算均方根误差RMSE。脚本主体是一个大循环注意每次实验都要重新生成信号和噪声不能复用同一组数据。代码如下L 200; snr_list 0:5:20; rmse_list zeros(size(snr_list)); for snr_idx 1:length(snr_list) snr snr_list(snr_idx); err_sum 0; for trial 1:L % 重新生成数据和噪声 s sqrt(10^(snr/10)) * (randn(1, N) 1j*randn(1, N)) / sqrt(2); n (randn(M, N) 1j*randn(M, N)) / sqrt(2); X a * s n; % 网格初值 cost_grid arrayfun((th) neg_log_likelihood(th, X, M, d_over_lambda), theta_grid); [~, idx] min(cost_grid); theta_init theta_grid(idx); % 精细优化 theta_hat fminunc((t) neg_log_likelihood(t, X, M, d_over_lambda), ... theta_init, options); err_sum err_sum (theta_hat - theta_true)^2; end rmse_list(snr_idx) sqrt(err_sum / L); end disp(table(snr_list., rmse_list., VariableNames, {SNR_dB, RMSE_deg}));这里arrayfun遍历整个网格比for循环简洁但二者等价。err_sum累加的是角度误差平方最后除以L再开方得到RMSE。当你运行这段代码时会看到RMSE随信噪比升高而下降但在高信噪比段会进入一个平台这个平台通常是由网格扫描分辨率或优化容差导致的而不是MLE本身的问题。4.2 克拉美罗界的数值实现最大似然估计的理论下界是CRB。单信源ULA的CRB有闭式表达式可以直接计算。这里我给出一个数值实现方便和蒙特卡洛结果画在同一张图上function crb compute_crb(theta_deg, M, N, snr, d_over_lambda) theta deg2rad(theta_deg); m (0:M-1).; % 导向矢量对角度的一阶导数 a exp(1j * 2 * pi * d_over_lambda * m * sin(theta)); a_der 1j * 2 * pi * d_over_lambda * m .* cos(theta) .* a; snr_linear 10^(snr/10); % 单信源CRB公式 crb 1 / (2 * N * snr_linear * (norm(a_der)^2 - abs(a*a_der)^2 / M)); endCRB公式里的分母由两部分组成导向矢量导数的能量norm(a_der)^2以及投影到导向矢量方向后剩下的部分。注意这里的a*a_der是复数内积abs()取模后平方。利用恒等式a*a M可以简化成a_der在导向矢量正交方向的能量。这个值越大CRB越小说明阵列对该方向的角度敏感度越高。4.3 扫信噪比时的性能图表在matlab里通常用semilogy画RMSE和CRB曲线纵轴用对数刻度。下面是一个快速绘图代码crb_list arrayfun((s) compute_crb(theta_true, M, N, s, d_over_lambda), snr_list); figure; semilogy(snr_list, rmse_list, o-, LineWidth, 1.5); hold on; semilogy(snr_list, sqrt(crb_list), r--, LineWidth, 1.5); xlabel(SNR (dB)); ylabel(RMSE / sqrt(CRB) (deg)); legend(MLE RMSE, sqrt(CRB)); grid on;用semilogy而不直接用plot是因为RMSE在低信噪比可能是好几度在高信噪比会小于0.1度线性坐标下会压扁。sqrt(crb_list)是把CRB转成角度标准差与RMSE同量纲才能对比。如果MLE结果在高信噪比区域没有贴着CRB曲线通常问题出在初值网格没有覆盖全局最优或者优化器在平坦区域提前停止。下面是典型输出格式实际数值取决于你的随机种子和参数设置但趋势是一致的SNR (dB)RMSE (deg)sqrt(CRB) (deg)01.2e-18.1e-254.9e-23.6e-2101.5e-21.3e-2158.7e-37.9e-3可以看出当信噪比大于5dB时MLE已经接近克拉美罗界继续增加信噪比误差下降速度变缓。这个现象在波达角估计中非常典型。5. 工程落地初值、多起点搜索与参数校验最后这部分分享几个我从调试这个脚本里总结出来的实用技巧特别是当你把单信源扩展成多信源或者低信噪比场景时这些小改动会直接决定算法是否收敛。5.1 多起点搜索解决局部极值问题集中似然函数在低信噪比时会出现很多毛刺单次fminunc容易停在旁瓣位置。我的办法是取网格扫描中最小的三个局部极小值点作为初值分别运行fminunc最终取代价最小的结果。实现时不需要同时开线程顺序跑就行因为一维优化很快。大致伪代码如下[~, idxs] islocalmin(cost_grid); local_min_idx find(idxs); [~, sort_idx] sort(cost_grid(local_min_idx)); best_starts theta_grid(local_min_idx(sort_idx(1:min(3, length(sort_idx))))); best_theta NaN; best_cost inf; for start best_starts [th, val] fminunc((t) neg_log_likelihood(t, X, M, d_over_lambda), start, options); if val best_cost best_cost val; best_theta th; end endislocalmin是MATLAB R2017b之后提供的函数它会返回逻辑索引。取出局部极小值后按代价排序只取前三个。这个逻辑与全局最优点可能被局部极小值包围时特别有效。注意fminunc返回的val要与网格扫描的代价一致因为代价函数已经取了负号所以这里用最小值比较是合理的。5.2 高精度网格和优化容差的配合网格扫描步长和优化容差不是无关的。如果你把网格步长设置为0.1度但StepTolerance仍然是1e-8后续优化会做很多次微小迭代其实没有意义。我一般让网格步长在0.5度到1度之间然后让fminunc在最优网格附近做局部精化。反过来如果网格步长太粗比如2度遇到靠近0度的曲线平坦区初值误差会让fminunc多跑几十步甚至漂移到另一个旁瓣。经验值半波长阵元间距、8阵元、单信源时网格步长1度足够阵元数增加到16时可以放大到1.5度。5.3 用零度附近的角度测试初始化一个很实用的调试手段是把theta_true设成0度查看网格扫描代价曲线。因为0度附近sin(theta)对角度变化最敏感但代价函数关于0度对称如果初值落在正负两侧优化结果可能都是对的。利用这个特性可以验证你的梯度方向和坐标定义是否正确如果0度时估计值经常跳变到±90度说明导向矢量的相位因子符号写反了或者d_over_lambda用了负数。这种错误在MLE脚本里非常隐蔽因为MUSIC可能对这种符号不敏感但fminunc对梯度方向是敏感的。另外当你要把这段脚本用于实测数据时别忘了先校准阵列的幅度和相位一致性。仿真里我们假设每个阵元通道响应完全一致实际系统里会有通道失配这时最大似然估计会失真但不会完全失效。一个临时对策是在估计前对X做幅度归一化让每个阵元的平均功率相同然后用归一化后的数据走同一套似然函数。这个技巧能救回很多硬件误差但也只能作为权宜之计严谨的做法还是用校准源估计出各通道增益相位。本文还有配套的精品资源点击获取