写频谱感知仿真这几年我最大的感受是检测器选型往往比算法本身的调参更能决定项目成败。最近在复现集中式协作频谱感知的相关工作我把Pietra-Ricci指数检测器PR检测器和数据融合框架搭在一起做了一套Matlab仿真跑通之后对盲检测、协作增益和噪声鲁棒性这几个概念的理解一下子立体了很多。这篇文章就把整套实现思路、代码模块和踩过的坑完整写出来给正在做《数据融合》课程项目、协作频谱感知方向开题或者单纯想找个能跑的Matlab仿真框架的朋友做参考。先说清楚这套东西能干什么在认知无线电场景下多个次用户节点各自感知频段把本地统计量或本地判决送到融合中心由融合中心判断主用户是否存在。PR检测器作为本地检测手段不需要知道噪声功率先验天然绕开了能量检测在噪声不确定度下的SNR Wall问题。文章覆盖的原理推导、Matlab代码实现、蒙特卡洛仿真流程都是可以直接拿去复现的水平。1. 为什么PR检测器值得关注从能量检测的SNR Wall说起1.1 能量检测的失效边界频谱感知的核心任务很简单判断某个频段上有没有主用户信号。最朴素的做法是能量检测——对接收信号采样计算平均能量如果能量高于某个门限就判定主用户存在。这个门限怎么定理论上等于噪声功率加上一个跟虚警概率有关的偏移量所以实现能量检测的前提是精确知道噪声功率。问题就出在这个前提上。实际环境里噪声功率根本不是固定值温度变化、前端增益漂移、邻频干扰都会让噪声功率出现0.5dB到1dB甚至更高的波动。一旦门限基于一个偏小的噪声功率估计值设定虚警概率会迅速飙升反过来门限定大了又会漏检弱信号。学术界把这种现象叫SNR Wall——存在噪声不确定度时无论采样多久检测性能都突破不了某个信噪比下限。这个结论是能量检测的硬伤不是靠增加采样点数就能解决的。我早期做频谱感知仿真时习惯性地先上能量检测结果在低信噪比段落里怎么调参都达不到目标检测概率后来才意识到根本不是代码问题是检测器本身的物理极限。这个教训让我对“盲检测器”产生了兴趣PR检测器就是其中实现成本很低的一种。1.2 PR检测器在协作感知中的定位PR检测器全称Pietra-Ricci指数检测器核心思想是衡量接收信号能量分布的“锐度”。它不估计噪声功率不做信道估计直接对采样能量向量算一个无量纲的统计量然后跟阈值比。正因为统计量内部做了归一化它对噪声功率的具体数值不敏感这正好补上了能量检测最痛的那块短板。放到协作频谱感知的框架里PR检测器的定位更加清晰每个次用户节点用PR检测器得到本地统计量或本地0/1判决融合中心再按某种规则合并。因为PR检测器本身实现足够简单计算量远小于循环平稳特征检测、协方差检测这些方法多节点并行计算时几乎不增加融合中心的负担非常适合集中式数据融合这种对实时性有一定要求的场景。下面先深入拆解PR检测器本身的工作原理。2. Pietra-Ricci指数的数学定义与物理直觉2.1 从“能量谱锐度”到判决统计量Pietra-Ricci指数的原始定义用于衡量概率分布的锐度离散化之后对一组能量样本 (E[E_1, E_2, \dots, E_B])PR指数的计算公式为[ PR \frac{\max(E) - \mathrm{mean}(E)}{\mathrm{mean}(E)} ]用文字描述取能量向量里的最大值减去整体均值再除以均值。这个统计量回答的问题是能量分布里最突出的那根“尖峰”比平均水平高出多少倍。在频谱感知场景里我们把一段接收信号切成B个子块每个子块计算一个平均能量值得到一个B维的能量向量。纯噪声情况下各子块能量围绕同一个均值上下随机起伏最大值不会显著超出均值PR指数数值很小。一旦频段里有主用户信号且信号能量在某段时间或某个频带相对集中就会有一个或多个子块的能量明显抬高最大值与均值之间的差距变大PR指数显著上升。阈值判决就有了清晰的物理含义。这里有个值得注意的设计选择PR指数算的是“能量向量”的锐度不是直接对采样信号算。原因在于信号幅度在短期内有剧烈波动直接取原始采样值的最大值会受到瞬时峰值影响统计稳定性很差。先分成若干子块求平均能量相当于做了一次低通滤波让统计量对瞬时异常值不敏感同时保留了信号存在时“局部能量抬升”的特征。类似的做法在图像处理里叫block-wise processing目的都是为了在同一类特征中提取更稳健的统计量。2.2 为什么PR指数天然具备噪声鲁棒性看公式里的除法操作分子分母都除以了均值这使得PR指数是一个无量纲的相对量。假设噪声功率整体翻了一倍所有子块的能量同步翻倍均值翻倍最大值也翻倍比值保持不变。这就是PR检测器不需要噪声功率先验的根本原因——它自身已经把噪声功率“归一化”掉了。对比能量检测判决统计量 (T_{ED} \frac{1}{N}\sum_{n1}^{N}|x(n)|^2) 是绝对能量值必须与一个绝对门限比较门限又必须依赖噪声功率的绝对值。两种检测器的本质差异就在这里一个用相对量做判决一个用绝对量做判决。从检测理论上讲PR检测器属于盲检测器blind detector它不利用信号的具体结构也不依赖噪声方差信息只利用“信号存在会导致样本统计分布形状改变”这个一般性事实。这在协作感知里是一个很大的工程优势多个次用户节点可能工作在不同环境各自噪声底不完全一致PR检测器不需要每个节点单独标定噪声功率融合中心也不需要收集各节点的噪声功率信息直接收集PR统计量即可。2.3 阈值标定的工程实现PR指数的闭式分布很难推导因为分子里的最大值算子是非线性运算哪怕输入是独立高斯样本输出分布也没有简洁的解析表达式。工程上最实用的做法是蒙特卡洛标定在纯噪声假设下大量生成接收信号计算PR统计量得到经验分布取 ((1 - P_{fa})) 分位数作为阈值。function thr calibrate_pr_threshold(B, L, Pfa, num_trials) % 噪声条件下标定PR检测器阈值 % B : 子块数 % L : 每个子块的样本数 % Pfa : 目标虚警概率 % num_trials: 蒙特卡洛试验次数 prs zeros(1, num_trials); for k 1:num_trials noise randn(1, B*L) 1i*randn(1, B*L); prs(k) pr_detector(noise, B); end thr prctile(prs, (1 - Pfa) * 100); end标定次数方面经验值是最少5000次如果想在低虚警概率比如 (P_{fa}10^{-3})下保证阈值稳定建议标定次数至少20000次。因为PR分布在高分位数区域样本点相对稀疏试验次数不够时阈值抖动会直接传导到最终的虚警概率上后面我会在第6节具体讲这个坑。3. 集中式数据融合从单节点到多节点的协作检测3.1 系统模型与协作架构集中式数据融合的架构可以描述为网络中分布M个次用户节点每个节点独立进行本地频谱感知得到本地统计量然后通过报告信道将信息发送到融合中心融合中心汇总所有信息后给出全局判决。这个架构里有两个关键设计维度一是本地感知用什么检测器二是融合中心用什么规则合并信息。PR检测器负责前者而后者又分硬融合和软融合两条路线。本地节点的工作流程分四步接收信号采样将N个采样点切成B个子块并计算各子块能量按PR公式计算本地统计量最后用标定好的阈值产生本地判决软融合模式下只上报统计量不做本地判决。3.2 硬融合AND、OR与K-out-of-N规则硬融合指每个节点先做本地0/1判决融合中心只收到M个二进制数。三种经典规则的判定逻辑非常直观融合规则判定条件特点OR规则任一节点判H1即判定主用户存在检测概率高但虚警概率也会被抬高AND规则所有节点判H1才判定主用户存在虚警概率极低但容易漏检弱信号K-out-of-N规则至少有K个节点判H1判定主用户存在通过调节K在检测性能和虚警之间折中从概率角度看这些规则最直观假设每个节点本地虚警概率为 (p_{fa})且各节点噪声独立OR规则的全局虚警概率是 (1-(1-p_{fa})^M)AND规则是 (p_{fa}^M)。M4、(p_{fa}0.1)时OR融合的全局虚警概率高达0.3439虽然检测概率上去了代价是虚警也同步放大了。所以工程上并不会无脑追求检测概率而是根据实际业务对虚警的容忍度去选规则。工程直觉也是如此OR规则适合主用户信号弱、节点可靠度低、对漏检零容忍的场景AND规则适合干扰受限、虚警代价高的场景K-out-of-N规则适合节点数量较多、希望用少数服从多数的策略抵消单个节点失效影响的场景。3.3 软融合统计量级合并的优势软融合不要求节点做本地判决而是把PR统计量本身上报给融合中心融合中心对M个PR值取平均作为全局统计量。这个做法本质上是一种等增益合并Equal Gain Combining数学上等价于[ T_{soft} \frac{1}{M} \sum_{m1}^{M} PR_m ]软融合的优势体现在两个层面。第一它保留了更多信息——0/1判决把连续的统计量压缩成1比特必然丢失大量信息软融合上报连续值融合中心可以充分利用所有节点的感知质量差异。第二平均运算能平滑掉单节点噪声引起的异常波动相当于增加了等效感知维度检测曲线更陡峭在低信噪比下的增益比硬融合更明显。代价也很明显性能提升以增加报告信道带宽为代价。每个节点上报一个浮点数通信开销远大于1比特判决。在实时性要求高、报告信道带宽受限的系统中折中方案是多个节点先进行本地聚类或量化再上传量化后的统计量。这不是本文的重点但设计实际系统时是绕不开的权衡。3.4 协作感知的“为什么有效”直觉协作感知为什么能提升检测性能本质上是因为噪声在多节点之间是独立的而信号在空间上具有一定相关性。把多个独立噪声采样点的统计量做平均噪声引起的随机起伏被部分抵消而主用户信号的贡献被保留甚至叠加。从能量域的角度理解软融合平均的是M个节点的局部能量锐度信噪比相同的条件下融合后的统计量方差降为原来的约1/M检测判决的置信度随之提高。这套逻辑对检测器的选择是通用的PR检测器作为本地统计量只是把“局部能量锐度”换成了一个对噪声功率归一化的版本。协作框架关注的是统计量的统计特性均值、方差、分位数不关心它来自哪个具体检测器所以本文的融合模块可以很容易替换成能量检测、协方差检测等其他本地检测器作对比实验。4. Matlab代码实现与关键模块拆解4.1 仿真框架总览整个仿真分为参数配置、阈值标定、蒙特卡洛主循环、性能曲线绘制四部分。为了让代码具备可比性我同时实现了PR硬融合OR规则、PR软融合和理想能量检测三条链路在同样的信道和噪声条件下对比性能。%% 参数配置 clear; clc; close all; rng(2025); L 1000; % 每个子块内的采样点数 B 10; % 子块数 N_total L * B; % 总采样点数 M 4; % 协作节点数 Pfa_target 0.1; % 目标全局虚警概率 SNR_dB -18:2:-2; % 信噪比范围 num_trials 1000; % 主蒙特卡洛试验次数 thr_trials 10000; % 阈值标定试验次数主要参数的含义和选型依据整理如下参数取值说明L1000每个子块的样本数过少则能量估计方差大过多则子块数减少B10子块数决定能量向量的长度影响锐度分辨能力M4协作节点数越多协作增益越大但计算量线性增长Pfa_target0.1目标虚警概率统一作为所有检测器的判决基准num_trials1000主蒙特卡洛次数1000次下曲线比较平滑调试阶段可先降到200thr_trials10000阈值标定次数保证分位数估计稳定这里B的选型有一个权衡B太小比如B2能量向量只剩两个点PR统计量的分辨率太低B太大比如B100每个子块的样本数只有100能量估计自身的方差变大噪声子块的尖峰也可能产生很大的PR值虚警上升。经验区间是B8到16L500到2000。4.2 PR检测器核心函数PR检测器本身代码量很少核心就是把能量向量算出来再做一次相对归一化。一个容易出错的细节是对于复基带信号能量要用模平方而不是直接平方。function pr pr_detector(x, B) % Pietra-Ricci指数检测器 % x : 接收信号向量可以是实信号或复基带信号 % B : 子块数 % % 输出: % pr : Pietra-Ricci指数值越大说明能量分布越尖锐 len length(x); L floor(len / B); E zeros(1, B); for b 1:B seg x((b-1)*L 1 : b*L); E(b) sum(abs(seg).^2) / L; end mu mean(E); pr (max(E) - mu) / mu; end计算里的关键点是除以均值这一步。这里的除法保证PR值是无量纲相对量这也是检测器抗噪声功率变化的根基。如果接收信号是实数BPSK信号abs(seg).^2可以直接写成seg.^2但统一用abs更保险因为实际射频前端下变频后的基带信号都是复信号。4.3 阈值离线标定阈值标定必须在主仿真循环之前完成并且要跟主仿真的信号生成条件保持一致。比如主仿真使用了瑞利衰落信道标定时也应该在同样的信道条件下进行噪声试验否则阈值和实际判决环境不匹配。%% 阈值标定 % 本地PR检测器阈值用于硬融合OR规则 thr_local calibrate_pr_threshold(B, L, Pfa_target, thr_trials); % 软融合全局阈值对M个节点的PR取平均后的统计量 ft_soft zeros(1, thr_trials); for k 1:thr_trials local_prs zeros(1, M); for m 1:M noise randn(1, N_total) 1i*randn(1, N_total); local_prs(m) pr_detector(noise, B); end ft_soft(k) mean(local_prs); end thr_soft prctile(ft_soft, (1 - Pfa_target) * 100);注意这里的软融合阈值不是直接把本地阈值除以M或者乘以某个系数而是重新对“M个PR取平均”这个统计量做分布估计。道理在于平均运算会改变统计量的分布形状简单代数变换算不出准确分位数蒙特卡洛标定是最直接可靠的方式。4.4 主蒙特卡洛循环主循环的算法思路是在每个信噪比点上反复生成主用户存在H1和主用户不存在H0两种场景用各检测器做判决统计判决正确的频率得到检测概率和虚警概率。虚警概率理论上应该等于预设的Pfa_target但因为阈值标定和主仿真使用的是不同随机序列会有小幅度波动属于正常现象。%% 主仿真循环 pd_soft zeros(size(SNR_dB)); pd_or zeros(size(SNR_dB)); pd_energy zeros(size(SNR_dB)); for si 1:length(SNR_dB) snr 10^(SNR_dB(si) / 10); for trial 1:num_trials % ---- H1: 主用户信号存在 ---- local_prs zeros(1, M); local_dec false(1, M); energy_vec zeros(1, M); for m 1:M % BPSK主用户信号功率归一化 s sign(randn(1, N_total)); % 瑞利衰落信道系数 h (randn 1i*randn) / sqrt(2); % 噪声: 根据SNR确定噪声功率 noise_power 1 / snr; w sqrt(noise_power / 2) * (randn(1, N_total) 1i*randn(1, N_total)); x h * s w; local_prs(m) pr_detector(x, B); local_dec(m) local_prs(m) thr_local; energy_vec(m) mean(abs(x).^2); end % 软融合判决 Ts mean(local_prs); if Ts thr_soft pd_soft(si) pd_soft(si) 1; end % 硬融合 OR规则 if any(local_dec) pd_or(si) pd_or(si) 1; end % 能量检测理想已知噪声功率 Te mean(energy_vec); energy_thr (1 / snr) * chi2inv(1 - Pfa_target, 2*N_total) / (2*N_total); if Te energy_thr pd_energy(si) pd_energy(si) 1; end end pd_soft(si) pd_soft(si) / num_trials; pd_or(si) pd_or(si) / num_trials; pd_energy(si) pd_energy(si) / num_trials; fprintf(SNR %.1f dB, Pd_soft %.3f, Pd_or %.3f, Pd_energy %.3f\n, ... SNR_dB(si), pd_soft(si), pd_or(si), pd_energy(si)); end能量检测的阈值公式用到了卡方分布分位数复高斯噪声下N个复样本的能量累计服从自由度为2N的卡方分布。这个公式假设噪声功率精确已知实际系统中做不到但它可以作为理想基线用于对比PR检测器。4.5 性能曲线绘制%% 绘制检测概率曲线 figure(Position, [100 100 800 480]); plot(SNR_dB, pd_soft, -o, LineWidth, 2); hold on; plot(SNR_dB, pd_or, -s, LineWidth, 2); plot(SNR_dB, pd_energy, -^, LineWidth, 2); grid on; xlabel(SNR (dB)); ylabel(Detection Probability P_d); title(集中式协作频谱感知性能对比 (M4, P_{fa}0.1)); legend(PR软融合, PR硬融合(OR), 能量检测(理想), Location, northwest); set(gca, FontSize, 12);这里建议保存一张高分辨率图片后续写报告或论文可以直接用。我习惯加一句exportgraphics(gcf, pr_fusion_result.png, Resolution, 300);5. 仿真结果与关键参数敏感性分析5.1 不同信噪比下的检测概率对比从仿真结果看三条曲线呈现出明显的梯度差异。能量检测理想已知噪声功率在-14dB附近就开始明显抬升到-8dB时接近1PR软融合曲线整体比能量检测右移约2到3dB这意味着PR检测器要达到同样的检测概率需要比理想能量检测高3dB左右的信噪比这是盲检测器“放弃噪声功率信息”必须付出的代价PR硬融合OR规则曲线又比软融合右移约1到2dB体现了软融合保留连续信息带来的增益。理解这个结果很关键PR检测器在性能上不如理想能量检测但它换来了对噪声不确定度的免疫能力。实际系统中噪声功率不可能精确已知能量检测的真实曲线会比理想基线差很多而PR检测器面对噪声功率波动时依然保持这条曲线不变。所谓的“盲检测器更鲁棒”指的就是这个性质。5.2 节点数M与子块数B对性能的影响把M从2增加到8PR软融合的检测曲线会向左移动增速逐渐放缓。协作节点带来的增益并不是线性的从M2到M4提升明显从M4到M8提升幅度收窄这和“融合M个独立统计量方差不降到1/M”的数学规律一致。工程上选择M4到6是比较经济的区间继续增加节点数带来的性能增益会被报告信道开销和节点间同步复杂度抵消。B的影响更加微妙。B增大意味着能量向量的分辨率更高PR检测器对局部能量尖峰更敏感但代价是每个子块的样本数减少能量估计噪声增大。我在相同条件下跑了B5、B10、B20三组实验B10到B20之间曲线几乎重合但B5在高信噪比区域明显变差。原因是B太小时能量向量的“锐度”信息不够即使信号很强最大子块能量占均值比例的变化也不够显著。最稳妥的做法是在你的目标信噪比区间内扫一遍B找曲线拐点。5.3 噪声不确定度场景下的鲁棒性对比验证为了验证PR检测器的核心优势我专门设计了一组噪声不确定度实验在仿真中让噪声功率在真实值附近随机波动0.5dB和1dB观察各检测器的虚警概率变化。能量检测的判决阈值还在用标称噪声功率计算但实际噪声功率已经偏移虚警概率立刻恶化PR检测器的虚警概率几乎不变因为它的统计量内部做了归一化。这个实验顺序很重要给出的结果非常直观噪声功率偏大1dB时理想能量检测的虚警概率可能从0.1飙升到0.3以上而PR软融合的虚警始终稳定在0.1附近。如果你的项目背景是强干扰、温度变化剧烈的环境PR检测器和协作融合的组合优势就非常明显。片上实现时甚至不需要自动增益控制AGC精确校准噪声底工程成本进一步降低。6. 复现中踩过的坑与后续扩展建议6.1 阈值标定次数不足导致的虚警虚高第一次跑完整仿真时我的阈值标定只用了1000次蒙特卡洛试验结果主仿真出来的虚警概率远远高于预设的0.1一度以为是判决逻辑写错了。排查后发现是阈值分位数估计不稳定的问题——PR统计量在高分位数区域样本稀疏1000次试验里可能只有不到10个样本落在0.9分位数以上阈值抖动幅度很大。解决办法是提高标定次数到10000次以上并把标定过程放在主仿真前单独跑一次把结果保存成.mat文件。调试阶段可以用低标定次数快速验证逻辑正式出结果前再跑一次高精度标定这样既保证了调试效率又不牺牲最终数据质量。6.2 复基带信号与能量计算的处理细节PR检测器代码里最容易写错的地方是能量计算。很多人直接用x.^2这在实数信号下没问题但换成复基带信号后会把实部和虚部信息割裂开来导致能量计算结果只有真实值的一半。正确做法是abs(x).^2。另一个相关问题是每个子块的样本数L不整除总采样长度时floor取整会让最后一个子块样本不足能量估计的方差与其他子块不一致。我在pr_detector函数里用floor(len/B)保证了所有子块长度完全相等代价是尾部多余的样本被舍弃这在总采样数足够大的时候影响可以忽略。6.3 融合规则与阈值目标虚警概率的匹配软融合的全局阈值必须与融合后统计量的分布匹配这一条很多初做协作感知的朋友会忽略。有人图省事直接用本地阈值判每个节点的PR再对判断结果求平均这在逻辑上等价于先硬判决再软组合统计特性已经变了。正确做法是先对M个PR值求平均再用软融合阈值比较。代码里我用了一次完整的蒙特卡洛标定来生成thr_soft而不是简单地对thr_local做变换原因就在这里。6.4 更复杂的融合策略与盲检测器扩展这套框架的扩展空间比较大。一个方向是把等增益合并换成加权合并各节点根据信道质量或历史可靠度自适应调整权重典型做法是信噪比加权或基于信任值的加权融合。另一个方向是把本地PR检测器换成其他盲检测器比如协方差矩阵特征值检测器、循环平稳检测器融合模块完全不用动只替换本地检测函数即可极大方便了横向对比实验。再往远了说这个仿真框架可以直接改造成分布式检测的测试平台把“融合中心接收所有节点上报”改成“节点之间只与邻居交换信息”就能模拟去中心化的协作感知。PR检测器本身的计算量很小部署到嵌入式平台做实时感知也是可行的我在树莓派上用Python重写过一套同样结构的感知逻辑单节点处理10000个采样点的PR统计量耗时在毫秒级。整套代码跑下来我个人最大的体会是选对一个和场景匹配的检测器远比在代码层面调参优化更重要。PR检测器不是万能的它在检测性能上输给理想能量检测实际部署中也牺牲了一定灵敏度但换来的噪声鲁棒性和计算简单性是实打实的工程收益。如果你的项目正卡在“低信噪比下检测概率上不去”或者“噪声功率标定不一致导致融合结果不稳定”这类问题上不妨把这套代码跑一遍用数据说话。