简介这是一份基于MATLAB 2021a的可见光通信仿真资源面向通信工程、光学无线传输领域的研究者与高年级本科生适合用于室内可见光通信链路建模与性能分析。资源以室内环境为场景通过配置LED布局、接收机视角、朗伯辐射阶数、墙面/天花板/地板反射率等参数输出误码率空间分布与信噪比空间分布帮助读者直观理解不同空间位置的光照覆盖与通信质量差异。压缩包共9个文件包含6个.m脚本、2张结果图片和1段avi操作录像整体仅3.06MB脚本涵盖主运行程序、可见光链路视距与非视距损耗计算、朗伯辐射函数等模块便于按流程复现、断点调试与二次修改。录像详细演示了MATLAB当前文件夹路径设置与运行全过程能够有效规避新手常见环境配置问题。这套资源已有926人学习适合想快速搭建可见光通信仿真实验、深入分析空间分布特性并验证算法效果的读者。1. 可见光通信仿真先做空间分布别只在单点上算误码率如果只是把接收机固定在天花板正下方可见光通信的信噪比和误码率用笔算公式几分钟就能出结果。但实际场景里 LED 灯具是离散安装在房间顶部的接收平面在桌面上信号强度、噪声基底和误码率在空间上会呈现明显的“亮区”和“暗区”。标题里的“误码率空间分布、信噪比空间分布”要求的是网格化接收区域上的性能热力图而不是单点性能值。这种仿真要解决三个具体问题信道增益怎么空间化、噪声怎么逐点计算、结果怎么画成可对比的分布图。本文按系统模型、信噪比分布、误码率分布、脚本组织和录像交付这条线展开代码可以直接在 MATLAB 里跑通适合做可见光通信课程设计、科研预研和灯具布局方案对比的工程师参考。2. 可见光通信信道模型与接收功率分布朗伯辐射和坐标网格是关键2.1 房间坐标与 LED 布局先定坐标系再写接收平面循环可见光通信仿真的第一步不是写公式而是把房间、LED 和接收平面统一到一个坐标系里。常见做法是以房间一角为原点x、y 轴沿地面z 轴向上。LED 是发射端安装在天花板接收机是一个光电检测器PD放在距地面 0.85 米左右的桌面上。接收平面不是单点而是用 meshgrid 生成的一个网格网格点密度直接决定分布图的精细程度和计算耗时。我一般先把以下参数集中定义方便后续修改布局做对比。参数默认值说明房间尺寸5 m × 5 m × 3 m典型办公室尺寸LED 数量4 个对称吸顶安装LED 坐标(1,1,3)、(4,1,3)、(1,4,3)、(4,4,3)单位米接收平面高度0.85 m桌面高度网格范围与步长x/y 从 0.5 到 4.5步长 0.1200×200 量级半功率角半角60°朗伯辐射阶数相关PD 视场角 FOV60°超过该角度不接收接收平面高度取 0.85 而不是地面是因为真实使用场景里终端放在桌面上。网格范围如果直接取 0 到 5 米墙角的点在 FOV 判定上经常出现极端值画出来的分布图边缘一圈没有参考意义所以常见做法是留出 0.5 米边距。2.2 朗伯辐射阶数与 LoS 信道增益单位、角度、FOV 三个坑可见光通信中最常用的是朗伯辐射模型。LED 的辐射方向性由朗伯辐射阶数 m 描述公式是[ m -\frac{\ln 2}{\ln(\cos \Phi_{1/2})} ]其中 (\Phi_{1/2}) 是半功率角半角。这个公式里的 (\ln) 是自然对数很多人在 MATLAB 里写作 log10 导致 m 变成负数这是第一个坑。光线到达接收面后直射链路LoS信道增益写为[ H(0) \frac{(m1)A}{2\pi d^2} \cdot \cos^m(\phi) \cdot \cos(\psi) \cdot \text{rect}\left(\frac{\psi}{\psi_{\text{FOV}}}\right) ]这里 d 是 LED 到接收点的距离φ 是 LED 出射角ψ 是接收点入射角。rect 函数表示只有当入射角小于 PD 视场角 FOV 时才有增益否则直接置零。第二个坑是角度计算用点积求余弦时要先用向量归一化否则 acos 传入超出 [-1,1] 的值会得到 NaN。第三个坑是单位cos 函数传弧度但表达式里不涉及显式角度转换容易忽略。2.3 接收功率分布的矩阵化实现接下来把上述公式写成矩阵化代码。先准备一个函数计算所有网格点到某个 LED 的增益和接收功率。注意这里用隐式扩展MATLAB R2016b 及以上不需要写双重 for 循环。% compute_vlc_channel.m % 输入led_pos(1x3), rx_grid_x(NxM), rx_grid_y(NxM), params(结构体) % 输出H(信道增益矩阵), Pr(接收功率矩阵) function [H, Pr] compute_vlc_channel(led_pos, rx_x, rx_y, params) % 提取参数避免函数体里到处写 params.xxx A params.pd_area; % PD 有效面积m^2 FOV params.fov_deg * pi/180; % 视场角转弧度 Phi params.half_power_deg * pi/180; % 半功率角半角转弧度 Pt params.led_power; % LED 光功率W % 朗伯辐射阶数注意是自然对数 m -log(2) / log(cos(Phi)); h_rx params.rx_height; % 接收平面高度 % 接收点坐标网格 [M, N] size(rx_x); % 把网格展成列向量方便统一计算 rx_vec [rx_x(:), rx_y(:), ones(M*N,1)*h_rx]; % 从 LED 指向接收点的向量 dx rx_vec(:,1) - led_pos(1); dy rx_vec(:,2) - led_pos(2); dz rx_vec(:,3) - led_pos(3); d_vec [dx, dy, dz]; d sqrt(dx.^2 dy.^2 dz.^2); % 距离m % LED 法线方向为 (0,0,-1)向下 % 出射角余弦 向量与 LED 法线的点积 / 距离 cos_phi (-dz) ./ d; % 等价于 cos(出射角) % 接收面法线方向为 (0,0,1)向上 % 入射角余弦 向量与接收面法线的点积 / 距离 cos_psi dz ./ d; % 超出 FOV 的接收点信道增益置 0 in_fov cos_psi cos(FOV); % 计算信道增益 H_vec (m1) * A ./ (2*pi*d.^2) .* cos_phi.^m .* cos_psi; H_vec(~in_fov) 0; % 接收功率 W Pr_vec Pt .* H_vec; % 恢复成网格形状 H reshape(H_vec, M, N); Pr reshape(Pr_vec, M, N); end这段代码里最关键的是 cos_phi 和 cos_psi 的符号。LED 向下发射出射方向与 LED 法线方向z 轴负方向的夹角余弦用 -dz/d 计算没有问题。接收面朝上入射方向与接收面法线z 轴正方向的夹角余弦是 dz/d。有些教材推导时把两个余弦都写成 |dz|/d这在水平接收面时结果相等但如果以后扩展倾斜接收面就会出错建议从一开始就按向量法线来写。多个 LED 叠加时接收功率是各 LED 贡献的直接相加。写一个脚本调用上面的函数即可% 初始化房间与网格 params.led_power 10; % 单颗 LED 光功率 10W params.pd_area 1e-4; % PD 面积 1 cm^2 params.fov_deg 60; params.half_power_deg 60; params.rx_height 0.85; x 0.1:0.1:4.9; y 0.1:0.1:4.9; [rx_x, rx_y] meshgrid(x, y); [rows, cols] size(rx_x); led_pos [1 1 3; 4 1 3; 1 4 3; 4 4 3]; Pr_total zeros(rows, cols); for led_idx 1:size(led_pos, 1) [~, Pr] compute_vlc_channel(led_pos(led_idx,:), rx_x, rx_y, params); Pr_total Pr_total Pr; end % 绘制接收功率分布 figure; pcolor(rx_x, rx_y, Pr_total*1e3); shading interp; colorbar; xlabel(x (m)); ylabel(y (m)); title(接收功率分布 (mW));功率分布图上能明显看到 LED 正下方的热点和房间角落的暗区。如果暗区面积过大说明 LED 布局不合理。这里网格步长 0.1 米200×200 量级的网格在普通电脑上运行时间可以接受比逐点双重循环快两个数量级以上。3. 信噪比空间分布散粒噪声和热噪声分开算SNRdB 图才可信3.1 光电转换与噪声模型影响信噪比分布的两个噪声源接收功率分布算出来后下一步是把它转换成电信号和噪声。PD 把光功率转换成光电流即 (I_{\text{ph}} R \cdot P_r)其中 R 是响应率单位 A/W。接收端噪声主要有两个来源散粒噪声和热噪声。散粒噪声来自光电流的随机性包括信号光和背景光两部分双边功率谱密度下的方差为[ \sigma_{\text{shot}}^2 2q \cdot (I_{\text{ph}} I_{\text{bg}}) \cdot B ]其中 q 是电子电荷 (1.602\times10^{-19}) C(I_{\text{bg}}) 是背景光电流B 是接收机带宽。热噪声来自前端放大电路常见近似表达式取第一项[ \sigma_{\text{thermal}}^2 \approx \frac{8\pi k T}{G_{\text{ol}}} \cdot \eta \cdot A \cdot I_2 \cdot B^2 ]其中 k 是玻尔兹曼常数T 是温度(G_{\text{ol}}) 是开环电压增益η 是 PD 电容相关因子A 是 PD 面积(I_2) 是噪声带宽因子。第二项与 B³ 成正比在高速系统里才明显低速仿真可以忽略。噪声参数默认值说明背景光电流 (I_{\text{bg}})200 μA室内荧光灯/日光背景响应率 R0.53 A/W典型 PIN PD 在 450nm带宽 B10 MHz低速 OOK 系统温度 T300 K室温开环增益 (G_{\text{ol}})10典型跨阻放大器噪声带宽因子 (I_2)0.562矩形脉冲对应的值注意 (I_{\text{bg}}) 的取值。背景光电流是常数不随信号位置变化。在 LED 正下方信号强背景光占比小在角落信号弱散粒噪声被背景光主导这就导致 SNR 分布比功率分布衰减得更快。3.2 电 SNR 定义与 dB 分布图线性功率比转 dB 后才有对比意义电 SNR 的定义是信号光电流的均方值与总噪声方差之比[ \text{SNR}e \frac{(R \cdot P_r)^2}{\sigma{\text{shot}}^2 \sigma_{\text{thermal}}^2} ]这里用的是光电流的平方对应电功率。有些资料里用平均光功率计算得到的是平均 SNR有些用 OOK 峰值功率算差别正好是 3 dB。后面的误码率公式里我把这个 SNR 当作平均电 SNR 使用写代码时保持这一个约定整条链路才不会乱。% 计算信噪比空间分布 q 1.602e-19; k_B 1.381e-23; R 0.53; % A/W B 10e6; % Hz I_bg 200e-6; % 背景光电流 A % 散粒噪声 I_ph R * Pr_total; sigma_shot2 2 * q * (I_ph I_bg) * B; % 热噪声第一项近似 T 300; G_ol 10; eta 0.9; I2 0.562; sigma_therm2 (8*pi*k_B*T/G_ol) * eta * params.pd_area * I2 * B^2; % 总噪声 sigma_total2 sigma_shot2 sigma_therm2; % 电 SNR线性值与 dB 值 SNR_lin (R * Pr_total).^2 ./ sigma_total2; SNR_dB 10 * log10(SNR_lin); % 空间分布热力图 figure; pcolor(rx_x, rx_y, SNR_dB); shading interp; colorbar; clim([10 50]); % 固定刻度范围便于对比不同布局 xlabel(x (m)); ylabel(y (m)); title(信噪比空间分布 (dB));绘图时有两件事要注意。第一是使用 pcolor 加 shading interp 而不是 surf因为 pcolor 是俯视热力图surf 的高度信息在这里没有意义。第二是固定 colorbar 的范围即 clim 函数。如果不固定每次重新计算布局后颜色映射会自动拉伸视觉上会误导对比。3.3 布局对信噪比空间分布的直接影响单灯与四灯的分布差异用同一个脚本把 led_pos 改成单个 LED 放在房间中心 (2.5, 2.5, 3)其余参数不变会看到 SNR 分布从四灯时的“四个峰值”变成“中心一个峰、边缘梯田状下降”。单灯方案中心点 SNR 比四灯方案高因为所有功率集中在一处但边缘角落 SNR 会掉到 20 dB 以下。四灯方案虽然峰值略低但整体覆盖更均匀。这就是做空间分布仿真最重要的产出不是某一处的 SNR 数值而是覆盖均匀性的直观判据。实际工程里调整布局时我会把每个候选布局的 SNR_dB 矩阵存下来统计低于某个门限的像素点占比再结合分布图决定取舍。这些统计在下一章算误码率时可以直接复用。4. 误码率空间分布OOK 的 BER 公式、Q 函数与等值线边界4.1 OOK 误码率公式与 Q 函数误码率和误信率在 OOK 下是同一个数可见光通信最简单的调制是 OOK即开关键控。光强代表二进制 1无光代表二进制 0。在 AWGN 信道下OOK 的误码率理论公式为[ P_e Q\left(\sqrt{\frac{\text{SNR}_e}{2}}\right) ]其中 (Q(x)\frac{1}{2}\operatorname{erfc}(x/\sqrt{2}))SNR_e 是上一章算出的平均电 SNR。这里分母里的 2 是 OOK 信号平均功率和峰值功率的折算因子照着这个式子用就不会出问题。如果误把公式写成 (Q(\sqrt{\text{SNR}_e}))相当于高估了 1.5 dB 左右的性能画出来的分布图整体偏乐观。还要说清一个概念OOK 每个符号只承载 1 个比特所以误码率等于误信率误比特率。热词里常搜的“误码率和误信率的大小关系图”在这里的最直观结论就是对 OOK 二者相等对 PAM4、QAM 这类多进制调制误码率大于误信率需要按格雷映射换算。本文默认只做 OOK遇到高阶调制时另写一套映射函数。4.2 误码率网格计算与范围裁剪高 SNR 区域的 BER 值先做饱和处理把 SNR_lin 带入 BER 公式得到误码率空间分布。问题在于 SNR 高的时候 BER 会小到 1e-20 以下直接画图时颜色条会被极小值拉扁整个图看起来像只有两种颜色。常见的处理方式是把 BER 矩阵做上下限裁剪下限取 1e-8上限取 0.5无信息时的误码率然后画 log10 坐标。% 计算误码率空间分布 SNR_safe max(SNR_lin, 1e-10); % 防止 sqrt 里出现极小值 BER qfunc(sqrt(SNR_safe / 2)); % OOK 误码率 % 裁剪到 [1e-8, 0.5]避免色标被拉扁 BER_clipped min(max(BER, 1e-8), 0.5); % 用对数色标画分布 figure; surf(rx_x, rx_y, log10(BER_clipped), EdgeColor, none); view(2); % 俯视图 colormap(parula); colorbar; clim([-8 -1]); xlabel(x (m)); ylabel(y (m)); title(误码率空间分布 (log10 BER));使用 qfunc 是 MATLAB 通信工具箱的函数如果机器里没有安装通信工具箱可以自己定义BER 0.5 * erfc(sqrt(SNR_safe / 2) / sqrt(2));这个替代写法在纯 MATLAB 基础环境里就能运行。色标单位是 log10(BER)clim 设为 [-8 -1]对应 BER 从 1e-8 到 0.1。0.5 被裁剪到 -0.301 但不参与色标范围所以暗区显示为最深色不会破坏整体对比。4.3 用等值线提取覆盖边界1e-6 与 1e-3 两条线画在空间图上误码率分布的色块图能看出大致趋势但工程上更需要“哪些区域可用、哪些区域不可用”的边界。做法是在 pcolor 图上叠加两条等值线一条是 1e-6对应无 FEC 的可靠通信门限一条是 1e-3对应带前向纠错后仍可恢复的门限。figure; pcolor(rx_x, rx_y, log10(BER_clipped)); shading interp; colormap(parula); colorbar; clim([-8 -1]); hold on; % 提取等值线并标注 [~, h1] contour(rx_x, rx_y, log10(BER_clipped), [log10(1e-6) log10(1e-6)], ... LineColor, w, LineWidth, 1.5, ShowText, on); [~, h2] contour(rx_x, rx_y, log10(BER_clipped), [log10(1e-3) log10(1e-3)], ... LineColor, k, LineWidth, 1.5, LineStyle, --, ShowText, on); legend([h1, h2], BER 1e-6, BER 1e-3); xlabel(x (m)); ylabel(y (m)); title(误码率空间分布与覆盖边界);注意等值线取的是 log10(BER) 矩阵所以等值线数值要写成 log10(1e-6)直接在原矩阵上画阈值 1e-6 是不行的。白色实线标出可靠通信区域黑色虚线标出 FEC 门限。两张图叠在一起房间的“通信可用区”一目了然。这个覆盖边界图是我每次调整 LED 布局后必看的第一张图比单纯看 BER 数值直观得多。5. 参数化仿真脚本与矩阵化提速避免仿真发散和“算得慢”5.1 仿真参数结构体一次定义一个默认配置布局对比时只改坐标前面章节的代码把参数散落在脚本各处改一个参数要翻好几处。做布局对比时常见做法是写一个初始化函数 set_params把全部参数集中放在结构体里每次跑仿真只改需要对比的字段。function params set_params() % 默认仿真参数集中管理 params.room [5 5 3]; % 房间尺寸 m params.led_pos [1 1 3; 4 1 3; 1 4 3; 4 4 3]; % 每个 LED 坐标 params.led_power 10; % 单颗 LED 光功率 W params.pd_area 1e-4; % PD 面积 m^2 params.rx_height 0.85; % 接收平面高度 m params.fov_deg 60; % PD 视场角全角 params.half_power_deg 60; % LED 半功率角 params.background_current 200e-6; % 背景光电流 A params.bandwidth 10e6; % 接收机带宽 Hz params.responsivity 0.53; % PD 响应率 A/W params.grid_step 0.1; % 空间网格步长 m end对比单灯和四灯布局时主脚本只需要写params set_params(); params.led_pos [2.5 2.5 3]; % 单灯居中 [SNR_dB, BER_log] run_vlc_sim(params); % 主仿真函数 % 四灯布局 params set_params(); % 恢复默认 [SNR_dB, BER_log] run_vlc_sim(params);主仿真函数 run_vlc_sim 把前面信道计算、噪声计算、误码率计算串起来返回分布矩阵。这个结构的价值是参数变更时不容易漏改某处也方便后面做录像时在循环里反复调用。5.2 循环仿真 vs 矩阵广播网格变大后的性能差距是数量级的初学 MATLAB 时容易把信道计算写成三重循环遍历所有网格点和所有 LED每层循环里算一次增益。这样写 200×200 网格加 4 个 LED 就是 16 万次标量计算加上角度求反三角和条件判断单次运行要数十秒。矩阵化写法把所有接收点当成列向量一次完成全部计算通常能快一到两个数量级。% 对比两种写法的耗时 % 循环写法示意 tic; H_loop zeros(length(x), length(y)); for i 1:length(x) for j 1:length(y) for k 1:4 % 标量计算信道增益并累加 end end end t_loop toc; % 矩阵化写法 tic; % 上面 compute_vlc_channel 中的向量化代码 t_vec toc; fprintf(循环耗时 %.2f s矩阵化耗时 %.2f s\n, t_loop, t_vec);这个对比在生产脚本里非常值得做。网格从 200×200 变成 500×500 时矩阵化版本内存占用会升高但计算时间仍在可接受范围循环版本则基本不可用。内存方面500×500 网格展成 25 万×3 的矩阵只有几 MB 量级完全不用担心。5.3 仿真发散与异常结果排查单位、FOV、NaN 这三处先看仿真“发散”在静态信道模型里一般不是数值爆炸而是结果明显不合理。最常见的三类异常和排查顺序如下。第一是接收功率出现 NaN。几乎都是角度计算时 acos 越界或者 led_pos 的前两维坐标与接收平面坐标重叠导致某点距离为零。排查办法是在 compute_vlc_channel 里检查 d 是否含零if any(d(:) 0) error(存在接收点与 LED 重合请检查坐标); end第二是 SNR 全部偏低或偏高 60 dB 这类倍数级偏差。先查功率单位LED 光功率如果用 mW 而 PD 面积用 m²接收功率会差 30 dBSNR 差 60 dB。常见做法是全部使用国际单位 W、m、m²画图时再转换单位。第三是 BER 恒为 0.5相当于所有点都判断为随机猜测。这种情况几乎都是 FOV 判定写反了条件比如把 in_fov 判断成入射角小于某个阈值时置零导致所有点都被清零。在脚本里插入如下自检可以快速定位% 自检房间中心距离 LED 最近的点SNR 应明显大于其他位置 if SNR_dB(round(size(SNR_dB,1)/2), round(size(SNR_dB,2)/2)) 25 warning(房间中心 SNR 异常偏低检查 FOV 或噪声参数); end仿真环境里没有真实硬件的“实际值”自检的意义在于锁定相对关系的合理性。中心点 SNR 高于角落、LED 正下方出现局部峰值、BER 与理论单点计算一致这三个特征验证通过后整套脚本才进入录像交付阶段。6. 仿真操作录像与结果复现把空间分布仿真包交付成可验证的产物6.1 用 VideoWriter 录制网格收敛动画标题里“含仿真操作录像”的交付物最稳定可靠的做法是用 VideoWriter 把仿真过程录成视频而不是单纯用录屏软件录操作。录屏软件的缺点是容易录进无关桌面内容画面压缩格式不统一。用 VideoWriter 可以在仿真脚本内部逐帧生成图像保证录制的就是仿真结果本身。% 生成网格收敛动画并录制成 MP4 v VideoWriter(ber_distribution.avi); v.FrameRate 5; % 每帧代表一个网格粒度 open(v); grid_steps [0.5 0.2 0.1 0.05]; % 逐步细化网格 for step grid_steps params.grid_step step; [~, BER_log] run_vlc_sim(params); figure(Visible, off); % 设置不可见窗口加快录制速度 pcolor(rx_x, rx_y, BER_log); shading interp; colormap(parula); clim([-8 -1]); title(sprintf(网格步长 %.2f m, step)); colorbar; frame getframe(gcf); writeVideo(v, frame); close(gcf); end close(v);录制时固定 clim 范围是这套流程里最重要的细节。如果不固定每个网格粒度计算出的 BER 最大值与最小值不同MATLAB 会自动重映射颜色视频播放时会产生亮暗闪烁观感极差也容易让看视频的人误判仿真相位变化。固定范围后视频呈现的是同一色标下网格从粗糙到精细的收敛过程信息密度高很多。6.2 脚本执行录像的常见做法与输出物的组织除了结果动画完整交付还需要一份展示“如何修改参数、如何运行脚本、如何在命令行查看结果”的操作录像。MATLAB 本身不自带屏幕录制功能常见做法是用系统录屏工具录制命令行窗口和编辑器窗口录制前先确认脚本可以一次跑通不要边录边排错。一个可复现的仿真包建议包含set_params.m、run_vlc_sim.m、compute_vlc_channel.m 三个脚本一份 BER 与 SNR 分布图输出以及一段网格收敛动画。把这三个脚本理清楚后即便换一个人用不同的 MATLAB 版本打开也只改参数函数里几个数值就能复现所有分布图。对常见的那类“把整个任务交给 codex 之类的工具直接生图”的做法这里的经验是脚本能生成单张图容易但要想让生成过程可回放、可对比边界还是要靠 VideoWriter 和控制好色标范围这套组织方式。录像文件体积上步长 0.1 米的网格画面对应的 AVI 文件一般几百 MB条件允许时优先输出 MP4 格式并设置较低帧率需要在参数对比时再单独截图保留细节。本文还有配套的精品资源点击获取