做认知无线电相关仿真的朋友应该都清楚频谱感知是整个系统的地基。我前阵子在做协作频谱感知项目时一开始用的是经典能量检测后来换了Pietra-RicciPR指数检测器做集中式数据融合效果比预期好了不少尤其是在信噪比低于-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 融合规则与门限确定融合中心用EGCT_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曲线我固定了这几个参数采样点数N1024感知用户数M4SNR-10dB蒙特卡洛次数numMc10000。这里SNR是指单个噪声样本信噪比用信号平均功率和噪声功率的比值来计算。跑出来的ROC曲线趋势很典型。PR检测器的曲线在低Pfa区域明显高于能量检测器。举个例子Pfa设为0.1时能量检测的Pd大概是0.62左右PR检测器能到0.78左右。差别主要来自PR对噪声不确定性的鲁棒性在仿真里我故意给噪声功率加了正负2dB的随机扰动能量检测的性能立刻掉下来PR的损失小很多。检测器Pfa0.01时的PdPfa0.1时的Pd能量检测0.310.62PR指数检测0.450.78这个结果符合我预期也是我最后把PR写到项目结论里的依据。表格里的具体数值依赖随机种子但相对关系和趋势是稳定的。4.2 检测概率随SNR的变化再看检测概率随SNR变化的曲线。横轴SNR从-20dB梯度到0dB纵轴是检测概率固定Pfa0.1。低信噪比区间也就是-15dB到-8dB这段PR的优势最明显。在-12dB附近PR检测器的Pd能到0.5能量检测只有0.32左右。到了-5dB以上两者都接近1差距反而看不出来了这时候说明信道条件足够好选哪种检测器差别不大。这个现象很好理解。高信噪比下信号能量远大于噪声任何基于能量的统计量都能分开H0和H1。低信噪比下统计量的高阶特性才会体现差异。4.3 感知用户数对融合增益的影响我还跑了一组用户数变化的实验。M分别取1、2、4、8、16SNR固定-12dB其他参数不变。结果如下用户数M单用户PREGC融合PR10.180.1820.180.2640.180.4180.180.55160.180.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严重偏离预设门限集与测试集混用单独生成门限标定样本统计量出现NaNm2过小导致除零分母加eps或归一化跑完ROC曲线抖动剧烈蒙特卡洛次数不足提到10000次以上多用户融合无增益所有用户信道完全相关检查信道系数是否独立生成高SNR下性能不升反降数值溢出或信号削波样本归一化检查幅值范围parfor不生效变量切片读写混乱用随机数流隔离或改用普通for6. 后续扩展方向PR指数检测器目前用的是全局二阶矩和四阶矩没有利用信号本身的循环平稳特性。如果主用户信号是OFDM或者单载波循环前缀这类带周期相关性的信号可以引入循环谱域的处理思路把PR指数扩展到频域子带上理论上能进一步对抗窄带干扰。融合规则这块也可以继续挖。我之前拿EGC做了基线后面可以试最大比合并或者选择性合并也就是选择统计量最优的K个用户参与融合牺牲少量性能换取信令开销的大幅下降。还有一点是关于门限标定。蒙特卡洛标定虽然准确但是在线系统里很难实时跑出大量H0样本。工程落地时可以考虑用解析近似式或者查找表来替代蒙特卡洛把门限计算时间从秒级压到毫秒级。这个方向我后面准备专门用一篇来做对比分析。我在实际跑这个项目时有一个比较深的体会检测器性能再强也架不住融合结构设计不合理。好的融合规则能成倍放大检测器的优势反过来如果用户选择策略粗糙所有统计量一股脑全送到融合中心反而可能被深衰落用户拖累。仿真里可以加一个简单的信噪比筛选门低于某个阈值就直接弃用这个用户的感知数据效果会比你预想的还要明显。再分享一个写Matlab仿真的小习惯。所有随机数生成之前先把随机种子固定下来rng(2024);这样每次跑出来的结果完全可复现改一个参数重跑时能精确对比到底是谁影响了性能。对发论文、写报告、复现同学的结果都特别关键。我见过太多跑出来的结果自己都复现不了的案例根子就在随机种子没控制住。这套代码整体跑一遍从参数配置到出图半小时内能完成。如果你在调参或者复现的过程中有卡壳的地方优先检查噪声功率换算和门限标定这两块90%的异常结果最后都能追溯到这两处。