
简介MATLAB参数已知条件下的GLRT信号检测仿真面向通信、雷达等需要掌握广义似然比检验原理的研究生、工程师及高校相关专业学生。围绕随机信号检测中零假设与备择假设的判定问题资源以完整工程代码呈现似然函数建模、对数似然比计算、参数估计与阈值比较等核心步骤并专门针对检测门限设置给出Q函数与逆Q函数的实现便于控制虚警概率。包内共4个文件以3个M文件为主分别承担仿真主流程、概率密度绘制与Q函数求解等任务辅以1张BMP结果图直观展示检测性能。整个压缩包仅7KB结构清晰、轻量易读无需大型工程环境即可快速运行。已有1107人学习下载。通过研读源码和结果图可举一反三迁移到其他参数已知的检测场景既能深化GLRT理论理解也能巩固MATLAB信号处理实操能力。1. 参数都已知的 GLRT 检测先把问题锁定在二元判决上接收机每天都在回答一个问题这一帧观测里目标信号到底在不在。雷达回波判读、通信帧头检测、频谱感知本质上是同一个二元判决问题。GLRT广义似然比检验是这类问题的通用框架而“参数都已知”是最容易被跳过、但又最值得先跑通的一种情形——此时 GLRT 直接退化为 Neyman-Pearson 检测器判决统计量就是匹配滤波器的输出。这个仿真标题真正想让你掌握的是在 MATLAB 里把虚警概率、检测概率、门限三者一次跑成对应关系并能为将来释放某个未知参数留好代码接口。它适合正在做信号检测仿真、不想只在公式里推演结果的工程师和研究生。2. 参数全已知时 GLRT 退化成什么从似然比到匹配滤波2.1 二元假设与似然比GLRT 的骨架设观测向量 x ∈ Rᴺ噪声 w ~ N(0, σ²I)。两个假设写作H0x w只有噪声H1x s w信号 s 的幅度、频率、相位全部已知似然比检验的判决规则是 Λ(x) p(x|H1)/p(x|H0) γ。GLRT 的通用写法里H1 的参数如果未知要先对参数求极大似然估计再代入似然比。但标题给出的前提是“参数都已知”所以 s 就是一个确定的向量似然比里没有任何需要最大化的未知量。对高斯噪声把似然函数展开并取对数可以得到一个非常干净的结果ln Λ (sᵀx)/σ² − E/(2σ²)其中 E sᵀs这个式子说明判决完全等价于把 sᵀx 和某个门限比较。sᵀx 就是匹配滤波器的输出而除以 σ² 只是一个不影响判决逻辑的定标因子。换句话说参数全部已知时GLRT 不再需要“广义”的那一步最大化剩下的就是一个固定系数的线性检测器。2.2 参数全已知时GLRT 就是 NP 检测器Neyman-Pearson 准则做的事情是在给定虚警概率 Pfa 的前提下最大化检测概率 Pd。当 H0 和 H1 的分布完全已知似然比检验就是达到这个最优的判决方式。所以参数已知的 GLRT 与 NP 检测器是同一件事只是名字不同。工程上保留“GLRT”这个叫法是有实际意义的仿真代码的骨架按 GLRT 写之后如果想把幅度 A 改成未知参数只需要把固定的 s 换成 A·s再对 A 做一次一维搜索或闭式最大化主循环和门限标定逻辑都不用动。很多教材把这种退化情形仍然叫 GLRT是因为框架一致而不是真的做了广义最大化。初学者常犯的错误是在 A 明明已知的时候仍然对幅度做 max得到一个和理论分布对不上的统计量导致 Pfa 标定失败。判断标准很简单如果你的统计量里出现了估计值那它就不再是本节这个“参数已知”的检测器。2.3 把 GLRT 统计量写成一个 MATLAB 函数function t glrt_stat_known(x, s, sigma2) % 参数全部已知时的GLRT检测统计量 % x : N x 1 观测向量 % s : N x 1 已知信号能量 E s*s % sigma2: 已知噪声方差 t (s * x) / sigma2; % 实数高斯噪声下的充分统计量 ends*x是匹配滤波输出/ sigma2不影响判决结果但保留它能让后面的理论门限公式与概率分布直接对应。注意这里假设噪声是白高斯且方差已知如果噪声是色噪声统计量要换成 sᵀR⁻¹x即先做预白化再匹配参数已知的 GLRT 框架会自然地给出这个形式这也是为什么把检测器写成独立函数而不是嵌在仿真循环里更利于扩展。3. 用 MATLAB 搭一个参数已知的 GLRT 信号检测仿真3.1 信号模型把“参数都已知”落实到具体数值仿真里“已知”的含义是信号幅度 A、频率 f0、初始相位 φ、噪声方差 σ² 在生成数据时给定检测器也使用同一组真值。这样做的意义是能够验证理论门限和检测概率真实系统里的参数来自标定或估计那属于未知参数 GLRT 的范畴不在本节讨论。fs 256; % 采样率 N 64; % 快拍长度 t (0:N-1) / fs; A 1.0; % 幅度已知 f0 20; % 频率已知 phi 0; % 相位已知 s A * cos(2*pi*f0*t phi); % 已知信号 E s * s; % 信号能量 sigma2 0.01; % 噪声方差已知 SNR_dB 10*log10(E / sigma2); % 检测信噪比参数之间的换算关系值得注意A 和 σ² 同时放大或缩小相同倍数时SNR_dB 不变检测概率理论上也不变这是后面做自检的一个好抓手。f0 的选取要避开采样率的一半否则余弦信号在离散采样下会产生混叠仿真结果会和理论曲线系统性偏离。3.2 蒙特卡洛仿真主循环虚警概率与检测概率分开统计信号检测仿真的核心是蒙特卡洛试验重复生成随机噪声、做判决、统计频率。虚警概率要在 H0 条件下统计检测概率要在 H1 条件下统计两者必须使用独立的数据生成过程不能混用同一批随机数。Pfa_target 1e-3; % 目标虚警概率 th sqrt(E / sigma2) * norminv(1 - Pfa_target); % 理论门限 M 1e5; % 蒙特卡洛次数 det_h0 false(M, 1); % H0 下的判决结果 det_h1 false(M, 1); % H1 下的判决结果 for k 1:M w sqrt(sigma2) * randn(N, 1); % 生成噪声 x0 w; % H0只有噪声 x1 s w; % H1信号加噪声 det_h0(k) glrt_stat_known(x0, s, sigma2) th; det_h1(k) glrt_stat_known(x1, s, sigma2) th; end Pfa_sim mean(det_h0); % 仿真虚警概率 Pd_sim mean(det_h1); % 仿真检测概率门限 th 来自下一章要详细推导的理论分布这里先直接使用。norminv(1 - Pfa_target)是标准正态分布的 1−Pfa 分位数。循环里的两个判决互相独立不共享随机数所以 Pfa_sim 和 Pd_sim 的统计误差不会相互传染。仿真结束后把 Pfa_sim、Pd_sim 与目标值、理论值对照偏差在 1e-5 量级对 M1e5说明实现正确。3.3 信号检测仿真的参数表与设置原则参数含义示例值设置原则N快拍长度64决定信号能量 E 与频率分辨率N 越大同样 SNR 下 Pd 越高A信号幅度1.0与 σ² 一起决定 SNR单独调 A 等价于调 SNRf0信号频率20 Hz小于 fs/2避开过零附近可减少离散化误差sigma2噪声方差0.01与 A 配合得到目标 SNR_dBPfa_target目标虚警概率1e-3决定门限蒙特卡洛次数 M 至少要大于 10/PfaM蒙特卡洛次数1e5Pfa 越小需要的 M 越大否则虚警点太少统计抖动剧烈提示Pfa 取 1e-3 时M 至少取 1e5这样 H0 下大约能统计到 100 次虚警Pfa_sim 才稳定。只跑 1e4 次的话Pfa_sim 可能落在 0.5e-3 到 1.5e-3 之间很容易误判代码有 bug。4. 门限推导与 ROC 曲线用 MATLAB 画图验证仿真没写错4.1 理论门限与检测概率的闭式解参数已知的检测统计量 t sᵀx/σ² 在 H0 和 H1 下都服从高斯分布这是它能用闭式解验证的根本原因。推导如下H0 下E[t] 0Var[t] E/σ²即 t ~ N(0, E/σ²)H1 下E[t] E/σ²Var[t] E/σ²均值等于方差记 d sqrt(E/σ²)则 H0 下 t ~ N(0, d²)H1 下 t ~ N(d², d²)。门限由虚警概率反推th d · norminv(1 − Pfa)检测概率的闭式解是Pd 1 − normcdf((th − d²)/d) normcdf(d − norminv(1 − Pfa))如果用 Communications Toolbox 的 qfunc 表示等价于 Pd qfunc(qfuncinv(Pfa) − d)。两个形式计算结果一致选一个用即可。这个式子也说明参数已知 GLRT 的检测概率只由 d sqrt(E/σ²) 决定和信号的具体波形无关——波形只通过能量进入性能。4.2 用 MATLAB 画图工具验证 ROC 曲线ROC 曲线描述的是 Pd 随 Pfa 的变化关系是信号检测仿真最常用的验证手段。下面的代码把理论 ROC 和蒙特卡洛点画在同一张图上Pfa_vec logspace(-4, 0, 21); % 21 个 Pfa 刻度点 Pd_mc zeros(size(Pfa_vec)); M 2e4; % 每个点单独跑蒙特卡洛 for i 1:numel(Pfa_vec) th sqrt(E/sigma2) * norminv(1 - Pfa_vec(i)); hits 0; for k 1:M x s sqrt(sigma2)*randn(N,1); % 只统计 H1 if (s*x)/sigma2 th hits hits 1; end end Pd_mc(i) hits / M; end d sqrt(E / sigma2); Pd_theory normcdf(d - norminv(1 - Pfa_vec)); semilogx(Pfa_vec, Pd_theory, -, LineWidth, 1.5); hold on; semilogx(Pfa_vec, Pd_mc, o, MarkerSize, 5); grid on; xlabel(P_{fa}); ylabel(P_d); legend(理论 ROC, 蒙特卡洛, Location, southeast); title(参数已知 GLRT 的 ROC 曲线);这段代码把每个 Pfa 刻度点上的蒙特卡洛试验拆开跑理论曲线用 semilogx 画在对数横坐标上能同时看到低虚警区和高虚警区的行为。理论线和散点之间的偏差来源有两条一是 M 有限导致的统计抖动Pfa 越小抖动越大二是 Pfa_vec 最右端接近 1 时logspace 取点过密曲线尾部对数值敏感。4.3 仿真和理论对不上的三个常见原因蒙特卡洛次数不足。Pfa1e-4 时用 M1e4H0 虚警期望只有 1 次Pfa_sim 几乎必然为 0 或某个随机小值。至少保证 M 10/Pfa并且看 H1 侧的 Pd 时要让 M 20/(1−Pd)。检测器用的信号与数据生成器用的信号不一致。比如数据里用了 A1 的 s检测器里却用了 s/norm(s)门限公式里的 E 就对不上。把 s 和 E 作为同一个变量的派生量能避免这类问题。方差与标准差混用。sigma2 是方差代码里生成噪声要写 sqrt(sigma2)*randn门限里用的是 sqrt(E/sigma2)。任何一处写成直接乘 sigma2统计量分布整体偏移Pfa 会偏离目标一个数量级。5. 从实数到复信号参数已知 GLRT 的 I/Q 扩展5.1 复基带信号的统计量I/Q 两路合并实际接收机做数字下变频后观测是复基带信号I/Q 两路都携带噪声。模型变为 x_c s_c w_c其中 w_c 的实部和虚部独立各服从 N(0, σ²/2)这样每个复采样点总方差才是 σ²。对复高斯分布重推似然比得到的充分统计量是T Re(s_cᴴ x_c) · sqrt(2 / (σ²E))这个统计量在 H0 下服从标准正态分布在 H1 下均值变成 sqrt(2E/σ²)。也就是说复信号情形下的检测增益是 d² 2E/σ²比同功率实数信号高一倍对应约 3 dB 的增益来源是 I/Q 两路噪声在匹配滤波时被平均掉了。5.2 复信号门限修正与一个容易踩的坑s_c A * exp(1j*2*pi*f0*t); % 复基带已知信号 E_c real(s_c * s_c); % 复信号能量实数标量 w_c sqrt(sigma2/2) * (randn(N,1) 1j*randn(N,1)); % 复噪声 x_c s_c w_c; T real(s_c * x_c) * sqrt(2 / (sigma2*E_c)); % 归一化统计量 th norminv(1 - Pfa_target); % H0 下 T~N(0,1)门限就是分位数 detected T th;最容易踩的坑是复噪声生成时把 sqrt(sigma2/2) 写成 sqrt(sigma2)。这样实部方差变成 σ²T 在 H0 下的分布不再是标准正态门限按 norminv 算出来的 Pfa 会翻倍。反过来如果统计量忘了乘 sqrt(2/(sigma2*E_c)) 的归一化因子门限就必须改回带 E 和 σ² 的形式两种写法只能选一种。验证方法是固定 Pfa_target1e-3跑 M1e5 次复信号蒙特卡洛统计检测器输出大于 th 的频率结果会落在 1e-3 附近。这套复信号版本的骨架后续扩展到未知幅度或未知噪声方差的 GLRT 时只需要替换统计量分子的估计部分归一化分母保留原样。本文还有配套的精品资源点击获取