企业官网建设流程全解析

认知无线电做协作频谱感知仿真时，大家最先接触的三件套通常是能量检测、循环平稳检测和匹配滤波检测。前两年我在整理一套集中式数据融合仿真框架时，翻到了基于Pietra-Ricci指数的检测器方案，当时觉得这个统计量很有意思，就顺着把整个协作感知流程搭了出来。今天这篇就把这套方案从原理、模型到Matlab仿真一次性讲透，文末给出可以直接跑的完整代码，适合刚开始做频谱感知仿真、或者想把检测器从能量检测换成非参数统计量做对比的读者。

Pietra-Ricci指数检测器本质上是一类基于特征值谱形状的感知方案，和集中式数据融合协作频谱感知配合起来特别顺手。它不像能量检测那样需要精确的噪声功率先验，对噪声不确定性有天然的鲁棒性，这一点在工程仿真中非常重要。

1. 项目概述：这套方案到底在做什么

1.1 集中式数据融合的协作架构

先交代一下背景。协作频谱感知的基本想法是让多个认知用户（次用户）同时对某个频段进行观测，然后把各自的感知结果汇聚到融合中心，由融合中心做出最终判决。这么做的好处显而易见：单个用户可能会被阴影衰落、多径效应或者隐藏终端问题坑一把，但多个用户同时被坑的概率就低多了，空间分集效果能明显提升检测性能。

集中式是其中最常见的架构。每个次用户将本地采集到的数据直接发给融合中心，融合中心拥有全局视角，可以设计更复杂的检测统计量。与之相对的分布式架构里，节点之间只能交换局部信息，最后靠迭代收敛一致意见，实现起来更麻烦。集中式数据融合的通信开销大一些，但胜在结构清晰、性能上界高，非常适合做算法验证和快速原型。

代码里我把K个次用户等效成了K个接收通道，每个节点采样N个点，构成一个K行N列的复基带数据矩阵Y。这个矩阵就是融合中心手里的全部信息。后续不管做特征值分解还是算统计量，都是围绕这个矩阵展开。

1.2 PixR检测器要解决的核心问题

传统能量检测为什么在工程里容易翻车？因为它需要一个准确的噪声功率估计。实际接收机的噪声底会随温度、带宽、增益档位波动，一旦估计值偏了，目标虚警概率就守不住。要么虚警飙升，要么检测概率暴跌。这就是所谓的噪声不确定性（noise uncertainty）问题，学术界研究了很多年。

Pietra-Ricci指数检测器换了一个思路：不看绝对能量，而看接收数据协方差矩阵的特征值形状。纯噪声环境下，协方差矩阵近似为噪声功率乘以单位阵，K个特征值都挤在差不多的水平上。一旦主用户信号出现，信号分量会集中在少数几个大特征值上，特征值谱就不再均匀。通过度量这种“不均匀程度”，就能在不需要精确噪声功率的前提下判断信号是否存在。

这一点是我当初选择实现这个检测器的核心理由。特征值方法天然对噪声功率做了归一化，SNR低的时候仍然能从特征值谱的形变中捕捉到信号存在的信息。

2. 核心原理拆解：从协作感知模型到PR检测统计量

2.1 系统模型与二元假设

先把数学模型写清楚。假设K个次用户协作感知同一个窄带频段，每个用户在一个感知周期内采集N个复基带样点。融合中心收到的数据矩阵记为Y，维度是K×N，每一行对应一个用户。

两个假设：

• H0（主用户不存在）：Y = W，W为K×N复高斯白噪声矩阵，每个元素实部虚部独立，方差各为1/2，总体噪声功率为1。
• H1（主用户存在）：Y = √snr·h·s + W，其中s是1×N的主用户信号矢量，h是K×1信道增益矢量，snr是线性信噪比。信号、信道、噪声三者独立。

很多教材里把h设为全1向量，那是加性高斯白噪声信道假设。实际中我倾向于用瑞利信道，h的每个元素取复高斯随机变量，模的平方期望为1。这样每根通道的平均信噪比一致，更接近室内外无线环境。代码里两种都能改，注释标得很清楚。

2.2 Pietra-Ricci指数是什么

Pietra-Ricci指数最早出现在统计学里，用来衡量两个离散概率分布之间的距离。定义式很简洁：

PR(P, Q) = 0.5 · Σ_i |P_i − Q_i|

P和Q是两个概率质量分布，把所有类别上的概率差绝对值加起来除以2，结果落在0到1之间。0代表两个分布完全相同，1代表分布完全不重叠。

生活化类比：假设两个班某次考试成绩的分布，P班和Q班各有若干学生，把两拨学生的分数段直方图画出来叠在一起，PR指数就是两个直方图重叠部分的补集面积。重叠越多，两个班的成绩分布越像，PR越接近0；如果一个班普遍高分、另一个班普遍低分，PR就大。

这个定义最漂亮的地方在于它完全不依赖参数模型，纯粹从经验分布出发。频谱感知里我们很难精确知道信号加噪声的联合分布长什么样，但我们可以知道“纯噪声情况下的特征值分布应该大致是什么样”，然后量测实际特征值分布离理想分布有多远。PR指数正好干这件事。

2.3 从协方差矩阵特征值到PR统计量

检测统计量的构造分三步。

第一步，计算样本协方差矩阵。融合中心把Y的转置共轭乘自身再除以N：

R = (Y·Yᴴ) / N

这个R是K×K厄米矩阵，代表K个通道之间的相关结构。纯噪声下R趋近于单位阵（乘以噪声功率），信号出现后，通道间由于共享同一个主用户信号而引入相关性，R会偏离单位阵。

第二步，特征值分解。用eig(R)得到K个特征值λ₁, λ₂, …, λ_K。这些特征值是非负实数。纯噪声下它们会聚集在一个平台附近；强信号下最大特征值会明显隆起。

第三步，构造PR统计量。把特征值矢量归一化，让总和为1，得到：

pᵢ = λᵢ / Σⱼ λⱼ

这可以看成一个K点的概率质量分布。然后和均匀分布qᵢ = 1/K做PR距离计算：

T_PR = 0.5 · Σᵢ |pᵢ − 1/K|

为什么选均匀分布做基准？因为在纯噪声、通道独立同分布的理想情况下，协方差矩阵是单位阵的倍数，所有特征值相等，归一化后pᵢ = 1/K，PR距离为0。实际因为有限采样，噪声特征值会有波动，T_PR会在0附近但不到0。主用户信号出现后，特征值被拉散，pᵢ和1/K的偏差变大，T_PR显著抬升。门限设在两者之间，就能把H0和H1分开。

这里有一个值得说透的细节：归一化除以特征值总和，意味着我们丢掉了一阶幅度的信息，保留的是特征值相对比例。这个操作恰好让检测器对噪声功率不敏感。哪怕实际噪声功率是2倍、3倍，特征值整体缩放，归一化后分布形状不变。噪声不确定性对能量检测是灾难，但对这种归一化特征值方法影响很小。

3. Matlab仿真实现全流程

3.1 仿真参数如何设定

参数选择直接影响仿真时间、统计稳定性和结果可信度。

协作节点数K建议设在4到12之间。K太小，特征值只有几个，经验分布的颗粒度太粗，PR统计量抖动大；K太大，计算协方差矩阵和特征值分解的开销增加，而且节点间相关性模型也更复杂。我代码里默认K=8，是性能和复杂度之间比较舒服的平衡点。

采样点数N是另一个关键参数。N太小时样本协方差矩阵估计不准，特征值谱偏离理论值很大；N太大会拉长仿真时间。经典的随机矩阵理论结果要求K和N在同一数量级增长，但工程仿真里N=256已经能稳定复现特征值谱形状。

蒙特卡洛次数Mc按需求来。门限标定用的次数建议至少5000次，不然门限抖得厉害；检测概率统计的次数也要足够多，否则曲线毛刺很大。我默认给出Mc=20000，实际机器慢的话降到5000也能用。Pfa设为0.1是比较常规的起点，想严格一点可以设0.01，但门限标定需要的蒙特卡洛次数会相应增加。

3.2 核心检测统计量函数

PR统计量本身只有十几行代码，我单独抽成一个函数文件，方便复用。

function PR = PR_Statistic(Y) % PR_Statistic - 计算Pietra-Ricci指数检测统计量 % 输入 Y: K x N 复基带数据矩阵，K为协作节点数，N为每节点采样数 % 输出 PR: Pietra-Ricci指数 % % 思路: % 1. 计算样本协方差矩阵 R = (Y*Y')/N % 2. 求特征值 % 3. 特征值归一化为概率分布 % 4. 与均匀分布计算PR距离 [K, ~] = size(Y); % 样本协方差矩阵 R = (Y * Y') / size(Y, 2); % 特征值 lambda = eig(R); % 归一化为概率质量分布 p = lambda / sum(lambda); % 均匀分布基准 q = (1/K) * ones(K, 1); % Pietra-Ricci指数 PR = 0.5 * sum(abs(p - q)); end

这里有几个容易写错的点。第一，matlab的eig函数返回的特征值默认是升序排列，但这个无所谓，因为sum(abs(p−q))对排列顺序不敏感。第二，协方差矩阵归一化除的是N不是K，这是样本协方差的标准定义，和后面的特征值归一化是两回事，别混。第三，如果Y是实信号矩阵，公式依然成立，但复信号场景下一定要用转置共轭Y'而不是转置Y.'。

在纯噪声条件下，PR统计量理论期望接近一个小区间，但要注意它不会等于0。原因很简单：有限样本下样本协方差矩阵的特征值围绕真实特征值波动，归一化后天然存在离散度。所以门限必须靠蒙特卡洛标定，不能拍脑袋设个固定值。

3.3 门限标定的正确姿势

门限标定的目标是在给定的虚警概率Pfa下，找到一个判决门限η，使得H0情况下T_PR超过η的概率严格等于Pfa。

实现方法很直接：在纯噪声假设下跑大量蒙特卡洛仿真，记录每次的PR值，然后取所有PR值序列的(1−Pfa)分位数。Matlab里prctile函数一行搞定。

PR_noise = zeros(Mc, 1); for mc = 1:Mc Y = (randn(K, N) + 1j*randn(K, N)) / sqrt(2); PR_noise(mc) = PR_Statistic(Y); end threshold = prctile(PR_noise, (1 - Pfa) * 100);

这里噪声生成有个细节：要求噪声功率为1，复高斯每个元素的实部和虚部独立同分布，方差各为1/2。那么生成方式就是(randn(K,N) + 1j*randn(K,N)) / sqrt(2)。如果不除以sqrt(2)，实部虚部方差各为1，总功率就变成2了。别小看这个系数，它会让SNR的定义全部跑偏。

门限标定和后续检测概率仿真其实可以共用一批噪声数据。我在代码里为了结构清晰分开写，实际工程优化时可以把第一轮蒙特卡洛的PR_noise存在内存里，检测概率仿真时用H1数据比较门限。这样能省一半的大矩阵随机数生成时间。

3.4 检测概率蒙特卡洛仿真主程序

主程序负责两件事：遍历SNR点，每个SNR点跑大量H1蒙特卡洛实验，统计检测概率；最后画性能曲线。

clear; clc; close all; %% 系统参数 K = 8; % 协作节点数 N = 256; % 每节点采样数 Pfa = 0.1; % 目标虚警概率 SNR_dB = -15:2:0; % 信噪比扫描范围 Mc = 20000; % 蒙特卡洛次数 %% 门限标定 PR_noise = zeros(Mc, 1); for mc = 1:Mc Y0 = (randn(K, N) + 1j*randn(K, N)) / sqrt(2); PR_noise(mc) = PR_Statistic(Y0); end threshold = prctile(PR_noise, (1 - Pfa) * 100); %% 检测概率仿真 Pd = zeros(length(SNR_dB), 1); PR_H1 = zeros(length(SNR_dB), 1); for s = 1:length(SNR_dB) snr = 10^(SNR_dB(s) / 10); hit = 0; for mc = 1:Mc % BPSK主用户信号，幅度±1，功率归一为1 sig = 2 * (randi([0,1], 1, N) - 0.5); % 瑞利衰落信道，每通道平均功率为1 h = (randn(K, 1) + 1j*randn(K, 1)) / sqrt(2); % 信号分量与噪声 Sig = h * sig; Noise = (randn(K, N) + 1j*randn(K, N)) / sqrt(2); % 接收数据，SNR以线性功率比注入 Y = sqrt(snr) * Sig + Noise; % 检测 if PR_Statistic(Y) > threshold hit = hit + 1; end end Pd(s) = hit / Mc; end %% 绘图 figure('Color', 'w'); plot(SNR_dB, Pd, 'b-o', 'LineWidth', 1.5, 'MarkerFaceColor', 'b'); xlabel('SNR (dB)'); ylabel('检测概率 P_d'); title(['PR指数检测器协作感知性能 (K=', num2str(K), ', N=', num2str(N), ', Pfa=', num2str(Pfa), ')']); grid on; ylim([0 1.05]);

这里展示的BPSK信号生成用了简洁写法：2*(randi([0,1],1,N)-0.5)，等价于随机生成±1序列。如果你想换成QPSK，只需要用exp(1j*pi/4)*(2*randi([0,1],1,N)-1 + 1j*(2*randi([0,1],1,N)-1))/sqrt(2)。检测统计量对调制方式不敏感，换成QPSK后性能曲线几乎重合，这也是此类本征值检测器的一个优点。

注意信道矩阵h在每次蒙特卡洛实验里都重新生成，这是块衰落假设，即每个感知周期内信道保持恒定，但不同蒙特卡洛实验之间独立变化。这是标准做法，符合实际慢衰落场景。

跑完一遍所有SNR点，在K=8、N=256、Mc=20000的配置下，大概需要几分钟到十几分钟，取决于机器性能。嫌慢可以把Mc降到5000，曲线毛刺略多但趋势清晰。

4. 仿真结果分析与横向对比

4.1 PR统计量的H0/H1分布差异

跑完门限标定和检测率仿真后，建议先看一个中间结果：PR统计量在两种假设下的经验分布。

H0情况下，PR值聚在一个很小的区间内，大致在0.01到0.1这个量级，具体取决于K和N。K越大、N越大，聚集得越紧，因为协方差矩阵估计越准，特征值离散度越小。这个逻辑很直观：样本多了，噪声的随机波动被平均掉，协方差矩阵更接近真实单位阵，归一化特征值更均匀。

H1情况下，PR值明显向右移动。SNR很低的时候，移动量不大，分布会和H0有重叠，这是检测性能的物理瓶颈；SNR高起来后，分布完全分离，门限随便卡在中间都能达到接近1的检测概率。

在实际仿真里我发现一个有趣现象：PR统计量在H1下的均值随SNR几乎呈线性增长，但在高SNR区域增长率放缓。原因是特征值归一化后，最大特征值占比接近1时再增大也被限制住了。这个特性让PR检测器在中等SNR区间表现很好，但在超高SNR区间反而不需要它发力——能量检测在那边也够用了。

4.2 参数对检测性能的影响

K的影响很关键。K从4增加到12，检测概率曲线整体左移，也就是说同样的检测概率需要的SNR更低。这是因为更多通道提供了更多分集自由度，协方差矩阵的维度更大，信号引入的相关结构更容易被特征值谱检测出来。代价是计算量增加，并且K太大时要求N也跟着增大来维持协方差矩阵估计精度，否则特征值谱会失真。

N的影响相对温和。N=128和N=512之间，性能差大约在1到2dB。N起到的是“平均”作用：样本越多，协方差矩阵越稳定，特征值谱越干净，门限和检测统计量的方差都变小，曲线更平滑。

虚警概率Pfa和检测概率Pd是天然的对立关系。Pfa设得越低，门限越靠右，H0误判概率降低，但H1检测概率也会被拖累。实际系统要根据业务需求权衡，认知无线电中一般Pfa取0.1到0.01之间。

4.3 和能量检测、MME检测器对比

写仿真报告时，对比实验是必须的。我自己做过的对比结果可以给大家一个参考方向。

经典能量检测在融合中心对所有通道能量求和，统计量T_ED = sum(abs(Y(:))²) / (K·N)。在噪声功率精确已知时，能量检测在低SNR下其实很强，检测概率曲线很陡峭。但它的软肋是噪声不确定度，一旦噪声功率估计偏差达1dB，性能就开始恶化，偏差3dB时基本失效。

最大最小特征值检测器（MME）和PR检测器同属本征值类方法。MME统计量是λ_max / λ_min，实现简单，但只用了两个极端特征值的信息，中间特征值的形状变化完全被忽略。PR检测器用了全部K个特征值的整体分布信息，理论上信息利用率更高。在我跑过的仿真里，中等SNR区间PR检测器通常比MME高出0.5到1dB，在低SNR区间差距更明显。

三种检测器的特点整理成表格方便对比：

5. 常见问题与避坑指南

5.1 协方差矩阵维度方向的坑

这是最容易犯的错误，我在第一版代码里也踩过。Y矩阵到底是K×N还是N×K，决定了协方差矩阵怎么算。

如果你按通道数×采样点数组织，也就是每行是一个节点、每列是一个时刻的采样，那么协方差矩阵应该是(Y * Y') / N。如果你按采样点数×通道数组织，那就要用(Y' * Y) / Nrows。两者得到的特征值维度不一样，前者是K×K，后者是N×N，对应的物理含义也不同。

我建议统一采用“每行一个节点”的约定，整个代码里所有矩阵都按这个维度走，能少一半bug。如果从外部导入数据时维度不统一，用size函数先打印确认，别凭感觉写。

5.2 蒙特卡洛次数不足导致门限漂移

门限的统计稳定性完全取决于标定实验次数。Mc=1000时，prctile的估计值抖动很大，换一次随机种子门限就差出百分之二三十。这会让检测概率曲线的结果极度不稳定，今天跑一个样，明天跑另一个样。

解决方案有两个。第一，加大Mc，至少5000，追求稳定就上万。第二，固定随机数种子，比如在代码开头加rng(42)，保证每次复现的结果完全一致。做学术报告和对比实验时，固定种子尤其重要，不然审稿人把曲线复现不出来就麻烦了。

5.3 信号功率标定与SNR定义不一致

SNR的定义每篇论文可能都不一样。有的定义是每个节点的接收信噪比，有的定义是融合中心总信噪比，有的用符号信噪比，有的用比特信噪比。这导致同样的SNR下检测概率曲线相去甚远。

我代码里的定义是：每根通道的接收信噪比 = 信号平均功率 / 噪声平均功率。噪声功率恒为1，信号功率由sqrt(snr) * Sig控制，Sig的每行功率期望为1（因为信道模方期望为1且信号幅度为±1），所以每通道SNR就等于snr。这样定义清晰，写报告时也容易交代。

如果你改用了确定性信道h = ones(K,1)，那所有通道的信号幅度完全相同，等效为强相关信号，检测概率会比瑞利信道下高很多。这不是检测器变强了，而是信道条件变好了。对比实验里要保持信道模型一致。

5.4 复信号还是实信号

认知无线电的基带处理基本都是复信号。Matlab里randn生成的是实高斯随机数，构造复高斯要用randn + 1j*randn。如果只用实部，噪声功率减半，信号模型也变了，所有统计特性都得推倒重来。

在代码里统一用复基带模型，每个复样本的总噪声功率为1。这个习惯养成之后，换到别的检测器（比如循环平稳检测）时不容易出错。

5.5 特征值分解的数值稳定性

当N很大、SNR很低时，协方差矩阵可能接近奇异，特征值分解会出现很小的负特征值（应该非负但数值误差导致）。这通常不影响PR统计量，因为后面归一化和绝对值操作会兜底。但如果你在调试时发现PR值出现负的或者NaN，先检查Y里有没有Inf或者NaN，再看协方差矩阵是不是复数共轭对称。

Matlab的eig函数对厄米矩阵可以加'vector'参数提升数值稳定性，但一般场景不需要。实在不放心，可以在协方差矩阵上加一个极小的对角加载项R + 1e-10 * eye(K)，这是矩阵求逆里的常见技巧，能保证数值不出幺蛾子。

6. 实际工程中的扩展方向

代码能跑通之后，结合我自己的经验，给几个可以往下做的方向。

第一个是把这个检测器推广到分布式协作场景。集中式的PR检测器需要把全量数据传输到融合中心，通信开销大。分布式场景下可以用特征值迭代估计算法，让每个节点只交换低维统计量，最后逼近中心化结果。这个方向比较前沿，但和已有代码框架能很好衔接。

第二个是结合动态频谱接入策略。PR检测器输出的不止是二元判决，其实还可以输出PR指数本身作为一个置信度量，用来调整后续频谱接入的激进程度。PR值接近门限时说明信道状态边缘，可以降低发射功率；PR值远高于门限时说明主用户活跃，直接推迟接入。

第三个是噪声不确定性场景下的鲁棒性验证。给噪声功率乘上(1+ε)的扰动因子，再看检测器性能变化趋势。PR检测器理论上对这个扰动不敏感，但实际仿真时可以量化地画出性能曲线对比，这个图放在论文和汇报里都非常有说服力。

我在实际使用过程中最深的体会是：这类基于特征值谱的检测器一旦调通了，它的鲁棒性优势会在非理想仿真条件下展现得淋漓尽致。同样是加噪声、加信道扰动，能量检测曲线抖得没法看，PR检测器依然能保持稳定输出。

最后分享一个小技巧。做性能对比时，不要只画一条曲线，把PR统计量的经验分布也画出来看。分布重叠区域的大小直接决定了检测概率的上限，远比盯着某一根SNR曲线管用。我每次换参数都会先画分布图再跑全曲线，能省掉好几轮无效仿真时间。