简介: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)。两个假设写作:
- H0:x = w,只有噪声
- H1:x = 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 量级(对 M=1e5)说明实现正确。
3.3 信号检测仿真的参数表与设置原则
| 参数 | 含义 | 示例值 | 设置原则 |
|---|---|---|---|
| N | 快拍长度 | 64 | 决定信号能量 E 与频率分辨率;N 越大,同样 SNR 下 Pd 越高 |
| A | 信号幅度 | 1.0 | 与 σ² 一起决定 SNR;单独调 A 等价于调 SNR |
| f0 | 信号频率 | 20 Hz | 小于 fs/2;避开过零附近可减少离散化误差 |
| sigma2 | 噪声方差 | 0.01 | 与 A 配合得到目标 SNR_dB |
| Pfa_target | 目标虚警概率 | 1e-3 | 决定门限;蒙特卡洛次数 M 至少要大于 10/Pfa |
| M | 蒙特卡洛次数 | 1e5 | Pfa 越小,需要的 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] = 0,Var[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 仿真和理论对不上的三个常见原因
- 蒙特卡洛次数不足。Pfa=1e-4 时用 M=1e4,H0 虚警期望只有 1 次,Pfa_sim 几乎必然为 0 或某个随机小值。至少保证 M > 10/Pfa,并且看 H1 侧的 Pd 时要让 M > 20/(1−Pd)。
- 检测器用的信号与数据生成器用的信号不一致。比如数据里用了 A=1 的 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_target=1e-3,跑 M=1e5 次复信号蒙特卡洛,统计检测器输出大于 th 的频率,结果会落在 1e-3 附近。这套复信号版本的骨架,后续扩展到未知幅度或未知噪声方差的 GLRT 时,只需要替换统计量分子的估计部分,归一化分母保留原样。
本文还有配套的精品资源,点击获取