简介面向蜂窝通信与随机几何研究者的PPP仿真代码包聚焦泊松点过程PPP与蒙特卡洛方法在基站部署、干扰及覆盖分析中的应用。资源包含完整的MATLAB仿真脚本可用于生成PPP分布、执行随机试验并评估不同基站密度下的网络性能指标适合通信专业学生、科研人员及网络规划工程师用于算法验证与方案对比。压缩包共33个文件以28个m脚本为核心辅以3个fig图表与2张jpg结果图整体仅135KB便于快速下载与二次开发。已有347人学习使用。代码覆盖多种典型场景如蜂窝重叠覆盖分析、依赖型PPP建模、用户密度影响、频谱效率与干扰变化等并附有结果图像与绘制脚本可直接对照输出进行复现。通过研究这些文件可深入理解PPP随机几何与蒙特卡洛仿真相结合的方法并据此优化自身网络建模流程。1. PPP 蒙特卡洛蜂窝仿真这套脚本到底能帮你算什么做蜂窝网络性能分析的人迟早会遇到一个反直觉的坎基站位置根本不能假设成均匀网格。真实网络里宏站、微站、室分混在一起PPP泊松点过程随机几何把基站看成平面上随机撒的点再用蒙特卡洛跑几千次取平均覆盖概率、干扰、ASE 这些指标才挨得上真实的网络。这套 MATLAB 脚本正是干这个的。它适合三类人通信方向做毕设的研究生、刚摸到随机几何门槛想跑通一条链路的入门者以及需要给论文补仿真图的工程师。它能解决的是“怎么把 PPP 和蒙特卡洛落到可复现的代码上”这个具体问题而不是给你一个炫酷的黑匣子——里面的坑和细节本文都会拆开讲。2. 文件清单与角色分工33 个文件里谁是入口、谁是绘图、谁是工具解压后会看到 .m 脚本和 .fig 图形文件混在一起。这里要提醒一句这套资源不是按“一个 main 函数带一堆子函数”的标准工程结构组织的更像是一堆互相引用工作区变量的脚本集合。直接双击 run 大概率会报错正确做法是先搞清楚文件角色再决定执行顺序。2.1 文件分类先分清主仿真、ASE 分析、绘图和辅助工具我把压缩包里的 33 个文件按用途拆成五类对照着看就不会迷路。类别文件作用主仿真入口PPP_dep_sim.m、PPP_Intra_dep.m、test.m、example.m生成 PPP 布点、跑蒙特卡洛主循环、算基础指标ASE/干扰分析ASE_theta_5dB_lambda_u.m、intra_dep_ASE_.m、inter_dep_ASE_theta_5dB.m算面积频谱效率扫用户密度或 SIR 阈值绘图脚本fig6.m、fig7_mu.m、fig8.m、fig8_pu.m、fig9.m、fig12_a.m、fig12_b.m、fig13_c.m、fig13_c1.m、fig14.m、fig15.m、fig16.m、fig17.m复现论文里的性能曲线图形结果fig14.jpg、fig14.fig、fig14_zoom.fig、fig15_zoom.fig、fig14_zoom.jpg已导出的成品图、可继续编辑的 fig 图、局部放大细节工具函数count.m、multiply.m、btitle.m、brainavi.m、overflap.m、overflap1.m、picture.m、picture_intra_dep.m计数、矩阵乘法、图标题、动画/绘图辅助从文件命名能看出一个关键信息这套代码研究的是基站和用户之间的“空间依赖”dependence。PPP_dep_sim.m 和 PPP_Intra_dep.m 分别对应两套依赖模型intra_dep 可以理解为小区内用户与基站的相关性inter_dep 是小区间干扰的依赖结构。如果你做的是异构网络或者 D2D 场景这两套脚本可以直接改造成不同层之间的干扰模型。fig14.jpg 和 fig14_zoom.jpg 是作者最终的论文出图fig14.m 是重绘脚本。我建议把 fig14.m 的输出当作整套代码的“验收标准”——你的链路跑通后画出来的曲线应该和这张 jpg 形态一致。如果对不上优先怀疑参数单位或边界处理而不是代码逻辑。2.2 从文件名反推仿真链路theta_5dB 和 lambda_u 是什么意思这套资源的文件名信息量其实很大。ASE_theta_5dB_lambda_u.m 拆开读ASE 是面积频谱效率单位是 bps/Hz/km²theta_5dB 是 SIR 判决阈值 5 dBlambda_u 是用户密度。所以这个脚本干的活是固定 SIR 阈值 5 dB扫描不同的用户密度画出 ASE 变化曲线。这是蜂窝网络随机几何里最经典的一张参数扫描图。inter_dep_ASE_theta_5dB.m 里的 inter_dep 说明它对比的是考虑基站与用户空间相关性前后的 ASE。讲白了就是“用户独立均匀撒点”和“用户跟着基站热点聚在一起”两种极端模型下ASE 差多少。这种对比在论文里通常呈现为三条曲线无依赖模型、有依赖模型、解析解。fig7_mu.m 里的 mu 大概率是对数正态阴影衰落的均值参数fig8_pu.m 里的 pu 可能是用户占比或功率相关参数。这类命名虽然不规范但足够说明作者在做参数敏感性分析每个 fig 脚本固定其他参数只扫一个关键变量。这套资源的正确打开方式不是“找入口”而是“找图号”。你想复现哪张图就找到对应 fig 脚本看它引用了哪些变量再回溯到主仿真脚本把这些变量算出来。第三、四章我会按这条链路把 PPP 生成、蒙特卡洛主循环、SIR 到 ASE 的计算讲清楚。提示fig12_a.m 和 fig12_b.m 是成对出现的通常一个扫基站密度、一个扫用户密度对比同一指标。改参数时留意这些脚本内部有没有重复跑蒙特卡洛循环如果有先把它内部的迭代次数调小确认逻辑再放大量级否则一个 fig 脚本可能跑半小时。3. 从 PPP 生成到蒙特卡洛框架两块地基代码怎么搭压缩包里没有单独的 generatePPP.m 函数常见的做法是把 PPP 生成直接写进主循环。我习惯单独抽一个函数出来调试和换密度都方便。PPP 的数学定义不复杂在区域 A 内点的数量服从泊松分布且这些点在 A 内均匀独立分布。要生成它只需要两步先抽个数再撒点。3.1 为什么蜂窝仿真非得用泊松点过程先放结论PPP 是“最随机”的空间点过程它不假设任何空间相关性。真实基站部署受地形、人口密度、频谱管理影响但随机几何里第一条路就是把基站当 PPP 来算因为这样干扰的拉普拉斯变换有闭式解覆盖概率能写出解析式。对每平方公里 10~50 个宏站的密度范围PPP 模型和实测覆盖概率能对上这在大量文献里反复验证过。相对网格部署PPP 的优点是理论可导。单层蜂窝网络、用户关联最近基站的场景下SIR 覆盖概率可以写成Pc(θ) 1 / (1 (2/α) · θ^(2/α) · ∫_{θ^(-2/α)}^{∞} 1/(1 u^(α/2)) du)这个解析式是判断蒙特卡洛代码对不对的标尺。路径损耗指数 α4 时覆盖概率随阈值 θ 的衰减能被这个式子精确预测。你的仿真结果和它差超过 5%基本可以断定是单位换算或边界处理出了问题不是在算法层面纠结的时候。生成 PPP 前先确定强度 λ单位是个每平方米。基站密度通常按“个每平方公里”给换算时要除以 1e6这个换算最容易在蒙特卡洛循环里埋雷。标准生成做法是先按泊松分布抽样点的个数再在仿真区域内均匀撒点。function [pts, n] generatePPP(L, lambda_km2) % L: 仿真区域边长(米) % lambda_km2: 基站密度(个/平方公里) lambda lambda_km2 / 1e6; % 换算成 个/平方米 n poissrnd(lambda * L^2); % 点的数量服从泊松分布 pts L * rand(n, 2); % 在 L×L 区域内均匀撒点 end这段代码有三个参数要盯紧。L 建议取 1000 到 2000 米太小了边界效应明显太大了单次 realization 的撒点数过多内存和循环时间都上来了。lambda_km2 是密度名义值注意 1e6 的换算系数别丢。poissrnd 需要统计工具箱如果没有可以用逆变换法近似抽样但收敛慢不推荐在蒙特卡洛循环里用。提示poissrnd 在部分工具箱版本里不支持向量化输出几千次循环调用会有明显开销。我一般会先把所有 realization 需要的总数一次性算出来再按批切分省掉重复调用的损耗。3.2 蒙特卡洛主循环每个 realization 都要重新撒点蒙特卡洛的核心是“大量重复独立实验取平均”。蜂窝仿真的每一次 realization 必须重新生成基站和用户位置不能用同一组点重复算那样会低估空间相关性带来的方差。主循环的标准骨架长这样N_real 2000; % 实现次数, 论文里常见 1000~5000 theta_db 5; % SIR 阈值 5 dB theta 10^(theta_db/10); % 阈值转线性值 alpha 4; % 路径损耗指数, 城区一般取 3~4.5 P_tx 1; % 发射功率归一化 N0 1e-9; % 噪声功率, 量纲必须与信号一致 L 1000; % 仿真区域边长(米) lambda_b 20; % 基站密度 个/km² cov_count 0; % 统计超过阈值的用户数 for iter 1:N_real bs generatePPP(L, lambda_b); n_ue 100; % 每个格子固定撒 100 个用户 ue L * rand(n_ue, 2); for i 1:n_ue d2 sum((bs - ue(i,:)).^2, 2); % 到所有基站的距离平方 [dmin2, idx] min(d2); % 最近基站作为服务基站 sig P_tx * dmin2^(-alpha/2); % 有用信号功率 d2(idx) inf; % 排除服务基站 intr sum(P_tx * d2.^(-alpha/2)); % 累加所有其他基站的干扰 sir sig / (intr N0); % 含噪声的 SINR if sir theta cov_count cov_count 1; end end end Pc cov_count / (N_real * n_ue); % 覆盖概率估计几个细节值得展开。d2 存的是距离平方所以算功率时是 d2.^(-alpha/2)等价于 d^(-alpha)对应路径损耗模型 l(r)r^(-α)。最近距离用 min 找这是假设用户关联最近基站——现实中是最大接收功率关联在发射功率一致时两者等价但如果做异构网这里要改成按最大 SINR 或最大功率选站。干扰累加前把服务基站的距离置为 inf这个“排除自身”的步骤漏掉的话SIR 会虚高得离谱。噪声功率 N0 是最容易翻车的量。论文里常常假设热噪声功率谱密度 -174 dBm/Hz配 10 MHz 带宽折算是 -114 dBm然后还要跟信号功率量纲对齐。如果发射功率用 1 W噪声就得用 W 或 mW 换算混用单位会让覆盖概率出现“全 0 或全 1”两个极端。我建议先把量纲定成 mW把所有功率值统一换算后再进公式。最后一个坑是边界效应。L×L 正方形区域边缘的用户周围基站天然比中心少干扰被低估。常规解法是“保护带”撒点时用 (L2g)×(L2g) 的区域统计时只统计中心 L×L 里的用户g 取 200~500 米。这样边缘用户能看到完整的干扰环境覆盖概率曲线才平。第 5 章会再展开讲具体参数怎么设。4. 把 SIR 算成覆盖概率和 ASE复现 fig14 和 ASE 曲线的完整链路有了 PPP 生成和蒙特卡洛主循环下一步就是把 SIR 结果折算成论文里常见的覆盖概率曲线和 ASE 曲线。这是从“能跑”到“能出图”的关键一跳。很多新手卡在这一步不是仿真逻辑错而是统计口径和阈值换算出了问题。4.1 SIR 阈值、覆盖概率和解析对照覆盖概率的定义很直接随机抽取一个用户它的 SIR 大于某个阈值的概率。在蒙特卡洛里就是超过阈值的样本数除以总样本数。阈值 θ 用 dB 给时进比较前必须先转线性值即 10^(theta_db/10)。这个换算在第三章的骨架代码里已经出现过但值得单独拎出来强调5 dB 是 3.16 倍不是 5 倍。拿线性值 5 当阈值用覆盖概率会明显偏高和解析解对不上。我一般会把仿真结果画成“覆盖概率 vs SIR 阈值dB”的曲线横轴从 -10 dB 扫到 20 dB每次一个阈值跑一遍蒙特卡洛。每个阈值都是独立仿真不是在一组 SIR 样本上套不同阈值——后者虽然省时间但会引入相关性曲线比真实情况平滑。theta_db_list -10:2:20; Pc_list zeros(size(theta_db_list)); for t 1:length(theta_db_list) theta 10^(theta_db_list(t)/10); % 这里复用第三章的蒙特卡洛主循环, 统计 cov_count/N_total Pc_list(t) runMonteCarlo(theta, L, lambda_b, N_real); end semilogy(theta_db_list, Pc_list, o-)这段代码里 runMonteCarlo 就是把第三章主循环封装成函数输入阈值、区域边长、基站密度和实现次数返回覆盖概率。注意横轴是 dB 值纵轴用 log 坐标因为覆盖概率在阈值升高时衰减很快线性坐标下高频段曲线会被压扁。和解析解对比时把 3.1 节的闭式公式也画在同一张图上两条线贴得越近说明你的仿真实现越可信。4.2 ASE 的统计口径和 lambda_u 扫描ASE 是面积频谱效率物理含义是每平方公里每秒每赫兹能传多少比特。标准定义是用户密度乘上每个用户的平均可达速率。蒙特卡洛里常见的近似算法是对每个成功传输的用户记 log2(1θ)失败记 0求平均再乘上用户密度。这个近似的前提是“成功就按阈值速率传失败就不传”。实际系统会用自适应调制速率随 SIR 连续变化但论文里为跟解析解对齐普遍用这个简化口径。ASE_theta_5dB_lambda_u.m 这个脚本的逻辑就是固定 θ5 dB让 lambda_u 从 1 扫到 100 个/km²画出 ASE 曲线。lambda_u_list [1 2 5 10 20 50 100]; % 用户密度 个/km² theta_db 5; theta 10^(theta_db/10); for k 1:length(lambda_u_list) ASE_k computeASE(lambda_b, lambda_u_list(k), theta, N_real); ase_list(k) ASE_k; end semilogx(lambda_u_list, ase_list, o-)computeASE 的内部实现要点用户数用泊松抽样即 n_ue poissrnd(lambda_u/1e6 * L²)每个用户算 SIR成功则累加 log2(1theta)最后除以区域面积平方公里。注意用户密度和基站密度的单位都是个/km²区域面积换算成 km² 时是 (L/1000)²这个换算漏掉会让 ASE 差 6 个数量级。ASE 曲线的典型形态是倒 U 型用户密度太低时频谱资源闲置ASE 上不去密度太高时干扰主导成功概率骤降。峰值对应的密度就是网络的最优用户承载量。论文里这张图通常还有一条“无依赖模型”的对比曲线差异出现在高密度区——空间依赖让用户在基站附近聚集干扰结构变化ASE 曲线峰值位置会偏移。4.3 出图环节fig14.m 这类脚本的读法和改造建议fig14.m 这类脚本放在压缩包里是用来复现论文图的但直接 run 很可能报“变量未定义”。原因是这些 fig 脚本从主仿真的工作区取变量而不是自己重新计算。正确读法分三步先看脚本里引用了哪些变量名再回主仿真脚本找对应的计算位置最后把主仿真跑完后再执行 fig 脚本。我倾向于把 fig 脚本改造成函数把引用的变量变成入参。这样不用每次跑完整主流程才能画图调试效率高很多。改造时注意脚本里如果嵌了循环把 N_real 调成 200 先验证形状形状对了再拉满否则一次图要等半小时改个参数又得重来。提示fig14_zoom.jpg 和 fig14_zoom.fig 是局部放大图通常对应曲线在某个参数区间的细节。如果你的仿真曲线整体形态对、但局部有毛刺先把 N_real 加一倍再对比毛刺大概率是蒙特卡洛方差造成的不是代码逻辑问题。5. 避坑记录五个让仿真翻车的细节每条都是血泪经验这套包我前前后后跑过三遍前两遍跑出来的图和作者的 jpg 对不上最后定位到的全是细节问题。写在这章给后来者省掉这些弯路。5.1 直接 run fig 脚本先跑主仿真否则变量未定义现象双击 fig14.m 运行MATLAB 直接报错提示未定义函数或变量比如找不到 X、Pc 这类变量。原因fig 系列脚本假设工作区里已经有了主仿真算出来的变量它们只是“画图器”不是“计算器”。直接运行工作区是空的自然报错。解决先跑 PPP_dep_sim.m 或 test.m 生成基础变量再执行 fig 脚本。更稳妥的办法是把 fig 脚本改造成函数把需要的变量作为入参传入彻底摆脱对全局工作区的依赖。改造后的调用方式类似 fig14_reproduce(Pc, theta_db_list)一眼就能看清依赖关系。5.2 曲线毛刺多蒙特卡洛的玄学在次数不在算法现象覆盖概率曲线不光滑相邻阈值点之间跳来跳去像是被狗啃过。原因实现次数 N_real 太小或者每个 realization 里撒的用户数太少。100 次实现以内的蒙特卡洛方差能大到掩盖真实趋势500 次以上才勉强平滑。解决N_real 提到 2000~5000每个 realization 的用户数提到 100 以上。如果想减少计算量可以先用 500 次把逻辑跑通确认没问题再加量。曲线平滑度的估算参考覆盖概率的标准差大约正比于 sqrt(Pc(1-Pc)/N)Pc0.5 时要 2500 个样本才能做到 1% 的标准差。5.3 边界效应边缘用户干扰被低估覆盖概率虚高现象仿真出的覆盖概率比解析解高出一截尤其是阈值较低时偏差更明显。原因L×L 方形区域边缘的用户可视范围内的基站比中心用户少干扰被系统性低估。边缘越窄的区域这个偏差越严重。解决加保护带。撒基站时用 (L2g)² 区域统计覆盖时只看中心 L² 区域的用户。g 取 300 米左右对应典型宏站覆盖半径的量级。统计时用户坐标要限制在中心区域内否则边缘用户同样会污染统计结果。这个“撒大算小”的做法是随机几何仿真的标准配置。5.4 单位换算三连坑密度、dB、功率量纲现象覆盖概率要么接近 1要么接近 0曲线完全不是论文里那种中间过渡形态。原因三个单位问题叠加。基站密度 20 个/km² 直接当成 20 个/m² 用结果是基站密度爆炸、全覆盖5 dB 阈值没转线性直接用 5 参与比较信号功率用瓦、噪声功率用毫瓦分子分母量纲差 1000 倍。解决统一换算法则——密度一律先除以 1e6 转成个/m²dB 阈值用 10^(x/10) 转线性所有功率量纲统一成 mW。我在代码开头固定写一段换算注释把三个量纲列出来每次改参数前先过一遍这种低级翻车基本能根除。5.5 路径损耗指数 α 取太小干扰期望直接发散现象仿真跑出来的 SIR 普遍偏低加大基站密度也救不回来曲线形状完全走样。原因路径损耗指数 α 要大于 2干扰功率的期望才收敛。α 取 1.5 或 2 时远处基站的干扰贡献衰减太慢累加和发散SIR 被无限拉低。解决城区场景 α 取 3.5~4郊区取 2.5~3。这个参数物理含义很明确α 越小信号随距离衰减越慢基站间的干扰越强。蜂窝网络分析里 α2 的场景很少见如果论文或课程设计要求你扫 α把下限锁在 2.5 以上否则结果没有参考价值。6. 进阶把单点参数扩成参数扫描并用固定随机种子验证结果第三章和第四章的代码都是单点仿真一组参数出一条曲线。论文里真正需要的往往是参数扫描图——横轴是密度或阈值纵轴是性能指标。把单点仿真封装成函数后扫描就是一个循环加一个数组的事。这里给出一个完整的密度扫描套路代码可以直接嵌进 ASE_theta_5dB_lambda_u.m 这类脚本里。% 固定随机种子, 保证不同密度下对比公平 rng(42); lambda_u_list [1:1:10 20:10:100]; % 个/km², 稀疏到稠密 for k 1:length(lambda_u_list) [Pc_k, ase_k] computeASE(lambda_b, lambda_u_list(k), theta, N_real); pc_list(k) Pc_k; ase_list(k) ase_k; end % 保存中间结果, 避免重复跑仿真 save(ase_sweep_result.mat, lambda_u_list, pc_list, ase_list, theta_db);这里最关键的技巧是 rng(42)。固定随机种子后不同密度下的仿真共享同一套随机数序列曲线对比时才不会因为“这一组仿真恰好运气好”产生虚假差异。这是蒙特卡洛对比实验里最容易忽略的细节——不固定种子两条曲线之间的差距可能完全是噪声根本反映不了参数影响。computeASE 内部每个密度点都要跑满 N_real 次蒙特卡洛全部跑完可能要十几分钟。两个加速手段一是把 for 循环改成 parfor并行池开起来后速度提升接近核数倍数二是把每次 realization 的结果实时追加到一个数组中途中断也能看到已完成的阶段性数据不用全部推倒重来。从那以后我每次改密度或阈值参数都会强制先固定 rng再跑一组对比顺手把中间结果存成 .mat。这个习惯帮我在论文返修时省了大量重跑时间——审稿人要你补一条曲线你只需要改参数重新跑一次而不是对着旧结果从头猜。希望帮到你。本文还有配套的精品资源点击获取