简介这份MATLAB源码包围绕参数已知条件下的广义似然比检验GLRT信号检测问题编写适合信号处理、统计检测课程学习者与算法入门者参考。资源共4个文件压缩包仅7KB包含3个m脚本和1张结果图drawpd.m用于绘制两种假设下的概率密度函数Q.m与Qinv.m涉及Q函数及逆Q函数可用于计算检测阈值和误警概率结果.bmp直观展示了仿真输出的检测性能。目前已有1107人学习下载。通过研读和运行源码可系统掌握GLRT从模型建立、对数似然比构造到阈值判决的完整流程同时学习如何在MATLAB中实现统计检测仿真、绘制PDF曲线并评估检测效果。该小体积实例结构清晰、便于调试适合作为课程实验或自学GLRT算法的起步模板。1. 为什么参数已知时 GLRT 是信号检测的首选工具做信号检测的人大概都遇到过这样的场景噪声统计特性已知信号幅度、相位、到达时间也由先验信息给定唯一不确定的就是“此刻到底有没有信号”。这时如果还去套能量检测或者匹配滤波的固定判决门限往往会在低信噪比下损失检测性能而 GLRT广义似然比检验恰好能在这个前提下把检测问题转化为一个可解析的似然比比较过程。它不要求对未知参数做积分或贝叶斯平均而是用最大似然估计值代入似然比从而把复合假设检验简化成单边阈值判决。对做雷达目标检测、通信同步头捕获、生物医学信号事件判别的工程师来说这套 MATLAB 仿真代码的价值在于它把从假设建模、Q 函数查阈值到蒙特卡洛画曲线的完整链路都放在你面前既能验证理论曲线又能直接改成自己的数据格式。2. GLRT 检测原理与假设检验建模2.1 从 Neyman-Pearson 到 GLRT 的演进经典检测理论里当噪声是零均值高斯白噪声、信号波形完全已知时最优检测器是匹配滤波加固定阈值也就是 Neyman-Pearson 准则下的似然比检验。似然比定义为两个假设下观测数据概率密度函数的比值L(x) p(x|H1) / p(x|H0)当 L(x) 超过某个阈值时判 H1否则判 H0。此时阈值由虚警概率决定而虚警概率又依赖噪声方差和信号能量计算并不复杂。但实际工程里“完全已知”很少成立。更多情况是信号里还藏着几个未知参数比如直流偏置、正弦波的初相、脉冲信号的到达时刻。如果把这些参数当作随机变量去积分就是贝叶斯检测如果不知道它们的先验分布GLRT 就成了务实的选择先用最大似然估计把未知参数估出来再代入似然比。从原理上看GLRT 是给每个可能的参数值做一次匹配再取最大输出和门限比较性能虽然略逊于参数已知的最优检测但在未知参数维度不高时损失极小。参数全部已知的 GLRT 是这个框架的退化情形。这时最大似然估计退化成已知值GLRT 统计量就等价于标准似然比。但这个仿真项目特意保留了 Q.m、Qinv.m 这两个模块说明它不是简单调一次normcdf就完事而是把阈值计算、逆高斯分位数求解也拆成了独立函数方便以后往参数未知的方向扩展。2.2 参数已知情形下的似然比构造假设观测向量是 N 维的H0 下 x wH1 下 x s w其中 w 是零均值协方差矩阵为 σ²I 的高斯噪声s 是已知信号向量。两个假设下的概率密度分别为p(x|H0) (1/(2πσ²)^(N/2)) * exp(-||x||²/(2σ²)) p(x|H1) (1/(2πσ²)^(N/2)) * exp(-||x-s||²/(2σ²))两边取对数再相减得到对数似然比ln L(x) (1/(2σ²)) * (||x||² - ||x-s||²) (1/σ²) * (xᵀs - ||s||²/2)去掉与数据无关的常数项之后检测统计量就变成T(x) xᵀs也就是观测与已知信号的內积。这个结果和匹配滤波器是一致的。GLRT 在这里的“广义”体现在如果 s 里某些参数未知那么 T(x) 里会用这些参数的最大似然估计去替换从而得到T_glrt(x) max_θ xᵀs(θ)。在 MATLAB 里构造这个统计量不需要循环。假设噪声方差已知为sigma2信号向量为s观测量为x一行代码就能算出统计量T (x(:) * s(:)) / sqrt(sigma2 * (s(:) * s(:))); % 归一化相关系数形式归一化之后T 在 H0 下服从标准正态分布在 H1 下服从均值为sqrt(SNR)的正态分布。这样阈值可以直接从标准正态分位数得到也就是 Q 函数的逆。这也是项目里出现 Qinv.m 的核心原因。参数说明x(:)与s(:)都强制拉成列向量避免行向量转置错误除以信号能量开根号相当于把匹配滤波器归一化让统计量的方差恒为 1门限只和虚警率有关不再依赖信号幅度。2.3 阈值确定与 Q 函数的作用检测门限由虚警概率Pfa决定。定义 Q 函数为标准正态分布右尾概率Q(z) ∫_z^∞ (1/√(2π)) exp(-t²/2) dt如果希望虚警率不超过Pfa则门限应当满足λ Q⁻¹(Pfa)。在 MATLAB 里erfc函数和 Q 函数的关系是Q(z) 0.5 * erfc(z/√2)所以逆函数可以写成lambda sqrt(2) * erfcinv(2 * Pfa); % 由虚警率反推判决门限项目里的Qinv.m很可能就是封装了这一行变换。这样做的优势是门限独立于信号波形只要信号能量归一化做好同一套门限可以复用到任意波形检测。要注意的是如果噪声方差未知或者不是白噪声那么门限计算里还要引入噪声协方差矩阵的 Cholesky 分解做白化这是后续扩展的方向之一。下表总结了不同虚警概率下对应的门限值和等效 SNR 检测门限方便快速核对仿真参数是否设置合理PfaQinv(Pfa) 门限 λ所需 SNR(dB) 约 90% 检测率0.012.3263约 4.3 dB0.0013.0902约 6.0 dB1e-54.2649约 8.5 dB1e-75.1993约 10.4 dB这里的“所需 SNR”是在单次采样检测场景下把检测概率公式Pd Q(λ - sqrt(SNR))反推得到的。仿真时如果发现检测曲线和理论值偏差超过 0.5 dB首先要检查门限是否用了未归一化的统计量其次检查蒙特卡洛次数是否足够。3. MATLAB 仿真实现drawpd.m、Q.m 与 Qinv.m 拆解3.1 仿真框架与文件分工这套源码的文件结构非常典型适合作为检测仿真的骨架。Q.m和Qinv.m是数值函数库专门处理高斯 Q 函数及其逆运算drawpd.m负责绘制概率密度曲线用来直观展示 H0 和 H1 下统计量的分布结果.bmp是运行后的输出图一般包含检测概率曲线或 PDF 对比图。这种分工方式把一个完整的检测仿真拆成了“统计工具”和“场景脚本”后续做其他检测器时可以直接复用Q.m和Qinv.m。从工程角度看把 Q 函数单独抽出来是个好习惯。MATLAB 自带的normcdf或者erfc虽然也能用但erfc对大数值参数容易下溢而自实现时通常会做换元或对数域处理。实际项目里我一般会再包一层让Qinv支持向量输入方便一次批量计算多个门限。3.2 Q 函数及逆函数的数值实现标准正态分布的 Q 函数可以用补误差函数表达但为了数值稳定性我建议对极限情况做保护。一个实用的Q.m实现如下function y Q(x) % 标准正态分布右尾概率即 P(Z x) % 输入 x 可以是标量或向量 % 内部利用 erfc 精度高且支持向量化 x x(:); y 0.5 * erfc(x / sqrt(2)); % 对极端值做饱和处理避免向上/向下溢出产生 NaN y(x 8) 0; y(x -8) 1; endQinv.m则是对上式的逆运算输入概率值输出对应分位数function x Qinv(p) % Q 函数的逆给定右尾概率 p返回分位数 x满足 Q(x) p % 输入 p 应该在 (0,1) 区间内 p(p 0 | p 1) NaN; % 非法概率直接置 NaN便于调用方发现错误 x sqrt(2) * erfcinv(2 * p); end逻辑说明由于erfc在参数绝对值较大时可能出现数值饱和而实际仿真中门限通常不会超过 6两个函数的保护分支大多数情况下不会触发。但留着它可以避免蒙特卡洛循环里因为个别异常样本导致整个仿真崩溃。参数说明p建议取值范围在 1e-7 到 0.5 之间如果虚警率低于 1e-7建议改用对数域递推因为此时erfcinv的参数接近 2浮点分辨率可能不够。3.3 drawpd.m 中的概率密度绘制与可视化drawpd.m的核心作用是画出 H0 和 H1 两个假设下检测统计量的理论 PDF让使用者先看清楚“两个分布的重叠程度”再跑蒙特卡洛。假设我们用的是归一化内积统计量那么 H0 下 T ~ N(0,1)H1 下 T ~ N(√SNR, 1)。理论 PDF 可以直接用normpdf计算function drawpd(SNR_dB) % 绘制 H0 和 H1 下的统计量概率密度曲线 % SNR_dB信噪比单位 dB snr 10^(SNR_dB/10); t -4:0.01:8; % 统计量取值范围 pdf0 normpdf(t, 0, 1); % H0均值0方差1 pdf1 normpdf(t, sqrt(snr), 1); % H1均值 sqrt(SNR)方差1 plot(t, pdf0, b-, LineWidth, 1.5); hold on; plot(t, pdf1, r-, LineWidth, 1.5); xlabel(检测统计量 T); ylabel(概率密度); legend(H0: 纯噪声, H1: 信号噪声); title([PDF对比, SNR , num2str(SNR_dB), dB]); grid on; end这段代码的关键参数是统计量的均值。sqrt(snr)来自归一化内积的期望推导当信号能量为Es、噪声方差为σ²归一化后的 SNR 就是Es/σ²开根号后即为 H1 下统计量的均值偏移量。用-4:0.01:8作为横轴范围是因为门限通常在 2~5 之间且 H1 的均值在低 SNR 时仍在 2 附近横轴范围留够余量才能看出右尾差异。drawpd.m里的绘图语句不必和这个完全一致但核心逻辑一定是“把两个假设的分布画在同一坐标系同时标出判决门限”否则 PDF 图对参数选择的指导意义会大打折扣。3.4 主流程代码与参数设置把上面的模块串起来一个完整的参数已知 GLRT 检测仿真主流程如下%% 参数初始化 N 10; % 采样点数 Pfa 1e-3; % 期望虚警率 MC 100000; % 蒙特卡洛次数 SNR_dB 0:2:10; % 仿真信噪比范围 % 已知信号单位能量矩形脉冲 s ones(N,1) / sqrt(N); % 根据虚警率计算门限 lambda Qinv(Pfa); %% 蒙特卡洛仿真 for k 1:length(SNR_dB) snr 10^(SNR_dB(k)/10); amp sqrt(snr); % 信号幅度保证信号能量 SNR T0 zeros(MC,1); % H0 下统计量 T1 zeros(MC,1); % H1 下统计量 for m 1:MC w randn(N,1); % 标准高斯白噪声 x0 w; % H0纯噪声 x1 amp * s w; % H1信号加噪声 T0(m) x0(:) * s; % 投影到信号方向 T1(m) x1(:) * s; end Pfa_sim(k) mean(T0 lambda); Pd_sim(k) mean(T1 lambda); end逻辑说明外层循环遍历 SNR内层循环做蒙特卡洛。amp * s构造信号时由于s是单位能量向量amp的平方就是 SNR。门限lambda只依赖Pfa不随 SNR 变化这正体现了参数已知 GLRT 的恒虚警特性。mean(T0 lambda)是统计超过门限的样本比例用来验证虚警率是否落在 Pfa 附近mean(T1 lambda)就是检测概率。参数说明MC取 100000 时虚警率估计的方差约为sqrt(Pfa*(1-Pfa)/MC)对Pfa1e-3约为 0.0001足够看出和理论值的偏差。如果减小MC曲线会抖动建议不低于 20000。仿真结束后不要急着画图先打印Pfa_sim和Pfa的偏差。如果偏差超过 20%多半是门限计算用了norminv(1-Pfa)而没有注意左右尾方向或者统计量没有做能量归一化。4. 仿真结果解读与性能评估4.1 结果.bmp 中关键曲线的判读结果.bmp是这份源码运行后的输出图通常内容因不同版本而异但最常见的形式是“检测概率 Pd 随 SNR 变化”的曲线图横轴 SNR_dB纵轴检测概率同时可能叠加一条理论曲线。理论检测概率的计算方式如下Pd_theory Q(lambda - sqrt(10.^(SNR_dB/10)));这条公式来自 H1 统计量 T ~ N(√SNR, 1) 超过门限 λ 的概率。对比仿真曲线和理论曲线重点看三点低 SNR 区域Pd 小于 0.2是否贴合这反映门限计算是否正确中间区域Pd 在 0.3~0.9是否有水平偏移这反映信号构造的能量归一化是否出错高 SNR 区域是否平滑收敛到 1这里蒙特卡洛样本量不足会导致尾部抖动。如果结果.bmp里还有第二条曲线比如经验虚警率随 SNR 变化的平坦线那说明仿真同时验证了恒虚警特性。理想情况下 Pfa_sim 应该是一条水平直线幅度在 Pfa 附近随机波动。一旦发现 Pfa_sim 随 SNR 明显倾斜就要回头检查噪声产生是否用了有色噪声或者信号归一化时混入了噪声功率。4.2 检测概率与虚警概率的权衡参数已知 GLRT 没有自由调节的“灵敏度旋钮”唯一能改的就是虚警率 Pfa。虚警率每降一个数量级门限大约右移 0.5~1导致检测概率曲线整体右移。下表给出了在同一 SNR 下不同 Pfa 对 Pd 的影响数据基于 N10、SNR6 dB 的理论计算Pfa门限 λPd Q(λ - sqrt(SNR))1e-22.32630.9051e-33.09020.8401e-43.71900.7681e-54.26490.691可以看到为了把虚警从 1e-2 压到 1e-5检测概率损失约 0.2。工程上选择 Pfa 时要考虑后续处理的代价雷达一次扫描虚警太多会引起计算机饱和通信系统误同步会导致整个数据包报废。我一般习惯先定虚警率上界再反推所需 SNR用这个仿真代码先跑一遍理论曲线确认系统预算够不够而不是直接调门限碰运气。4.3 常见调试问题与坑第一个坑是统计量方向写反。x * s和s * x在实数域结果相同但如果信号里含有复数分量忘记取共轭转置就会得到错误的投影。建议统一写成x(:) * s(:)并且用(x * s)的实部作为判决量因为虚部只贡献噪声。第二个坑是门限单位不一致。Q 函数逆给出的分位数是标准正态尺度而有些实现里门限直接设为amp^2/2之类的能量阈值两者差了sqrt(SNR)倍。判断方法是把 SNR 设为 0 dB如果检测概率接近 0.5 而虚警率不等于 Pfa说明统计量或门限没有对准尺度。第三个坑和 Qinv 的数值特性有关。当Pfa小于 1e-8 时erfcinv的参数无限接近 2双精度计算会丢失有效数字。此时我一般改用对数域近似% 极小虚警率的对数域近似适用于 Pfa 1e-8 t sqrt(-2 * log(Pfa) - log(2*pi) - 2*log(lambda));这个近似在 Pfa 很小时误差可以忽略但不建议在常规范围内使用因为迭代求解更容易引入人为误差。仿真时如果发现 Pfa_sim 和设定值系统性偏差优先检查是否触发了这种数值边界。5. 把 GLRT 仿真扩展到实际应用场景5.1 参数失配时的鲁棒性验证实际工程中“参数已知”往往是理想假设。为了评估失配的影响可以改造主流程里的信号向量比如发射端信号是s_true接收端匹配滤波器用的是s_hyp两者之间存在频率偏移或多普勒失配。做法是在仿真循环里用s_true生成数据用s_hyp计算统计量s_true exp(1j*2*pi*fd*(0:N-1)); % 真实信号带多普勒频移 s_hyp ones(N,1)/sqrt(N); % 假设信号无频偏 % 生成数据用 s_true计算统计量用 s_hyp T abs(x * s_hyp); % 幅度检测相位未知时取模参数说明fd是多普勒频移与采样间隔的乘积。失配会导致输出 SNR 下降具体损耗因子是abs(s_true * s_hyp)^2 / (N^2)。你可以画一条“检测概率 vs 频偏”的曲线观察 GLRT 的性能悬崖在哪里。这个实验比单纯跑理想曲线更有工程价值因为同步误差是真实系统必然存在的。5.2 用蒙特卡洛仿真替换固定阈值标准 GLRT 的阈值由高斯假设解析给出但实际噪声可能带重尾或存在脉冲干扰。此时解析门限不再可靠更稳妥的做法是先在纯噪声条件下跑大量蒙特卡洛得到经验分布的 1-Pfa 分位数用这个分位数作为门限。实现如下MC_thresh 1000000; T_noise zeros(MC_thresh,1); for m 1:MC_thresh w randn(N,1) 0.2 * randn(N,1).^3; % 示例重尾噪声 T_noise(m) w * s; end lambda_emp quantile(T_noise, 1 - Pfa);这段代码的要点在于quantile得到的是经验门限不需要假设噪声闭式分布。代价是计算量增大但一旦噪声模型变化只需重新跑一次T_noise而检测部分的蒙特卡洛可以沿用同一个门限。注意MC_thresh要比1/Pfa至少大一个数量级否则尾部分位数估计偏差很大。5.3 从检测到估计结合 MLE 的进阶思路参数已知 GLRT 的下一步自然扩展是把信号幅度或相位作为未知量处理。做法是把检测统计量里的s替换成最大似然估计得到的s_hat。以幅度未知、波形已知为例幅度 MLE 是a_hat (x*s)/(s*s)对应的 GLRT 统计量变为T_glrt abs(x*s)^2 / (s*s); % 能量型统计量这个量在 H0 下服从指数分布H1 下服从非中心卡方分布。你可以在现有代码里加一个分支比较参数已知和幅度未知两种情况下检测概率的差异通常幅度未知会带来约 1.5~2 dB 的损失。这个扩展不需要改动Qinv.m只需把门限换成卡方分布的分位数chi2inv(1-Pfa, 1)。这样一套仿真代码就同时覆盖了参数已知与参数部分未知两类问题后续做自适应检测时可以直接复用这里的噪声生成和蒙特卡洛框架。本文还有配套的精品资源点击获取