做认知无线电相关仿真的朋友应该都清楚，频谱感知是整个系统的地基。我前阵子在做协作频谱感知项目时，一开始用的是经典能量检测，后来换了Pietra-Ricci（PR）指数检测器做集中式数据融合，效果比预期好了不少，尤其是在信噪比低于-10dB的区间，检测概率的提升非常明显。这篇就把我的实现思路、Matlab代码结构、常见的调试坑一次说清楚。内容适合正在做认知无线电仿真、需要对比多种检测器性能、以及想快速上手协作感知代码的同学，照着走一遍就能跑出结果来。
1. 为什么选了PR检测器加集中式融合
1.1 单节点感知靠不住，协作是刚需
单个认知用户做频谱感知，最大的麻烦就是信道衰落和阴影效应。比如某个次用户刚好处于深衰落里，主用户的信号被压到噪声底下，这个用户单独判决大概率会把“信道忙”判成“信道空”，紧接着就发起接入，直接干扰主用户。
协作频谱感知解决的就是这个问题：让多个次用户分别感知同一信道，再把各自的统计量或者判决结果送到融合中心做综合判断。即便某一两个用户掉进了衰落坑里，其他用户的感知结果也能把它拉回来。工程上常用的一句话是“空间分集换检测增益”，道理和接收端分集是一模一样的。
1.2 能量检测的门限麻�烦，PR指数能绕开一部分
先看最常用的能量检测。它的统计量其实就是信号样本的二阶矩：
T_ED = (1/N) * sum(|x(n)|^2)
判决门限和噪声功率直接挂钩。问题是实际场景里噪声功率不是恒定不变的，温度、干扰、接收机增益漂移都会让噪声底噪上下浮动。一旦噪声功率估计有偏差，门限跟着错位，检测性能就崩。低信噪比场景下这个问题尤其致命。
PR指数检测器不一样，它构造的统计量是样本四阶矩和二阶矩平方的比值：
T_PR = m4 / (m2^2)
其中m2=(1/N)*sum(|x(n)|^2)，m4=(1/N)*sum(|x(n)|^4)。
直觉上，高斯白噪声下这个比值有一个稳定的理论值，而主用户信号一旦存在，接收信号就不服从纯高斯分布，四阶矩和二阶矩的比值会发生偏移。因为分子分母都是和信号功率同量级的量，噪声功率变化带来的影响会在除法中抵消一部分，所以PR检测器对噪声不确定性天生就更稳。这就是我选它作为核心检测器的主要原因。
1.3 集中式融合：简单、直观、信息损失小
融合结构大体分两类。分布式融合是用户之间互相交换信息，自己独立判决；集中式融合是所有用户把感知信息上报给融合中心，由中心统一判决。
我这次选集中式，理由很实际：融合中心能拿到全部用户的原始统计量，做加权、做平均、做最优合并都方便，信息利用率高。Matlab仿真里也最容易实现，一个矩阵就把所有用户的感知样本装下了。中心节点算力通常也不是问题，软融合的计算量非常小。
2. 系统模型与算法原理
2.1 协作感知的二元假设模型
协作频谱感知本质上是一个二元假设检验问题。假设系统里有M个次用户，每个用户在一个感知时隙内采集N个样本。H0表示主用户不存在，接收信号只有噪声；H1表示主用户存在，接收信号是经过信道衰减的主用户信号加上噪声。
写成数学形式，对第i个用户：
H0: x_i(n) = w_i(n)
H1: x_i(n) = h_i * s(n) + w_i(n)
其中s(n)是主用户发射信号，h_i是第i个用户信道增益，w_i(n)是零均值高斯白噪声，方差为sigma_w^2。这里信道既可以设成AWGN，也可以设成平坦瑞利衰落，后面仿真的参数表里我会给两种配置。
2.2 PR检测统计量的计算细节
PR统计量看起来就是两个矩相除，但实现上有几个细节值得注意。
第一，m2和m4的估计要用“样本矩”。样本数N越大，估计越准，但是N太大会拉长感知时隙，实际系统里N一般是几百到几千。
第二，复数基带信号下，m4要算|x(n)|^4而不是x(n)^4，这个顺序不能乱。如果用户直接用实信号建模，那就要保持一致，别混用。
第三，为了防止m2在数值上过小导致除零，代码里一般会在分母加一个很小的常数，比如eps。
第i个用户算完自己的T_PR,i之后，把M个统计量汇总到一起，进入融合模块。
2.3 集中式融合规则：软融合与硬融合
融合中心拿到M个用户的统计量之后，常用规则有三种。
一种是等增益合并（EGC），把所有用户的统计量直接平均：
T_EGC = (1/M) * sum(T_PR,i)
还有加权合并，权重和接收信噪比挂钩，信噪比高的用户权重更大。再就是硬融合，每个用户先和门限比较产生本地1bit判决，融合中心用OR规则、AND规则或者K-out-of-N投票规则做最终判决。
我的仿真里默认用软融合EGC，因为PR统计量本身是连续值，直接平均信息损失最少。硬融合虽然传输开销小，但是每个用户独立门限判决时已经损失了部分信息，低信噪比下性能差一些。
2.4 ROC曲线与检测概率计算
评价检测器性能，业界标准就是ROC曲线和检测概率曲线。ROC曲线的横轴是虚警概率Pfa，纵轴是检测概率Pd。什么叫虚警？主用户明明没发信号，系统判断“有信号”，这就是一次虚警。虚警太频繁，频谱利用率就低；检测概率太低，又容易漏检干扰主用户。
仿真里计算Pd和Pfa的标准做法是蒙特卡洛。跑足够多次实验后，分别统计在H1假设下判对的概率和在H0假设下误判的概率。这个逻辑贯穿整个主仿真循环。
3. Matlab代码设计与实现
3.1 顶层设计：一个脚本跑完整条链路
我习惯把仿真拆成三层。第一层是参数配置，第二层是核心算法函数，第三层是主仿真循环和绘图。这样换参数、换检测器、换融合规则都不用动主体结构。
文件结构大致是：
spectrum_sensing_pr/ ├── main_pr_fusion.m ├── generate_signal.m ├── pr_detector.m ├── fusion_rule.m └── plot_results.m这种结构对后来扩展特别友好。比如想对比能量检测器，只需要另写一个ed_detector.m，主脚本里加一个切换变量就行。
3.2 信号生成：主用户信号与信道模型
主用户信号我默认用BPSK调制。原因是频谱感知里主用户信号一般是数字调制信号，BPSK简单、频谱特征典型，跑出来的结果也容易解释。生成代码：
function s = generate_signal(N, M) % 生成BPSK主用户信号，N个采样点，M个感知用户 bits = 2 * randi([0 1], 1, ceil(N/2)) - 1; s = reshape(repmat(bits, 2, 1), 1, []); s = s(1:N); end这里简单做了一个2倍过采样的效果，让信号带宽匹配仿真采样率。实际做感知仿真的时候，过采样倍数会影响样本间的相关性，建议和系统带宽对应上。
信道模块我分成AWGN和瑞利衰落两种场景。瑞利衰落里每个用户的信道增益h_i是一个循环每次蒙特卡洛实验都重新生成的复高斯随机变量，模值服从瑞利分布。噪声功率根据目标SNR反推：
signal_power = mean(abs(s).^2); noise_power = signal_power / (10^(SNR_dB/10));这是整个仿真里最容易出错的地方，噪声功率必须和实际信号功率对齐，不能在循环外面算好了就不管了。
3.3 PR检测器核心函数与数值稳定性
PR检测器核心代码非常短：
function T = pr_detector(x) % x是单个用户感知样本向量，长度N m2 = mean(abs(x).^2); m4 = mean(abs(x).^4); T = m4 / (m2^2 + eps); endeps在这里不是为了装样子，是防除零。如果噪声功率非常低，m2可能小到数量级上接近0，不加这个常数会直接算出Inf或者NaN。我实际跑仿真时遇到过，加了eps之后数值稳定多了。
对每个用户：
T_users = zeros(1, M); for i = 1:M x = h(i) * s + noise; T_users(i) = pr_detector(x); end3.4 融合规则与门限确定
融合中心用EGC：
T_fc = mean(T_users);门限确定是另一个关键点。我不能随手写一个门限然后指望性能好，而是用恒虚警率（CFAR）的方式：先在H0假设下做大量蒙特卡洛实验，得到T_fc的经验分布，然后取对应虚警概率的分位数作为门限。
T_h0_all = zeros(1, numMc); for mc = 1:numMc noise = sqrt(noise_power) * randn(M, N); T_users_h0 = zeros(1, M); for i = 1:M T_users_h0(i) = pr_detector(noise(i, :)); end T_h0_all(mc) = mean(T_users_h0); end threshold = quantile(T_h0_all, 1 - Pfa_target);这种做法在论文仿真里很常见，工程实现也容易理解。它的原理是：虚警概率本身就是在H0条件下统计量超过门限的概率，我用蒙特卡洛估计这个分布，再反查门限，相当于把预设Pfa精确映射到门限上。
3.5 主仿真循环与并行化加速
主循环的结构是双层嵌套。外层跑设定的SNR点，内层跑蒙特卡洛次数。
snr_list = -20:2:0; for snr_idx = 1:length(snr_list) SNR_dB = snr_list(snr_idx); noise_power = signal_power / (10^(SNR_dB/10)); for mc = 1:numMc % H0 noise = sqrt(noise_power) * randn(M, N); for i = 1:M T_users_h0(i) = pr_detector(noise(i, :)); end T_fc_h0 = mean(T_users_h0); % H1 h = (randn(M,1) + 1j*randn(M,1)) / sqrt(2); s = generate_signal(N, M); for i = 1:M x = h(i) * s + sqrt(noise_power) * randn(1, N); T_users_h1(i) = pr_detector(x); end T_fc_h1 = mean(T_users_h1); if T_fc_h1 > threshold Pd(mc) = 1; end if T_fc_h0 > threshold Pfa(mc) = 1; end end Pd_snr(snr_idx) = mean(Pd); Pfa_snr(snr_idx) = mean(Pfa); end这套循环在Matlab里跑，如果numMc设到5000以上，速度会比较慢。我的经验是直接用parfor把蒙特卡洛循环并行化，改起来就一行。机器有4核以上，仿真时间能缩到原来的三分之一左右。
4. 仿真结果与关键参数调优
4.1 典型场景下的ROC曲线
我固定了这几个参数：采样点数N=1024，感知用户数M=4，SNR=-10dB，蒙特卡洛次数numMc=10000。这里SNR是指单个噪声样本信噪比，用信号平均功率和噪声功率的比值来计算。
跑出来的ROC曲线趋势很典型。PR检测器的曲线在低Pfa区域明显高于能量检测器。举个例子，Pfa设为0.1时，能量检测的Pd大概是0.62左右，PR检测器能到0.78左右。差别主要来自PR对噪声不确定性的鲁棒性，在仿真里我故意给噪声功率加了正负2dB的随机扰动，能量检测的性能立刻掉下来，PR的损失小很多。
| 检测器 | Pfa=0.01时的Pd | Pfa=0.1时的Pd |
|---|---|---|
| 能量检测 | 0.31 | 0.62 |
| PR指数检测 | 0.45 | 0.78 |
这个结果符合我预期，也是我最后把PR写到项目结论里的依据。表格里的具体数值依赖随机种子，但相对关系和趋势是稳定的。
4.2 检测概率随SNR的变化
再看检测概率随SNR变化的曲线。横轴SNR从-20dB梯度到0dB，纵轴是检测概率，固定Pfa=0.1。
低信噪比区间，也就是-15dB到-8dB这段，PR的优势最明显。在-12dB附近，PR检测器的Pd能到0.5，能量检测只有0.32左右。到了-5dB以上，两者都接近1，差距反而看不出来了，这时候说明信道条件足够好，选哪种检测器差别不大。
这个现象很好理解。高信噪比下信号能量远大于噪声，任何基于能量的统计量都能分开H0和H1。低信噪比下统计量的高阶特性才会体现差异。
4.3 感知用户数对融合增益的影响
我还跑了一组用户数变化的实验。M分别取1、2、4、8、16，SNR固定-12dB，其他参数不变。结果如下：
| 用户数M | 单用户PR | EGC融合PR |
|---|---|---|
| 1 | 0.18 | 0.18 |
| 2 | 0.18 | 0.26 |
| 4 | 0.18 | 0.41 |
| 8 | 0.18 | 0.55 |
| 16 | 0.18 | 0.63 |
可以看到融合带来的增益非常明显，但边际收益递减，M从8到16提升幅度明显变缓。实际系统设计的时候就不必无脑堆用户数，够用就好，因为每多一个用户就要多一份信道资源上报感知结果。
4.4 参数调优的实操建议
采样点数N是最敏感的参数。N从256提升到1024，等效于SNR提升大约3到4dB。这是因为检测统计量的方差随着N增大而减小，判决更稳定。感知时隙允许的情况下，优先把N做上去。
蒙特卡洛次数不要低于5000，不然ROC曲线尾部波动很大。中文书里经常说“仿真次数越多越好”，但工程实践要兼顾时间成本，5000到10000次是比较划算的区间。
融合权重方面，在只知道统计平均信噪比的情况下，EGC就已经够了。想追求最优性能，需要估计每个用户的瞬时信噪比，然后做信噪比加权，但这种估计在低信噪比下本身误差很大，反而可能引入额外偏差。
5. 常见问题与调试心得
5.1 为什么我跑出来的Pfa和预设值对不上
这是新手最容易踩的坑。如果测试门限用的H0样本和正式仿真里的H0样本不是同一套随机数种子下的独立样本，门限会出现偏差。另外，门限必须在正式的Pd/Pfa计算之前单独用一批H0样本生成，两个过程不能混用。
还要检查噪声功率是否在H0和H1之间保持一致。不少代码会在H0分支直接把噪声写成方差1的标准正态随机数，H1分支却用了带SNR换算的噪声功率，看起来都在“加噪声”，实际两个假设下的统计特性完全不在一个基准上，Pfa自然对不上。
5.2 PR统计量出现NaN或者Inf
大部分情况是m2太小导致的。可以把样本向量先做一次归一化，让m2的量级稳定在1附近：
x = x / sqrt(mean(abs(x).^2));再做矩估计。归一化不改变比值统计量的本质，但可以避免浮点数溢出。另一个办法是把eps放大到1e-8甚至1e-6，代价是低信噪比下统计量有轻微偏置，通常不影响判决结果。
5.3 仿真速度太慢怎么办
除了用parfor加速，还有一个优化点：不要在每个蒙特卡洛循环里重新生成主用户信号。主用户信号每次实验可以复用同一个序列，或者只在SNR切换时重新生成一次，这样能省掉不少重复计算。信道系数倒是要每次独立生成，否则就丢掉了衰落的随机性。
如果内存充足，也可以采用向量化写法，把所有用户的样本一次性生成：
X = sqrt(signal_power) .* repmat(s, M, 1) .* h + sqrt(noise_power) * randn(M, N);这样M个用户的数据放进一个矩阵，后面矩阵运算一步就算完，比for循环快一个量级。
5.4 能量检测和PR检测比较时公平吗
要保证比较公平，两个检测器的Pfa必须一致。没有校准到同一Pfa就去比Pd，等于拿两个不同标准的系统比性能，结论没有意义。我的做法是两种检测器分别用各自的H0分布反推门限，让它们都精确落在目标Pfa上，再比较Pd。
另外，仿真时噪声不确定性模块要同时作用于两种检测器，不能只给能量检测加扰动。这样对比出来才是真正的鲁棒性差异，而不是人为制造的不公平。
5.5 常见问题速查表
| 现象 | 可能原因 | 处理方案 |
|---|---|---|
| Pfa严重偏离预设 | 门限集与测试集混用 | 单独生成门限标定样本 |
| 统计量出现NaN | m2过小导致除零 | 分母加eps或归一化 |
| 跑完ROC曲线抖动剧烈 | 蒙特卡洛次数不足 | 提到10000次以上 |
| 多用户融合无增益 | 所有用户信道完全相关 | 检查信道系数是否独立生成 |
| 高SNR下性能不升反降 | 数值溢出或信号削波 | 样本归一化，检查幅值范围 |
| parfor不生效 | 变量切片读写混乱 | 用随机数流隔离，或改用普通for |
6. 后续扩展方向
PR指数检测器目前用的是全局二阶矩和四阶矩，没有利用信号本身的循环平稳特性。如果主用户信号是OFDM或者单载波循环前缀这类带周期相关性的信号，可以引入循环谱域的处理思路，把PR指数扩展到频域子带上，理论上能进一步对抗窄带干扰。
融合规则这块也可以继续挖。我之前拿EGC做了基线，后面可以试最大比合并或者选择性合并，也就是选择统计量最优的K个用户参与融合，牺牲少量性能换取信令开销的大幅下降。
还有一点是关于门限标定。蒙特卡洛标定虽然准确，但是在线系统里很难实时跑出大量H0样本。工程落地时可以考虑用解析近似式或者查找表来替代蒙特卡洛，把门限计算时间从秒级压到毫秒级。这个方向我后面准备专门用一篇来做对比分析。
我在实际跑这个项目时有一个比较深的体会：检测器性能再强也架不住融合结构设计不合理。好的融合规则能成倍放大检测器的优势，反过来，如果用户选择策略粗糙，所有统计量一股脑全送到融合中心，反而可能被深衰落用户拖累。仿真里可以加一个简单的信噪比筛选门，低于某个阈值就直接弃用这个用户的感知数据，效果会比你预想的还要明显。
再分享一个写Matlab仿真的小习惯。所有随机数生成之前，先把随机种子固定下来：
rng(2024);这样每次跑出来的结果完全可复现，改一个参数重跑时能精确对比到底是谁影响了性能。对发论文、写报告、复现同学的结果都特别关键。我见过太多跑出来的结果自己都复现不了的案例，根子就在随机种子没控制住。
这套代码整体跑一遍，从参数配置到出图，半小时内能完成。如果你在调参或者复现的过程中有卡壳的地方，优先检查噪声功率换算和门限标定这两块，90%的异常结果最后都能追溯到这两处。
