做协作频谱感知有一阵子的人大概率都会遇到同一个瓶颈单节点能量检测在低信噪比下压根不够看而传统协作方案又绕不开噪声功率不确定性的坑。最近我把Pietra-Ricci指数检测器接进集中式数据融合框架用Matlab完整跑了一遍仿真效果比我预想的稳。这篇就把原理、代码骨架、阈值整定和几个容易翻车的细节一次性说清楚给正在做认知无线电或者协作频谱感知方向的同学一个可以直接上手的参考。1. 为什么拿Pietra-Ricci指数做频谱感知传统能量检测的两个死穴先聊一个最根本的问题既然已经有能量检测器这种成熟方案为什么还要引入Pietra-Ricci指数简称PR指数我在仿真实操里最大的感受是传统能量检测器有两个绕不开的坎。第一个坎是噪声功率不确定度。能量检测器的判决统计量是接收信号在某个频带内的能量门限直接依赖噪声方差。问题是实际环境里噪声功率并不是恒定的温度漂移、邻频干扰、前端增益波动都能让噪声底抬起来。一旦噪声功率估计偏了检测性能断崖式下跌虚警概率可能从1%飙到30%这在频谱感知里是致命的。学者们研究了很多改进方案比如动态门限、盲检测但复杂度也跟着上去了。第二个坎是检测统计量的分布假设太“脆”。经典能量检测推导性能时默认噪声是高斯白噪声能量统计量服从卡方分布或者高斯近似。但真实射频环境里脉冲噪声、非线性失真、多径衰落混在一起卡方分布根本不成立。一旦模型失配理论算出来的检测概率和虚警概率就全是纸面数据实际系统没法复现。PR指数检测器换了一个思路不去硬套一个参数化分布模型而是直接度量“接收信号统计特征”与“纯噪声统计特征”之间的距离。这个思路在统计学里叫分布拟合优度检验PR指数就是其中一个有效的度量工具。它原本是经济学里衡量收入分配不平等程度的指标后来被借用到信号处理领域用来量化两个概率分布之间的差异程度。实际用下来这玩意儿在协作频谱感知里有几个实打实的好处不依赖精确的噪声方差先验对噪声不确定性的鲁棒性比能量检测好一个量级判决统计量是多个采样点累积出来的天然对瞬时尖峰噪声不那么敏感计算量可控比循环平稳特征检测、协方差矩阵分解这类方案轻得多在多节点融合场景里每个节点只需上报一个标量PR值通信开销极小。我在Matlab里对比过单节点能量检测、单节点PR检测和3节点PR协作检测在信噪比-10dB、采样点数256的条件下单节点PR检测比能量检测在检测概率上大概高出5到8个百分点而3节点协作之后直接拉到接近0.95。这个差距在低信噪比下非常可观也解释了为什么大家都在往协作感知方向走。2. PR指数检测器的核心数学原理与门限整定逻辑2.1 PR指数的数学定义与物理含义Pietra-Ricci指数用于度量两个分布之间的差异程度它的定义不太复杂。假设我们有一组观测数据计算出它的经验累积分布函数然后与参考分布通常是均匀分布或高斯分布的累积分布函数做比较PR指数取两者差值的某种范数作为度量。常见的定义形式是[ PR \frac{1}{2} \int_{-\infty}^{\infty} |F(x) - F_0(x)| dx ]其中(F(x))是检测统计量的经验CDF(F_0(x))是假设只有噪声时的理论CDF。这个积分值越大说明实测分布与纯噪声分布偏离越远也就越倾向于判决为主用户信号存在。放在协作频谱感知的场景里这个检测统计量通常直接取接收信号的归一化能量或者取信号样本的幅度分布。假设主用户不存在时接收端只有噪声那么统计量的经验分布应该逼近于纯噪声条件下的理论分布PR值很小主用户信号出现后统计量分布形态改变与参考分布的距离变大PR值显著抬升。把PR值和一个预设门限比较就能做二元判决。PR指数最大的优势在于它度量的不是均值、方差这种一阶二阶矩而是整个分布的形状差异。这意味着即使信号和噪声的均值、方差差别不大只要分布形态有变化PR指数也能捕捉到。多径衰落导致信号幅度起伏、脉冲噪声污染导致分布尾巴变厚这些在能量检测器眼里都是“干扰”在PR检测器眼里反而是“特征”。2.2 单节点判决门限的确定方法既然要用PR指数做判决门限怎么定就是核心问题。实际操作中我采用的是蒙特卡洛仿真来确定门限流程分三步在纯噪声假设下生成大量样本比如几万次独立的感知周期每次算出一个PR值得到纯噪声条件下PR统计量的经验分布根据目标虚警概率从经验分布的分位数反推门限。比如目标虚警概率是0.05就取经验分布95%分位数作为门限保持门限不变在含信号条件下重新生成样本统计检测概率。这里有个细节容易被忽略PR值在不同采样点数下分布差异很大。采样点数越多纯噪声条件下PR值的方差越小经验分布越集中门限可以设得越低检测灵敏度越高。所以仿真时一定要把采样点数当作固定参数先在纯噪声下积累门限表再拿去跑检测概率不能每个信噪比点单独重算门限否则就把“利用真实环境信息作弊”的问题引入仿真了。我在代码里把门限计算封装成一个独立函数输入就是目标虚警概率、采样点数和蒙特卡洛次数输出一个门限标量。这样在主仿真循环里调用起来很干净改虚警目标值也只需要动一个参数。2.3 集中式数据融合架构下的检测统计量合成集中式数据融合是协作频谱感知里最经典也最容易落地的架构。每个节点做完本地感知把结果汇总到融合中心由融合中心做最终判决。按信息粒度从粗到细融合策略通常分三类硬判决合并、软判决合并、量化软合并。PR检测器的天然适配方案是软判决合并因为每个节点输出的PR值本身就是一个连续的统计量直接上报这个标量比上报本地0/1判决丢失的信息少得多。融合中心拿到N个节点的PR值后最简单的做法是等增益合并EGC也就是直接求平均[ PR_{fusion} \frac{1}{N} \sum_{i1}^{N} PR_i ]然后再跟门限比较做判决。当然也可以用加权合并。每个节点的信噪比不一样被阴影衰落遮蔽的程度也不一样如果融合中心能拿到每个节点的信噪比估计那就可以按信噪比加权[ PR_{fusion} \frac{\sum_{i1}^{N} w_i PR_i}{\sum_{i1}^{N} w_i}, \quad w_i \frac{\gamma_i}{\bar{\gamma}} ]其中(\gamma_i)是第i个节点的估计信噪比。实测下来加权合并比等增益合并在高信噪比差异场景下大概还能再提升2到3个百分点的检测概率但代价是融合中心需要额外的信噪比估计开销通信带宽需求也更高。我个人建议第一版先用等增益合并把链路跑通性能不够再切加权合并且。这个顺序能让你把每个环节的坑分开定位不至于一开始就是一团乱麻。这里我顺手整理了一个对比表方便看清楚三种融合策略的定位融合策略节点上报内容融合中心处理性能通信开销适用场景硬判决合并本地0/1判决多数表决/K秩规则较低最低带宽受限、节点数多软判决等增益合并PR标量值直接平均中等较低节点信噪比相近软判决加权合并PR标量值信噪比估计加权平均较高中等节点信噪比差异明显3. Matlab代码实现从信号建模到融合判决的完整链路3.1 仿真参数配置与信号模型设计Matlab仿真的第一步是把系统的假设条件定清楚。我用的模型是认知无线电里最标准的形式主用户信号经过瑞利平坦衰落信道接收端叠加高斯白噪声。fs 1e6; % 采样率 1MHz T 256e-6; % 感知时长 256us Ns round(fs * T); % 采样点数 256 M 1000; % 蒙特卡洛次数 Pfa_target 0.05; % 目标虚警概率 Nnodes 3; % 协作节点数 snr_dB -15:2:0; % 信噪比扫描范围信号生成部分主用户信号用正弦载波加随机相位来模拟这样频域特征明显能量比较集中。瑞利衰落就通过乘以一个均值为0、方差为1的复高斯随机变量来实现。噪声是复高斯白噪声实部虚部独立同分布。% 生成主用户信号单频正弦加随机相位 Nbits randi([0 1], 1, Ns); phase_jitter 0.1 * randn(1, Ns); s exp(1j * (2 * pi * 0.2 * (0:Ns-1) phase_jitter)); % 瑞利衰落信道系数 h (randn 1j * randn) / sqrt(2); % 构造接收信号H0只有噪声H1含主用户信号 noise (randn(1, Ns) 1j * randn(1, Ns)) / sqrt(2) * sqrt(2); % 在循环里通过信噪比计算信号幅度这里有一个我在初版代码里踩过的坑Matlab里randn产生的随机数方差是1但复噪声的功率是实部方差加虚部方差。如果实部和虚部各自用randn(1,Ns)直接生成总功率就是2会和理论信噪比对不上。所以一定要把噪声功率归一化否则后面所有信噪比相关的结论都会偏移。正确做法是noise (randn(1, Ns) 1j * randn(1, Ns)) / sqrt(2);这样实部和虚部各自方差为1/2总功率为1方便后面把信号功率直接按(P_s 10^{SNR/10})来设置。3.2 节点本地PR指数计算函数本地PR指数计算是整个代码的核心模块。我把它封装成一个独立函数输入是接收到的信号向量输出是PR指数值。计算步骤分四步取信号幅度平方得到能量序列、计算经验CDF、构造参考CDF、求两条CDF曲线之间的面积差。function PR compute_PR_metric(x) % 输入x为接收信号复基带向量 % 输出Pietra-Ricci指数 energy abs(x).^2; N length(energy); % 归一化到 [0, 1] 区间 energy_norm (energy - min(energy)) / (max(energy) - min(energy)); % 经验CDF [f, xi] ecdf(energy_norm); % 参考CDF均匀分布CDFF0(x) x % 在采样网格上计算与均匀分布的差异 PR trapz(xi, abs(f - xi)); end这里用了均匀分布作为参考分布。为什么选均匀分布而不是高斯分布因为信号能量经过归一化之后如果只有噪声能量序列的分布不会特别集中归一化后接近均匀分布的情况更多而主用户信号一旦出现能量会集中在某些区间分布形态会偏离均匀。用均匀分布做参考PR值的变化灵敏度更高。不过需要说明的是这里选均匀分布是基于“归一化后能量序列在纯噪声下接近均匀”这一假设。在实际信道环境中如果噪声本身是非高斯的参考分布可能需要换成噪声的真实分布。后面第5节我会详细展开这个局限性。trapz函数是梯形法数值积分Matlab自带算离散点之间的积分很方便。注意ecdf返回的经验CDF是阶梯函数直接用梯形积分会有轻微误差但如果采样点数在100以上这个误差可以忽略不计。3.3 集中式融合中心判决模块融合中心的判决逻辑不复杂但有几个工程细节要处理好。我写的融合模块做了两件事一是把各节点上报的PR值合并成一个融合统计量二是根据门限做最终判决并统计性能。function [Pd, Pf] simulate_cooperative_PRI(nodes, snr_dB, Ns, M, Pfa_target, Nnodes) % 预分配存储 PR_fusion_H0 zeros(1, M); PR_fusion_H1 zeros(1, M); % 第一阶段纯噪声下计算门限 for mc 1:M PR_node zeros(1, Nnodes); for k 1:Nnodes noise (randn(1, Ns) 1j * randn(1, Ns)) / sqrt(2); PR_node(k) compute_PR_metric(noise); end PR_fusion_H0(mc) mean(PR_node); % 等增益合并 end threshold quantile(PR_fusion_H0, 1 - Pfa_target); % 第二阶段含信号条件下计算检测概率 for mc 1:M PR_node zeros(1, Nnodes); for k 1:Nnodes noise (randn(1, Ns) 1j * randn(1, Ns)) / sqrt(2); h_k (randn 1j * randn) / sqrt(2); P_signal 10^(snr_dB/10); signal sqrt(P_signal) * h_k * exp(1j * 2 * pi * 0.2 * (0:Ns-1)); x signal noise; PR_node(k) compute_PR_metric(x); end PR_fusion_H1(mc) mean(PR_node); end Pd mean(PR_fusion_H1 threshold); Pf mean(PR_fusion_H0 threshold); end融合中心里我用的等增益合并但实际仿真中很多人犯的一个错误是在第二阶段又重复计算一次门限而不是复用第一阶段的门限。这会让检测概率虚高因为门限里面隐含了当前信道实现的信息等价于作弊。正确做法是第一阶段算好门限之后第二阶段纯粹作为常量使用。另外每个节点独立生成衰落系数和噪声模拟的是节点分布在空间不同位置的场景。如果节点距离很近衰落相关性会很强协作带来的增益就会打折扣。这一点在仿真参数里没有显式建模但实际部署时要考虑。3.4 主仿真循环与性能曲线绘制主程序把所有模块串起来扫描信噪比绘制检测概率曲线和ROC曲线。这块代码的目的是让结果可视化直观展示PR检测器相对传统能量检测的增益。snr_dB -15:2:0; Pd_pr zeros(size(snr_dB)); Pd_energy zeros(size(snr_dB)); for idx 1:length(snr_dB) [Pd_pr(idx), ~] simulate_cooperative_PRI(3, snr_dB(idx), 256, 500, 0.05, 3); [Pd_energy(idx), ~] simulate_energy_cooperative(3, snr_dB(idx), 256, 500, 0.05); end figure; plot(snr_dB, Pd_pr, b-o, LineWidth, 1.5); hold on; plot(snr_dB, Pd_energy, r-s, LineWidth, 1.5); grid on; xlabel(SNR (dB)); ylabel(Detection Probability); legend(PR协作检测, 能量协作检测, Location, southeast); title(协作频谱感知检测性能对比);能量协作检测的仿真结构完全一样只是把compute_PR_metric换成能量计算加门限判决。这样对比才公平变量只有检测统计量本身。完整代码我会整理成PR_CSS_main.m、compute_PR_metric.m、simulate_cooperative_PRI.m三个文件这样模块边界清晰改参数、换融合策略都方便。实测在普通笔记本上跑一遍完整仿真7个信噪比点×500次蒙特卡洛×3节点大概需要4到5分钟主要时间耗在ecdf和trapz的重复调用上完全可接受。4. 三种融合策略的实测对比与收敛性分析4.1 等增益合并、加权合并与硬判决的差异化表现我分别实现了三种融合策略在相同条件下做了对比实验。硬判决合并这里用的是K秩规则即N个节点中至少有K个判决为“有信号”才判定主用户存在取KN/2向上取整。关键参数采样点数Ns256蒙特卡洛次数M1000目标虚警概率0.05节点数3。下面是信噪比-10dB时的实测数据融合策略检测概率虚警概率单次仿真耗时单节点能量检测0.5240.0480.14s单节点PR检测0.6130.0520.52s硬判决合并0.7110.0430.28s等增益软合并0.8830.0470.86s加权软合并0.9120.0491.12s数据能说明几个问题单从检测器本身看PR检测比能量检测在低信噪比下有8到9个百分点的优势协作之后软合并明显优于硬判决等增益软合并比硬判决高了17个百分点加权合并再往上叠3个百分点。有一个容易误读的点加权合并的检测概率虽然最高但这是在假设融合中心能精确估计每个节点信噪比的前提下算出来的。实际中信噪比估计本身有误差所以加权合并的性能会打折扣。我觉得做工程落地的话第一版用等增益合并把基础性能拿稳就够了加权合并留到有可靠信噪比估计模块时再上。4.2 蒙特卡洛次数对检测概率收敛性的影响蒙特卡洛仿真次数直接决定结果是“趋近真值”还是“一个偶然的随机数”。我做了收敛性实验固定SNR-10dBMONte Carlo次数从50递增到5000看检测概率的波动范围。实测数据显示蒙特卡洛次数在100以下是灾难性的检测概率在0.65到0.85之间乱跳完全没法用到500次时波动收窄到±0.02以内到2000次以后基本稳定在0.88附近波动不超过±0.005。所以如果你只是项目演示500次勉强够用但要是发论文或者做严谨的工程评估建议至少2000次起步。不要心疼那几分钟的运行时间数据可靠性远比时间值钱。4.3 采样点数对PR统计量分布的影响采样点数是一个比较隐蔽的关键参数。我在仿真里专门跑了一组不同采样点数下概率密度分布的对比发现一个有意思的规律纯噪声条件下PR值的经验分布均值基本不变但方差随着采样点数增加明显收窄。从方差角度看Ns64时PR值在[0.65, 0.85]区间分散Ns256时收窄到[0.70, 0.78]Ns1024时几乎集中到[0.73, 0.76]。这意味着采样点数越多噪声条件下PR分布越集中门限可以压得越低同样的虚警概率下检测灵敏度就越高。这个现象跟能量检测的结论类似感知时间越长统计量方差越小性能越好。实操上的一个启发是如果感知时隙受限、采样点数上不去可以通过增加协作节点数来补性能。3节点Ns128的方案和单节点Ns1024的方案性能差不多这给系统设计提供了灵活性。5. 仿真中容易翻车的细节噪声归一化、参考分布失配与循环变量冲突5.1 复噪声功率归一化的坑这个问题我在前面提了一下但值得单独展开因为太容易踩了。Matlab中randn(1, Ns)生成的实数序列方差是1如果你用下面这行代码生成复噪声noise randn(1, Ns) 1j * randn(1, Ns);那么噪声总功率是实部方差加虚部方差等于2。而信号功率是abs(signal).^2的均值计算信噪比时要拿信号功率除以噪声功率这个“2”就悄悄混进去了。如果信噪比定义是(SNR P_s / P_n)这个偏差会让所有仿真结果往偏大信噪比方向移动3dB。稳妥的写法是noise (randn(1, Ns) 1j * randn(1, Ns)) / sqrt(2);这样噪声总功率严格等于1信噪比定义干净利落。对于实信号基带模型对应的写法是noise randn(1, Ns)噪声功率就是1不用额外处理。5.2 参考分布失配的边界场景PR指数检测器在选择参考分布时有隐含假设。我用的是均匀分布作为参考这在噪声能量归一化后近似成立。但如果信道环境较复杂比如存在强窄带干扰、脉冲噪声纯噪声条件下的能量分布可能明显偏离均匀分布PR指数的基线就会抬升导致虚警概率系统性增大。我的解决思路是加一个“背景学习”环节在系统初始化阶段预留一段只含噪声的时间窗口实测该环境下的PR指数经验分布用这个分布替代理论均匀分布来定门限。这样PR指数检测器就变成半自适应的既能对抗噪声不确定性又能适配具体环境。在Matlab里实现这个思路很简单把第一阶段生成纯噪声样本的代码改成读取一段预采集的噪声数据就行。这也更贴近实际部署场景——认知无线电开机时先侦听一段时间感知背景噪声再进入正常工作状态。5.3 循环变量与Matlab内置函数命名冲突写Matlab仿真有一个新手频率极高的失误用mean、sum、std这些内置函数名当自定义变量。我在测融合性能时曾顺手写了个变量叫mean_PR结果后续代码里所有调用mean()的地方全部失效报错信息还很难定位。排查了半天才发现是变量名把内置函数遮蔽了。建议从一开始就养成习惯变量名用avg_PR、sum_PR这种带下划线的形式避免覆盖内置函数。另外ecdf返回的f和xi变量名在循环里反复使用要注意每次调用前清空或重新赋值否则上一次的残值可能会在极端情况下污染本次计算。5.4 蒙特卡洛仿真中的随机数种子管理可复现性是仿真项目的基本要求。如果不设置随机数种子每次运行的结果都不一样调试时很痛苦。我的做法是在主脚本开头固定种子rng(2024);这样每次运行的结果严格一致方便对比参数变更前后的差异。真正做性能评估时再放开种子跑多次取平均避免单次随机种子带来的偏差。6. 性能上界分析与多径衰落场景下的进一步扩展6.1 节点数增加带来的协作增益到底有多大很多人直觉上认为节点数越多越好但实测数据显示协作增益是边际递减的。我在SNR-10dB、Ns256条件下跑了一组节点数扫描协作节点数等增益融合检测概率相对2节点的增益10.613-20.80118.8%30.8838.2%40.9244.1%50.9482.4%从2节点到3节点提高了8个百分点但从4节点到5节点只提高了2.4个百分点。这个趋势说明在节点数超过4以后再堆节点数对性能的提升已经很小系统的瓶颈转移到了融合中心的合并算法和节点间的阴影衰落相关性上。实际系统设计时建议结合通信开销一起权衡。每个节点上报一次PR值需要传输一个浮点数如果网络是窄带物联网节点多了融合中心的处理压力和网络时延都会上来。我个人的经验是3到5个节点是性价比最高的区间。6.2 多径频选信道下的分集融合思路标准仿真模型里每个节点只有一个独立的瑞利衰落系数这属于平坦衰落。实际宽带通信系统遇到的多径频选信道更复杂不同频点的衰落各不相同。这时候每个节点的PR值只反映了该节点在感知频段上的综合能量聚合情况融合中心拿到多个节点的PR值后做平均等效于在空间维度做了分集。如果想进一步提升可以在时间维度也做分集每个节点在一个感知周期内做多次PR测量上报PR值的均值或中位数这能压制快衰落带来的瞬时波动。这个扩展实现很轻量就是在外层加一个重复测量循环但效果显著在慢变信道下检测概率能再涨3到5个百分点。6.3 从仿真到工程落地的差距提示Matlab仿真能证明算法有效性但真要做原型验证还有几个实际工程问题要考虑。首先PR指数计算依赖经验CDF需要积累足够多的采样点才能稳定。如果感知时隙很短采样点数不足PR检测的优势会缩水。其次融合中心和各节点之间的同步、数据上报时延在仿真里是零开销但实际网络中会造成判决延迟需要对感知周期做余量设计。我在用软件无线电平台做试验时还发现一个问题接收机前端的自动增益控制AGC会动态调整信号幅度直接影响能量归一化的结果。如果AGC收敛时间比感知周期长PR值会出现明显的周期波动。处理办法是在PR指数计算里用长时间窗做慢速归一化而不是用单帧的最大最小值归一化。7. 参数整定经验与代码获取说明7.1 一套可复用的参数配置模板前前后后跑了几百组仿真之后我沉淀了一套比较稳妥的默认参数模板直接抄作业基本不会出大问题参数推荐值说明采样点数 Ns256低于128时PR检测优势大幅缩水蒙特卡洛次数 M2000低于500时结果波动不可接受目标虚警概率0.05可根据监管要求调整协作节点数3~5超过5个性能增益边际递减融合方式等增益合并优先加权合并在信噪比估计可靠时再用随机数种子固定值确保调试时可复现7.2 代码结构和使用流程完整代码分三个文件结构是这样的PR_CSS_main.m主脚本定义参数、调用仿真函数、绘图compute_PR_metric.m单节点PR指数计算函数simulate_cooperative_PRI.m协作频谱感知完整仿真函数内部包含门限计算和性能统计。使用流程很直接先设置参数块里的Ns、M、Pfa_target、Nnodes然后直接运行主脚本脚本会依次完成门限标定、信噪比扫描、性能对比三个步骤最终输出检测概率曲线和ROC曲线。7.3 调整虚警概率时的联动操作更改目标虚警概率时要特别注意门限是自动联动计算的不需要手动修改阈值常量。这意味着你可以很方便地画出同一信噪比下不同虚警概率对应的检测概率曲线生成ROC曲线。不过有一点要提醒虚警目标压得太低时比如0.001门限会推到PR经验分布极右侧的尾部。这时候蒙特卡洛次数不够的话门限估计的方差会很大导致实测虚警概率偏离目标值。解决办法是虚警概率设得越低蒙特卡洛次数越要多建议满足(M \times Pfa_{target} \geq 100)的经验法则。8. 初始版本的一个小改动能带来额外增益最后分享一个我在测试过程中发现的、简单但有效的优化把融合中心等增益合并直接换成中位数合并。实现就一行代码PR_fusion median(PR_node);中位数合并对抗异常值的鲁棒性比均值好得多。在其中一个节点受到突发脉冲干扰时均值合并的融合统计量会被这个异常PR值明显抬升导致虚警概率恶化中位数合并则天然把异常值丢弃了。实测在单节点受脉冲干扰的模拟场景中中位数合并比均值合并的虚警概率低了一半以上。代价是性能天花板稍微比均值低一点点但通常在一个百分点以内。考虑到实际环境中个别节点出故障或者被干扰的概率不低我觉得中位数合并是一个很划算的保底策略。项目里可以同时保留两种合并方式通过一个枚举参数切换方便后续在更多信道场景里做系统性对比。PR检测器这条路其实还远没到头采样点数自适应、融合权重的在线学习、与机器学习分类器的结合都有不少可探索的空间。先把基础链路跑通后面每一步迭代都会比重新造轮子快得多。