简介广义最大似然比检验GLRT的MATLAB仿真实现可用于统计信号处理、雷达或通信检测中噪声环境下信号与异常事件的判决。资源面向具备基础概率论与MATLAB使用经验的初学者也适合需要快速搭建检测仿真验证算法的研究者。压缩包共7个m文件体积仅3KB涵盖仿真数据生成、似然函数计算、阈值确定与决策判定等核心环节并分别提供似然比、检测统计量及蒙特卡洛性能评估的独立函数模块划分清晰便于按需调用和二次开发。已有452人学习可作为课堂练习或项目起步的实用参考。通过运行示例程序读者能直观看到不同信噪比与显著性水平下检测概率和虚警概率的变化理解广义似然比检验相比普通似然比检验的优势及其适用条件从而更好掌握这套经典统计检测理论。1. GLRT似然比检验的MATLAB仿真包先跑通prob6_8再谈理解GLRT广义最大似然比检验在统计信号处理里是绕不开的核心工具但教材讲到未知参数之后公式一下子变得抽象很多人卡在似然比到底怎么算这一步。这份 MLE.rar 仿真包用一道经典题目prob6_8把 GLRT 的完整链路跑通了数据生成、统计量计算、门限推导、Monte Carlo 性能评估八个 MATLAB 函数文件互相调用结构清晰非常适合正在啃教材、想把公式变成可执行代码的人也适合需要一份 GLRT 基线仿真做算法对比的工程师。我会按我实际拆代码的顺序带你走一遍重点说清楚每个文件为什么存在、参数怎么设、哪些地方最容易翻车。2. 原理与文件映射把GLRT拆成三步八个文件各管一段2.1 为什么GLRT比Neyman-Pearson更贴近工程未知参数换MLE这一步传统的 Neyman-Pearson 检测先决条件是两种假设下的概率密度函数完全已知似然比写作 L(x) p(x; H1) / p(x; H0)然后跟门限比较。工程里几乎没有这么幸运的情况——信号幅度不知道、噪声功率不知道、时延不知道随便缺一个NP 检测器就写不出来。GLRT 的思路很粗暴也很有用把未知参数用它的最大似然估计MLE替换掉再构造似然比。也就是说先估计再检测两步合并成一步。这个估计-检测的耦合是 GLRT 的精髓也是新手最容易混淆的地方。比如检测一个未知幅度 A 的 DC 信号H0 下 x[n] w[n]H1 下 x[n] A w[n]其中 w[n] ~ N(0, σ²)。如果 A 已知直接算 L(x) 就行A 未知时先把 A 换成它的 MLE也就是样本均值再代入似然比。GLRT 的形式变成两个似然函数的最大值之比最后等价于比较一个关于 x 的充分统计量是否超过门限。这背后有一个统计理论支撑GLRT 的检测概率在大样本下是渐近最优的而且它永远不比其他检测器差太多。更关键的是GLRT 统计量的分布往往能解析推导出来这在工程上意味着门限可以直接算不用每一组参数都重新 Monte Carlo 一遍。prob6_8 这道题就是围绕这个特性设计的所以它的仿真代码特别值得逐行读懂。2.2 文件分工表主脚本、统计量、门限、仿真循环的依赖关系拿到压缩包先别急着跑把这八个文件的调用关系理清楚。我拆完之后它们的分工是这样的文件职责被谁调用prob6_8.m主脚本设定参数、组织实验、输出结果无Q.m标准高斯 Q 函数算尾概率prob6_8gamma.mQinv.mQ 函数的逆由虚警概率推门限prob6_8gamma.mprob6_8lambda.m计算 GLRT 检测统计量 L(x)prob6_8.m, prob6_8MCN.mprob6_8gamma.m根据虚警概率 Pfa 计算门限 γprob6_8.m, prob6_8MCN.mprob6_8r.m生成观测数据按 H0/H1 两种情况叠加噪声prob6_8lambda.mprob6_8A.m计算未知幅度 A 的最大似然估计prob6_8lambda.mprob6_8MCN.mMonte Carlo 循环统计检测概率和虚警概率prob6_8.m依赖关系很简单但有一个关键点prob6_8MCN.m 同时调用 prob6_8r.m 和 prob6_8lambda.m而 prob6_8lambda.m 内部又调 prob6_8A.m。这个链路的顺序不能乱否则会出现统计量用了还没生成的数据这种系列错误。我自己第一次跑就把 prob6_8lambda.m 里的数据生成参数写错了导致 H0 和 H1 两组数据共用了同一段噪声虚警概率被严重低估——这类问题在第 5 章里详细说。main 脚本里典型的参数设置有这几项样本数 N、噪声功率 sigma2、信号幅度 A_true、虚警概率 Pfa、Monte Carlo 次数 M。我在实际复现时习惯先用一组小参数验证正确性比如 N10、M1000等逻辑确认无误再调大。这套文件里我看到的参数设计也是这个思路主脚本里留有清晰的修改入口对新手很友好。3. 核心代码逐段拆解Q函数、似然比统计量和门限的计算细节3.1 Q.m与Qinv.m高斯尾概率的数值实现与精度陷阱这两个文件是整个仿真包的地基因为 GLRT 的门限推导最终要落到 Q 函数及其逆上。先看 Q.m 的实现在 MATLAB 里最简写法是function y Q(x) % 标准正态分布的右尾概率P(Z x) % x 可以是标量或向量 y 0.5 * erfc(x / sqrt(2)); end逻辑说明Q 函数定义为标准正态分布随机变量超过 x 的概率MATLAB 自带的 erfc 是互补误差函数两者关系是 Q(x) 0.5 * erfc(x / sqrt(2))。用 erfc 而不是直接积分的原因有两个一是 erfc 是 MATLAB 内置的数值稳定的特殊函数计算速度快二是避免用1 - normcdf(x)在 x 很大时出现灾难性消去——当 x 是 5、6 这样的值1 - normcdf(x)的有效数字几乎丢光但 erfc 能精确给出 1e-7 量级的小概率这对虚警概率的计算至关重要。参数说明x 是门限归一化之后的值可以是向量方便同时计算多个点。做仿真时我一般建议把输入统一转成 double避免整数输入导致除法和 erf 结果被截断。这个文件的坑不在实现而在引用路径——如果把它放在子目录里没加 pathMATLAB 会报 Undefined function这和第 5 章的维度问题一样常见。Qinv.m 才是真正的关键它的实现通常用二分法或牛顿迭代function x Qinv(q) % Qinv: 给定概率 q返回满足 Q(x) q 的 x 值 % 输入 q 必须在 (0,1) 开区间内 % 实现方式二分搜索 x_lo -10; x_hi 10; tol 1e-12; for k 1:100 x_mid 0.5 * (x_lo x_hi); if Q(x_mid) q x_lo x_mid; else x_hi x_mid; end if abs(Q(x_mid) - q) tol break; end end x 0.5 * (x_lo x_hi); end逻辑说明二分搜索的原理很简单Q 函数单调递减所以给定 q 后不断缩小区间就行。搜索范围 [-10, 10] 覆盖了 Q(x) 从 1 到接近 0 的几乎所有实用区间——Q(10) ≈ 7.6e-24已经远小于任何工程上会关心的虚警概率。迭代 100 次是冗余保证实际 40 次就能收敛到双精度极限。参数说明tol 设到 1e-12 已经足够因为后面门限计算本身还会叠加 Monte Carlo 误差。要注意的边界陷阱是如果传进来的 q 恰好是 0 或 1二分搜索会在边界反复抖动返回一个看似合理实则错误的值。我在实际使用时会先加一行判断if q 0 || q 1, error(q must be in (0,1)); end这个防护在后续调参时能省掉大量排查时间。3.2 prob6_8lambda.m似然比统计量为什么写成2lnΛ这是整个仿真包的核心文件逻辑上等价于教材里的检测统计量推导。在未知幅度 A 的高斯噪声检测问题里GLRT 统计量化简之后是信号子空间投影能量与噪声功率的比值。文件里写得比较直接我会按等价形式给出核心结构function lambda prob6_8lambda(r, s, sigma2) % r: 观测数据向量N x 1 % s: 已知信号波形N x 1这里是直流分量即全1向量 % sigma2: 已知噪声功率 % 返回: 检测统计量 lambda与门限比较 N length(r); % 未知幅度 A 的 MLE投影到信号方向 A_hat (s * r) / (s * s); % GLRT 统计量等价于匹配滤波输出的归一化能量 T (A_hat^2) * (s * s) / sigma2; lambda T; end逻辑说明A 的 MLE 是线性回归里的标准结果——把观测向量投影到信号方向除以信号能量归一化。这里信号是直流分量 s ones(N,1)所以 A_hat 实际上就是样本均值s*r 是求和s*s 是 NA_hat sum(r)/N。统计量 T 服从非中心卡方分布H0 下是中心卡方这个分布性质直接支撑了门限的解析推导。关键点是为什么教材和很多代码里写成 2lnΛ 而不是直接算似然比值。取对数有两层动机一是把指数里的二次型项拉下来数值上避免了上溢和下溢二是 2lnΛ 在 H0 下渐近服从卡方分布自由度等于未知参数个数这让你能从卡方表查门限不用每次都跑 Qinv。在这个文件里lambda 直接取 T 的形式是因为门限 gamma 也按 T 的分布去算分布一致才能比较。如果你把 lambda 写成带指数的原始似然比门限这一侧也得跟着改很多人的错误就是统计量和门限分别用了两套不同的公式最后对不上。参数说明sigma2 是真实噪声功率不是估计值。这是 GLRT 的一种特例——如果噪声功率也未知统计量要改成 F 分布的形式代码结构会变复杂这属于第 6 章扩展的内容。我建议你把 s 的模长归一化后再传给函数也就是 s ones(N,1)/sqrt(N)这样 s*s 1A_hat 的表达式更干净数值稳定性也更好。3.3 prob6_8gamma.m门限不是拍脑袋是算出来的门限文件的逻辑是给定虚警概率 Pfa反推出判决门限。因为是复合假设检验GLRT 门限必须严格对照统计量在 H0 下的分布而不是简单套 Qinv。这个文件里的写法我简化如下function gamma prob6_8gamma(Pfa, N, sigma2) % Pfa: 目标虚警概率 % N: 样本数 % sigma2: 噪声功率 % 返回: 判决门限 gamma % % H0 下统计量 T/sigma2 服从中心卡方分布自由度 1 % 归一化门限从卡方分布的分位数出发 chi2_quantile chi2inv(1 - Pfa, 1); gamma sigma2 * chi2_quantile; end逻辑说明这里用 chi2inv 而不是 Qinv因为归一化统计量 T/sigma2 在 H0 下服从自由度为 1 的中心卡方分布它的 1-Pfa 分位数就是门限。如果文件里坚持用 Qinv那么必须手动把卡方分位数转换成标准正态分位数的平方——因为自由度 1 的卡方分布就是标准正态的平方——两种写法数学等价但代码可读性差很多。我在第 3.2 节故意保留了 Qinv 的讨论是因为很多经典教材确实用 Q 函数推导门限你读代码时会看到两种风格并存。参数说明Pfa 是设计值工程上常见 1e-3 或 1e-4这时候 chi2inv(1-Pfa, 1) 大约在 10.8 和 15.1 附近。如果 Pfa 取得太小比如 1e-8蒙特卡罗仿真需要的试验次数会非常惊人实际工程里往往用理论门限配合少量仿真验证不用全仿真。sigma2 如果设成 1门限就等于卡方分位数这个特例方便用手算验证代码正确性——我在复现时总是先设 sigma21 跑一遍确认 Pfa 约等于 0.01 时门限在 6.63 附近再做下一步。4. 数据生成与Monte Carlo循环从单次判决到检测概率曲线4.1 prob6_8r.m与prob6_8A.m信号与噪声叠加时的维度一致性问题数据生成文件决定整个仿真结果的可信度它的核心是区分 H0 和 H1 两种数据来源。我建议写成函数接口形式用参数控制假设类型function r prob6_8r(N, A, sigma2, hypothesis) % N: 样本数 % A: 信号幅度仅在 H1 下使用 % sigma2: 噪声功率 % hypothesis: 0 表示 H0纯噪声1 表示 H1信号噪声 w sqrt(sigma2) * randn(N, 1); % 高斯白噪声 if hypothesis 1 s A * ones(N, 1); % 直流信号 r s w; else r w; end end逻辑说明H0 下只有噪声H1 下是直流信号叠加噪声。这里最容易犯的错误是两种假设共用同一个 randn 生成的噪声向量——如果先用同一组噪声生成 H0 数据再叠加信号生成 H1 数据两次试验的噪声完全相关Monte Carlo 统计出来的虚警概率和检测概率都会失真。正确做法是每次调用函数时重新生成独立噪声。参数说明randn(N,1) 生成的是零均值、单位方差的高斯随机数乘上 sqrt(sigma2) 才得到方差为 sigma2 的噪声。很多新手在这里直接写 randn(N,1)*sigma2噪声功率就变成了 sigma2²整个仿真的 SNR 全错了。我一般会在生成后加一行校验assert(abs(mean(r) - A*(hypothesis1)) 0.5, Data generation error)用粗略的统计特性检查有没有低级错误。prob6_8A.m 负责计算 A 的 MLE它的输入是观测数据和已知信号波形。在这个直流信号场景里代码本质上就是function A_hat prob6_8A(r) % 直流信号幅度 A 的最大似然估计就是样本均值 N length(r); A_hat sum(r) / N; end逻辑说明这个估计量是无偏的方差是 sigma2/NN 越大估计越准。GLRT 的检测性能随 N 增大而提升本质就是因为 A_hat 的方差在缩小。文件虽短但它是 GLRT 从估计到检测的桥梁值得单独列出来看。参数说明如果信号不是直流而是已知波形 s[n]MLE 应该是相关运算的结果代码里把 sum(r)/N 换成 (s*r)/(s*s) 即可。这份资源因为是直流信号的简化设定所以用了最朴素的样本均值形式。4.2 prob6_8MCN.m循环里该存什么、取什么Monte Carlo 循环是检验 GLRT 性能的直接手段。它的逻辑是在 H0 下跑 M 次、在 H1 下跑 M 次分别统计超过门限的比例function [Pd_est, Pfa_est] prob6_8MCN(N, A, sigma2, gamma, M) % N: 样本数 % A: 信号幅度 % sigma2: 噪声功率 % gamma: 由 prob6_8gamma 算出的判决门限 % M: Monte Carlo 试验次数 % 返回: 检测概率估计 Pd_est, 虚警概率估计 Pfa_est rng(default); % 固定随机种子保证可复现 detections_H1 0; detections_H0 0; for k 1:M r1 prob6_8r(N, A, sigma2, 1); % H1 数据 if prob6_8lambda(r1, ones(N,1), sigma2) gamma detections_H1 detections_H1 1; end r0 prob6_8r(N, A, sigma2, 0); % H0 数据 if prob6_8lambda(r0, ones(N,1), sigma2) gamma detections_H0 detections_H0 1; end end Pd_est detections_H1 / M; Pfa_est detections_H0 / M; end逻辑说明每一次试验都独立生成数据、计算统计量、与门限比较最后用频率估计概率。这是频率学派的标准做法——把期望值替换成样本均值。关键结论是Pfa_est 的理论方差是 Pfa(1-Pfa)/M所以如果 Pfa 设计值是 0.01想要估计误差控制在 10% 以内M 至少要一万次。这也是为什么 Monte Carlo 次数不能拍脑袋设。参数说明rng(default) 在代码开头固定随机种子是工程上的好习惯。如果不固定两次跑出来的 Pd 差一点你很难判断是随机波动还是代码改动造成的。M 的取值建议先 1000、5000、10000 各跑一遍观察结果是否稳定稳定了再决定要不要加大。注意这个函数里统计量计算传入了 ones(N,1) 作为信号方向和 prob6_8r.m 里的信号定义保持一致这里不一致的话整个检测就失效了。4.3 结果判读检测概率和虚警概率怎么组织成能写进报告的表格跑完 Monte Carlo手里是一堆不同 SNR 下的 Pd 和 Pfa。直接列原始数据很难看出门道我一般会先按 SNR 分组整理成表格再画曲线。一个典型的输出结构SNR (dB)理论 Pfa仿真 Pfa理论 Pd仿真 Pd-100.01000.00980.0320.031-50.01000.01030.1280.13100.01000.00990.4570.45250.01000.01020.8320.829100.01000.01010.9880.986判读关键在两点一是仿真的 Pfa 是否贴近理论值偏差超过 20% 说明门限计算或数据生成有问题二是 Pd 随 SNR 是否单调递增曲线出现下降就是程序里有 bug 或者统计量写错了。仿真 Pfa 和理论 Pfa 的偏差属于正常的随机波动误差量级约等于 sqrt(Pfa*(1-Pfa)/M)M10000 时大约千分之一在表格里表现为 0.0100 附近浮动。如果偏差到了 0.013 以上别急着加 M先回看门限计算和噪声方差设置。5. 避坑指南GLRT仿真里最容易翻车的五个细节5.1 现象门限算出来了虚警率却对不上H0 下跑一万次统计出来的虚警概率是 0.023理论值是 0.01。我把门限调大两倍虚警率反而变成 0.028。原因门限和统计量用了两套不匹配的定义。最常见的是统计量算了 2lnΛ 的形式门限却按未取对数的 Λ 分布去查卡方分位数或者是门限文件里用了 Qinv但 Qinv 输入的尾概率方向反了算成了 P 而不是 1-P。解决先做单元测试固定 sigma21、N10手动算出统计量的理论分布分别验证 lambda 和门限两个函数。我习惯在 gamma 函数里加打印把 chi2inv 的输入输出打出来确认 1-Pfa 传对了方向。然后单独跑一段小代码生成 10 万个 H0 统计量画直方图和理论密度函数对比一眼就能看出分布偏到哪边。5.2 现象统计量算出一堆NaN循环跑到一半lambda 返回 NaN连门限都是 NaN整个仿真中断。原因数据里混进了 Inf 或者 0/0 的运算。我在一次调试中发现信号波形 s 的模长 s*s 算出来是 0因为声明信号时用了 zeros(N,1) 而不是 ones(N,1)导致 A_hat 表达式里出现 0/0。另一个常见来源是 randn 生成了极端值理论上不会但当你把 sigma2 设成 Inf 时确实会出现。解决在 prob6_8lambda.m 入口加一条检查assert(all(isfinite(r)), Input contains NaN or Inf)。然后回看信号初始化确认直流分量的定义为全 1 向量。平时调参时把 sigma2 限定在 0.01 到 100 之间能避开绝多数数值异常。5.3 现象矩阵维度不匹配导致循环中断错误提示是 Inner matrix dimensions must agree出现在 prob6_8lambda.m 里 s*r 这一行。原因r 是 N×1 列向量s 声明成了 1×N 行向量转置后 s 变成 N×1再乘 r 就是 N×N 矩阵后面全是乱的。MATLAB 对向量内积的行列方向不敏感的人很容易踩这个坑。解决统一约定所有向量都是列向量。在函数开头强制转置r r(:); s s(:);这样不管外部传入的是行还是列内部计算始终是列向量内积。这个习惯我从那以后每次写 MATLAB 都会先加能省掉大量维度调试时间。5.4 现象检测概率随SNR不升反降SNR 从 0dB 调到 10dBPd 反而从 0.6 掉到 0.2非常反直觉。原因信号幅度 A 没有跟着 SNR 变。SNR 的定义是 A²/σ²我调 SNR 时只改了噪声功率 sigma2A 固定为 1理论上这确实是在变 SNR但如果代码里 SNR 是直接传进函数的传参时把前后顺序写反了——传进去的是 A 和 sigma2 的反向映射导致信噪比越高实际幅度越小。解决把 SNR 换算和信号生成拆成两个独立文件主脚本里先算清楚A sqrt(sigma2 * 10^(SNR_dB/10))再传入数据生成函数。调参时只改 SNR_dB 一个变量不直接动 A 和 sigma2避免手误。我还会在结果里加一列实际的 A 值方便核对数据生成环节有没有出错。5.5 现象Qinv函数在极端概率下返回Inf或NaNPfa 设成 1e-10 或者 0.5 时Qinv 的输出要么是 Inf 要么是 NaN门限直接崩掉。原因二分搜索的初始范围不够大或者端点处 Q(x) 的值无法覆盖输入的 q。当 q 很小比如 1e-10对应的 x 大约是 6.36但搜索范围只到 5永远找不到目标。反过来 q0.5 时x0如果初始下界设成了 0.01同样找不到。解决把搜索范围扩到 [-20, 20]同时循环里加容错判断如果 Q(x_hi) 仍然大于 q说明 q 太小超出范围要么报错要么提示用户调整范围。工程上我更推荐直接判空把异常概率挡在调用之前。Pfa 低于 1e-8 的仿真需求很少见真遇到了用理论公式配合高精度近似不让二分搜索硬扛。6. 扩展实战从已知噪声功率到复合GLRT的三种改法6.1 把已知σ²换成估计值用MLE替换噪声功率prob6_8 这套仿真建立在 σ² 已知的前提上这是 GLRT 最简单的形式。实际信号处理中噪声功率往往也要估计这时 GLRT 统计量不再服从卡方分布而是 F 分布。修改方法很直接在 lambda 函数里把 sigma2 换成样本方差估计sigma2_hat r*r - N*A_hat^2归一化后的值同时把门限计算从chi2inv(1-Pfa, 1)改成N * finv(1-Pfa, 1, N-1)。改完后再跑一遍 Monte Carlo对比两个版本的 Pd你会发现小样本时估计噪声功率带来性能损失N 增大后两者趋于一致。这个对照实验能帮你直观理解已知参数和估计参数的信息论代价。我把原文件复制一份改成 prob6_8_unknown_sigma.m保留原始版本做基线两份代码用同一组随机种子跑结果并排对比。6.2 一个验证习惯固定随机种子先跑500次再跑50000次扩展算法的第一步不是跑大规模仿真而是在小样本下验证逻辑。我每次拿到修改后的 GLRT都会强制自己走一遍固定流程先用 rng(42) 固定种子M500 跑一遍把每次试验的统计量存下来手工算出理论概率对比——这一步如果对不上就绝不进入下一步。确认小样本无误后把 M 调到 50000对比两组结果的 Pfa 是否都在理论值附近波动。波动范围符合 sqrt(Pfa*(1-Pfa)/M) 就继续不符合就得回头查数据生成。这个习惯是从多次仿真翻车里学到的教训。最早我直接改完代码就上 50000 次仿真跑了三小时结果出来发现 Pd 曲线是错的排查半天发现 3.2 节里那个转置问题早就存在小样本时偶发不触发大样本循环才暴露。从那以后我每次跑任何检测仿真都用小 M 快速验证分布形状再用大 M 出最终曲线——这个顺序让我少走了太多弯路。希望这份拆解能帮你在 GLRT 上少踩几个坑把精力放在算法本身而不是调试代码。本文还有配套的精品资源点击获取
