前阵子在做集中式协作频谱感知的仿真对比把几种基于分布距离的检测器都过了一遍。Pietra-Ricci指数检测器是其中让我印象比较深的一个思路非常朴素代码量不大性能却不输传统的能量检测尤其在低信噪比区间和噪声功率估计有波动时反而有优势。这篇文章就把从原理到Matlab实现、再到调试中踩过的坑完整梳理一遍。如果你在做认知无线电方向或者正被协作频谱感知的融合策略、检测统计量设计搞到头大这篇应该能给你一个可以直接抄作业的参考。1. 为什么要做协作感知又为什么选集中式融合与P-R指数1.1 单节点感知的三大硬伤频谱感知的本质是判断某个频段上是否有主用户信号。听起来简单但单节点做这件事在实际场景里有几个绕不开的问题。最典型的是隐藏终端问题。主用户信号可能被建筑物遮挡某个次用户恰好处于阴影区接收信噪比很低本地检测根本感知不到信号的存在。它把“没信号”的结论上报之后系统可能就认为这个频段是空闲的于是安排次用户接入结果次用户反而对正在通信的主用户产生了干扰。这个问题靠增加单个节点的灵敏度是解决不了的因为盲区是空间位置决定的。第二个问题是多径衰落。无线信道里接收信号是多个反射路径叠加的结果在某些频点上可能深衰落十几甚至几十dB。哪怕节点离主用户不远只要当前信道处于衰落谷底单次感知的判决质量就会非常差。第三个问题是噪声不确定性。接收机前端的热噪声功率并不是恒定不变的会随温度、器件状态、邻频干扰而波动。传统能量检测的门限完全建立在噪声功率的精确估计上噪声功率只要偏一点虚警概率和检测概率就会剧烈恶化。这三个问题叠加在一起单节点感知的可靠性很难满足实际系统的要求。协作感知就是把多个空间上分散的节点联合起来利用空间分集来对抗阴影和衰落同时通过多个节点的观测来抑制噪声不确定性的影响。1.2 集中式融合信息最完整实现最省心协作感知的融合架构大体分两类集中式和分布式。集中式架构下各次用户把感知数据通过控制信道送往融合中心由融合中心统一判决分布式架构下节点之间互相交换信息经过多轮迭代收敛到一致判决。从工程实现角度看集中式有明显优势融合中心能拿到全部节点的观测信息理论上可以获得最优检测性能全局判决只需要一次不需要迭代协商时延可控节点之间不需要建立复杂的网状通信关系控制信道只到融合中心拓扑简单。当然集中式也有代价所有数据都汇聚到融合中心控制信道带宽压力大融合中心是单点。这也是很多研究转向分布式的原因。但在“先把算法和检测器验证清楚”这个阶段集中式是最合适的起点——数据融合的逻辑直观检测统计量的设计空间大跟分布式相比少了一层迭代收敛的不确定性调试起来不容易被干扰因素带偏。1.3 为什么是P-R指数而不是传统能量检测能量检测是频谱感知的“默认选项”实现简单、不需要信号先验但它的判决门限是贴在噪声功率上的。噪声功率估计偏一点门限就失效了。这是能量检测在实用场景里最头疼的问题。P-R指数检测器走的是另一条路。它不看单个能量值是否超过某个绝对门限而是看“当前观测能量值的整体分布”和“纯噪声下的理论分布”之间的差异有多大。只要观测分布明显偏离噪声分布就判定存在信号。这种基于分布差异的检测思想对单个能量值的波动不敏感对噪声功率小幅偏移也有一定容忍性因为它是拿整条累积分布曲线在做比较。之前常被拿来比较的分布距离还有KL散度、巴氏距离、柯尔莫哥洛夫-斯米尔诺夫统计量等。KL散度不对称巴氏距离涉及根号运算实现稍繁KS统计量只取CDF最大差异对中段分布的整体偏离不敏感。P-R指数用的是两条CDF曲线之间夹的“面积”是一个全局性的差异度量计算上又比巴氏距离简单得多实测下来稳定性很好。2. Pietra-Ricci指数检测器的原理拆解2.1 P-R指数的数学定义和几何直觉Pietra指数最早出自经济学用来衡量收入分配的不平等程度把人口的累计占比和收入的累计占比画在一起得到洛伦兹曲线曲线和完全平等线之间夹的面积越大说明分配越不平等。把这个思想搬到检测问题里就是比较两个概率分布的重合程度。假设 (H_0) 下的概率密度是 (f_0(x))(H_1) 下是 (f_1(x))P-R指数定义为[ J(f_0,f_1)\frac{1}{2}\int_{-\infty}^{\infty}|f_0(x)-f_1(x)|dx ]这里前面乘的 (1/2) 是归一化因子因为两个密度函数差的绝对值积分最大是2除2之后 (J) 的取值范围就是 ([0,1])。当 (J0) 时两个分布完全重合说明观测数据和噪声分布没有差异大概率是 (H_0)当 (J1) 时两个分布完全不重叠说明观测数据已经彻底偏离噪声分布大概率是 (H_1)。从几何上看这个值就是两个概率密度曲线之间“不重叠部分的面积取一半”。两个分布离得越远重叠越少(J) 越大检测越容易。2.2 从分布距离到检测统计量检测器真正能用P-R指数的前提是我们在 (H_0) 下知道噪声分布的理论表达式。这是标准假设接收机的噪声建模为循环对称复高斯功率 (\sigma^2) 可以通过静默期估计得到。那么纯噪声下能量统计量的理论分布是可以精确写出来的。实际感知时我们把 M 个次用户上报的本地能量值收集起来构造一条经验CDF然后和理论噪声CDF计算P-R指数。如果主用户存在能量值的分布会整体右移经验CDF和理论CDF之间会拉开明显距离P-R指数变大如果只有噪声经验CDF会贴着理论CDF走指数很小。这里有个工程细节我很早就踩过坑理论定义用的是概率密度函数 (f)但在数值实现里密度函数不好估计。直方图要选bin宽度核密度估计要选核函数和带宽选不好就会引入人为偏差而且M个用户带来的M个能量值样本量通常不大直方图会很稀疏。换成CDF就不一样了经验CDF直接由排序后的样本点构造没有参数选择问题理论CDF又是现成的闭式表达式数值积分稳定得多。所以我在实现里用“CDF夹层面积”作为P-R指数的实际计算方式效果等同而且鲁棒。2.3 理论噪声能量分布卡方分布与自由度匹配每个次用户采集N个采样点本地能量定义为[ E_i\frac{1}{N}\sum_{n1}^{N}|y_i(n)|^2 ](H_0) 下 (y_i(n)w_i(n))(w_i(n)\sim\mathcal{CN}(0,\sigma^2))即实部和虚部都是独立的高斯变量方差各 (\sigma^2/2)。把能量表达式重新整理[ \frac{2NE_i}{\sigma^2}\sum_{n1}^{N}\left|\frac{\sqrt{2}w_i(n)}{\sigma}\right|^2 ]括号里每一项是一个归一化复高斯变量模平方服从自由度为2的卡方分布。N项求和之后(\frac{2NE_i}{\sigma^2}) 服从自由度为 (2N) 的中心卡方分布。所以在融合中心(H_0) 下单个次用户能量 (E_i) 的理论CDF是[ F_0(x)F_{\chi^2_{2N}}\left(\frac{2Nx}{\sigma^2}\right) ]Matlab里直接调用chi2cdf(2*N*x/sigma2, 2*N)就能算出来。自由度这里特别容易搞错。如果接收机是实信号模型每个采样点只贡献1个自由度N个点对应自由度为N的卡方分布如果是复基带模型每个复采样点包含实部和虚部两个维度自由度要翻倍成 (2N)。我的仿真里都是复基带信号所以用的是 (2N)。你要是把自由度写错理论CDF和实际噪声能量分布对不上P-R指数在 (H_0) 下也会偏大虚警完全失控。2.4 为什么这个思路适合集中式数据融合P-R指数检测器的数据需求非常“轻”融合中心只需要每个次用户上报一个能量标量 (E_i)不需要原始采样点也不需要本地判决结果。M个能量值就构成了一组样本足够构造经验CDF。相比把所有原始数据搬到融合中心的方案控制信道带宽压力小很多相比硬判决融合每个节点只发1bit判决结果又保留了能量值的幅度信息信息损失小。还有一个对协作场景特别友好的点P-R指数作为融合统计量时节点之间完全不需要同步波形只要各自在同一个感知时隙内完成采样把能量发过来就行。这对实际系统的时频同步要求非常宽松。3. Matlab仿真框架与核心代码实现3.1 仿真参数与系统模型我要仿真的是这样的场景一个主用户信号被M个次用户感知每个次用户在感知时隙内采集N个复基带采样点然后各发一个能量值给融合中心融合中心算P-R指数并与门限比较。主用户信号用QPSK调制符号序列随机生成功率通过SNR设定。信道先做最简单的AWGN暂时不考虑衰落这样可以先把检测器的核心行为验证清楚。后续要扩展瑞利衰落也非常容易只要在信号前乘一个信道系数。仿真参数如下参数取值说明次用户数 M8融合中心收到的能量样本个数每用户采样点数 N128单个感知时隙内的采样数噪声功率 (\sigma^2)1归一化信噪比 SNR-12~0 dB步进2dB目标虚警概率 (P_f)0.05门限设计基准蒙特卡洛次数5000保证门限和检测概率稳定M取8比较有代表性既能体现协作增益又不至于让控制信道压力太大。N取128对应一个很短的感知时隙和实际系统的感知开销量级一致。3.2 核心代码P-R指数计算函数先把最核心的函数写出来。它做的事是输入M个能量值、噪声功率、每用户采样点数输出P-R指数估计值。function J computePRIndex(E, sigma2, N, gridPoints) % computePRIndex 计算Pietra-Ricci指数 % E : M x 1 能量向量融合中心汇总的各用户本地能量 % sigma2 : 噪声功率 % N : 每个用户在一个感知时隙内的采样点数 % gridPoints: CDF积分网格点数默认1000 % J : P-R指数估计值 if nargin 4 || isempty(gridPoints) gridPoints 1000; end E E(:); xmax max(max(E) * 1.2, sigma2 * 3); x linspace(0, xmax, gridPoints); % 经验CDF: 每个网格点上统计能量小于该值的比例 F_emp arrayfun((xx) mean(E xx), x); % 理论噪声能量CDF: 2*N*E/sigma2 ~ chi2(2N) F_th chi2cdf(2 * N * x / sigma2, 2 * N); % 两条CDF曲线之间的面积即P-R指数估计 J trapz(x, abs(F_emp - F_th)); end网格积分上限我设为“最大能量值再延伸20%”和“理论噪声能量均值的3倍”两者取大。这么做是为了保证积分范围覆盖CDF从0到1的完整变化区间。如果上限只取到最大观测能量尾部理论CDF接近1的部分就没算进去积分结果会偏小。F_emp arrayfun(...)在M只有8个、网格1000个点时完全没问题。如果M上百个同时网格又很密可以改成更高效的排序插值方式但实际场景M一般不会那么大这个实现已经够快。3.3 核心代码门限设计与蒙特卡洛流程P-R指数在 (H_0) 下的理论分布很难写闭式所以门限我用蒙特卡洛仿真确定。做法很直接在纯噪声条件下跑大量实验得到P-R指数的经验分布然后取它的 (1-P_f) 分位数作为门限。% 参数定义 M 8; N 128; sigma2 1; Pf_target 0.05; numMC 5000; snrVec -12:2:0; % 预分配结果 Pd_PR zeros(size(snrVec)); Pd_ED zeros(size(snrVec)); for s 1:length(snrVec) snr 10^(snrVec(s) / 10); Ps snr * sigma2; % 信号功率先跑H0生成噪声数据得到J0分布取分位数得到门限J0 zeros(numMC, 1); for mc 1:numMC noise sqrt(sigma2/2) * (randn(M, N) 1i * randn(M, N)); E mean(abs(noise).^2, 2); J0(mc) computePRIndex(E, sigma2, N); end lambda_PR quantile(J0, 1 - Pf_target);然后跑H1主用户信号用QPSKdetectPR zeros(numMC, 1); detectED zeros(numMC, 1); % QPSK基带符号功率归一化到Ps baseSym exp(1i * (randi([0 3], N, 1) * pi/2 pi/4)); signal sqrt(Ps) * baseSym; for mc 1:numMC E zeros(M, 1); for k 1:M % 每个用户收到同样的主用户信号叠加独立噪声 y signal sqrt(sigma2/2) * (randn(N, 1) 1i * randn(N, 1)); E(k) mean(abs(y).^2); end % P-R指数检测 J computePRIndex(E, sigma2, N); detectPR(mc) J lambda_PR; % 等增益合并能量检测作为性能对比 T_ED mean(E); lambda_ED chi2inv(1 - Pf_target, 2 * M * N) * sigma2 / (2 * M * N); detectED(mc) T_ED lambda_ED; end Pd_PR(s) mean(detectPR); Pd_ED(s) mean(detectED); end能量检测的对比门限用了解析式。等增益合并后的平均能量在 (H_0) 下满足 (\frac{2MN T_{ED}}{\sigma^2}\sim\chi^2_{2MN})所以chi2inv(1-Pf, 2*M*N)可以直接给出理论门限。这里要注意的是能量检测门限和P-R检测门限都按 (P_f0.05) 设计才能保证对比公平。3.4 画ROC曲线的方法固定SNR比如取-6dB要画一条ROC曲线做法是把门限从低到高遍历统计对应的虚警概率和检测概率。代码可以在主循环基础上改一下snrFix -6; Ps 10^(snrFix/10) * sigma2; % H0下的J分布 J0 zeros(numMC, 1); for mc 1:numMC noise sqrt(sigma2/2) * (randn(M, N) 1i * randn(M, N)); E mean(abs(noise).^2, 2); J0(mc) computePRIndex(E, sigma2, N); end % H1下的J分布 J1 zeros(numMC, 1); signal sqrt(Ps) * exp(1i * (randi([0 3], N, 1) * pi/2 pi/4)); for mc 1:numMC E zeros(M, 1); for k 1:M y signal sqrt(sigma2/2) * (randn(N,1) 1i*randn(N,1)); E(k) mean(abs(y).^2); end J1(mc) computePRIndex(E, sigma2, N); end PfGrid 0.01:0.01:1; PdCurve zeros(size(PfGrid)); for g 1:length(PfGrid) th quantile(J0, 1 - PfGrid(g)); PdCurve(g) mean(J1 th); end plot(PfGrid, PdCurve, LineWidth, 1.5); xlabel(虚警概率 P_f); ylabel(检测概率 P_d); grid on;ROC曲线上的每个点都对应一个门限值实际系统只需要按目标 (P_f) 选一个工作点即可。4. 实验结果与性能分析4.1 直观案例H0与H1下的分布差异先看一个直观的例子。M8、N128、SNR-6dB时取一次蒙特卡洛实验把H0和H1两种情况下融合中心收到的8个能量值分别画经验CDF再叠上理论噪声能量CDF差异非常明显。纯噪声时经验CDF基本贴着理论CDF走P-R指数大约在0.03~0.08之间有信号时能量值整体抬高经验CDF明显往右偏P-R指数跳到0.2以上。门限设在0.12左右two类情况分得很开。这也解释了为什么P-R检测器在低信噪比下还能有不错的性能——它不是看单个能量够不够大而是看整条曲线的偏移量。4.2 检测概率随SNR变化在 (P_f0.05)、蒙特卡洛5000次的配置下P-R指数检测器和等增益合并能量检测的检测概率随SNR的变化大致如下SNR (dB)P-R检测 (P_d)能量检测 (P_d)-120.230.18-100.450.39-80.680.62-60.850.80-40.950.93-20.990.9901.001.00低信噪比区间P-R检测比能量检测高3~6个百分点这个差距在对性能敏感的场合是很有价值的。高信噪比下两者都收敛到1没什么区别。能量检测在低SNR时吃亏的地方在于门限直接卡在噪声能量分布的尾部只要少量噪声样本超过门限就会虚警所以门限必须留足够裕量这反而压低了实际检测概率。P-R检测的判决依据是全局分布形态个别高能量噪声点对CDF面积的影响相对有限门限不用打那么大的裕量。4.3 关键参数的影响M、N与SNR的折衷M和N对性能的影响规律不一样仿真时可以分开扫。M增大时融合中心拿到的能量样本更多经验CDF更平滑P-R指数的方差明显下降。比如M从4涨到12在-8dB时检测概率大概能提升10个百分点以上。但M增大会线性增加控制信道开销和融合计算的输入规模实际系统中M不能无限加大。N增大时每个用户本地能量估计更准无论是H0还是H1能量值的离散度都更小两条分布之间的重叠区变窄P-R指数更容易区分。N的代价是感知时隙变长主用户对频谱的使用周期不等人所以N也不能拍脑袋取大。实际系统里M和N需要一起权衡。有一次测试我把M从8降到4同时把N从128提到256检测性能几乎没变但控制信道开销减半。这说明在调节系统参数时两个维度是可以互补的。如果硬要找一个经验规律总样本数 (M \times N) 大致不变时检测性能基本保持稳定但融合中心和用户之间的通信开销随M线性增长。所以优先保证N够用再把M压到满足协作分集要求的最小值是比较合理的调参思路。4.4 与能量检测的ROC对比在SNR-6dB时P-R检测的ROC曲线整体在能量检测的左上方。取几个典型点(P_f)P-R检测 (P_d)能量检测 (P_d)0.010.700.620.050.850.800.100.920.88虚警概率要求越严格P-R检测的优势越明显。这一点很有意思在 (P_f) 很低的工作区能量检测的门限被推到噪声分布的极右侧脆弱的尾部区域一旦噪声功率稍有波动门限就失效了。P-R检测的低虚警工作区相对稳定因为CDF面积度量对尾部单个异常点没那么敏感。5. 常见问题与调试心得5.1 能量归一化最容易被坑的地方写仿真第一版时我因为没把能量除以N理论卡方CDF和实际经验CDF完全对不上H0下P-R指数就有0.5以上门限完全失效。原因很简单如果不除以N能量值随N线性增长而理论CDF的表达式是按归一化能量 (E_i\frac{1}{N}\sum|y|^2) 推出来的。能量定义、理论CDF、卡方分布的自由度三者必须严格配套改任何一处都要同步检查另外两处。一个排查技巧先单独把H0数据喂进去如果算出来的P-R指数平均值明显大于0.1基本可以确定理论分布和观测分布没有对齐。正常纯噪声下经验CDF应该贴合理论CDFP-R指数应该在0.05以下。5.2 卡方自由度的匹配问题复基带信号对应自由度 (2N)实信号对应 (N)已经反复强调过。这里再提醒一个细节如果你的接收机模型里有上采样、匹配滤波等效采样点数要按实际进入能量累加的样本数来算不是按符号数。比如发了128个QPSK符号但每个符号过采样4倍那实际能量累加是512个点自由度按 (2 \times 512) 算。我在一次实验里就因为这个多算了一倍自由度导致理论CDF偏窄检测阈值被整体带偏。5.3 M太小时经验CDF锯齿严重怎么办M8时经验CDF只有8个台阶虽然能用但会有波动。如果仿真条件限制M只能在4左右有几个补救办法一个感知周期内让每个用户发送多次能量测量值融合中心把这 (M \times T) 个值合并起来一起构建经验CDF。这是对性能改善最明显的手段增加了样本量而不需要增加用户数。对能量值做核平滑处理把离散的阶跃CDF变平滑。这会引入核宽度的选择问题调试时要注意平滑后的CDF不能把真实分布拉开。放弃CDF路径改用直方图加PDF定义来算P-R指数bin数取 (\sqrt{M}) 量级。但正如前面说的直方图对bin选择敏感除非样本量实在太小否则我不太推荐。5.4 门限稳定性与蒙特卡洛次数P-R指数的门限完全依赖H0下的蒙特卡洛仿真。如果次数太少门限的随机波动会直接传导到检测概率上同一组参数跑两次结果差别很大你会误以为是算法问题。我的经验是做性能趋势扫描时蒙特卡洛至少2000次如果要对标论文里的精确数值至少5000次在低虚警率工作点(P_f0.01) 以下评估时要更多因为此时需要H0分布的极右侧尾部样本太少那个分位数根本不可靠。另外每次仿真前都要确认随机数种子。虽然Matlab默认每次启动都会重新设种子但如果你在一个大循环里反复调门限而没有用rng(0)固定对比结果就会带额外噪声。5.5 计算复杂度与代码提速技巧P-R检测器的计算量几乎可以忽略融合中心只需要对M个值排序、在约1000个网格点上计算理论CDF、做一次梯形积分。复杂度是 (O(G \log M M \log M))G是网格点数。在我的测试里5000次蒙特卡洛、7个SNR点总共大约几秒就跑完了主要的耗时反倒在生成随机信号上。如果要做更大规模的参数扫描有几个提速技巧computePRIndex里的理论CDF可以预先算好。当 (\sigma^2) 和 (N) 固定时(F_0(x)) 只取决于x完全可以在循环外生成一个查表向量循环内只查表。经验CDF那一步M很小时用arrayfun没问题如果M到了几千改成先排序再用histcounts或interp1速度能提一个量级。如果同一SNR下需要跑大量蒙特卡洛把接收数据一次性生成成一个大矩阵再向量化算能量比一层层for循环快很多。% 向量化生成M个用户、N个采样点的噪声 noise sqrt(sigma2/2) * (randn(M, N, numMC) 1i * randn(M, N, numMC)); E squeeze(mean(abs(noise).^2, 2)); % M x numMC这样一次就把所有蒙特卡洛实验的能量都算出来了后面循环只处理判决部分速度提升非常明显。5.6 一个容易忽略的细节信号功率归一化生成QPSK信号时如果直接用exp(1i * phase)每个符号的模都是1平均功率就是1。此时信号功率就是1乘以sqrt(Ps)才能得到想要的信噪比。但如果你用了成形滤波、插值或者信号是某个调制库生成的实际平均功率不一定等于1必须先测一下再归一化。更稳的写法是signal exp(1i * (randi([0 3], N, 1) * pi/2 pi/4)); signal signal / sqrt(mean(abs(signal).^2)) * sqrt(Ps);这样不管信号怎么生成最终功率都严格等于 (P_s)SNR的设定就不会跑偏。6. 后续还可以怎么扩展目前这套代码里我用的是AWGN信道、QPSK信号这是最干净的基础验证场景。实际系统里肯定会遇到衰落信道包括瑞利衰落和阴影衰落。扩展起来也不难给每个用户的信号前乘一个信道系数 (h_k)检测性能会整体掉一些但因为M个用户是空间分散的协作分集会把这些衰落平均掉一部分。另一个值得试的方向是噪声功率估计有偏置的情况。你把理论CDF里的 (\sigma^2) 和实际噪声真实功率故意设成不一致分别测一下P-R检测和能量检测的虚警概率会看到很有意思的现象这也是P-R检测相对能量检测最有说服力的优势场景。代码本身我已经整理好核心函数和主仿真脚本都在上面了参数也可以直接改。整个框架往下扩展成分布式协作、加权重融合、自适应门限都不会太费劲。实际做下来我的体会是P-R指数这种“比较两条曲线夹层面积”的思路比想象中要稳而且在真正干扰杂乱的场景里因为不依赖某个瞬时观测值抗波动能力确实比传统门限检测要好。如果你也在做协作感知的融合算法强烈建议把它加进对比列表里试一把。
