基于MATLAB的CRI显色指数计算：从SPD光谱到Ra的完整流程

简介针对照明设计与光学研究中的光谱功率分布SPD与显色性指数CRI计算需求这套MATLAB程序为照明工程师、LED研发人员及光学专业学生提供了轻量工具。代码通过解析光谱测量数据自动完成波长筛选、噪声过滤和功率分布计算并输出标准显色指数Ra与补充指标R9便于快速评估光源的色彩还原能力。资源包共2个文件包含一个MATLAB脚本文件.m和一个频谱数据文件.asv压缩包仅1KB结构精简、易于修改。目前已有767人学习下载适用于博物馆照明、摄影、植物生长灯及医疗照明等需要精准显色控制的场景。使用者可直接运行脚本处理自定义光谱数据结合输出结果比较不同光源的SPD与CRI从而辅助LED芯片选配或灯具光色优化兼顾科研分析与工程应用是照明设计中的实用小工具。1. 用 MATLAB 算显色性 CRI首先要把 SPD 光谱当成一手数据而不是成品做过光源检测的人都会有同感积分球或光谱仪导出的 SPD 文件打开之后是一排波长、一排功率值看着就一长串数字。但显色性 CRI 的整个计算逻辑恰恰是从这串 SPD 出发的先得到光源的三刺激值 XYZ再按相关色温匹配参考光源最后拿 14 个标准色样在两个光源下比色差折算出 Ri 和 Ra。仪器自带的软件能直接给结果但如果你想自己评估一批光源、批量跑数据或者准备一份能放进汇报 PPT 的光谱对比图绕不开自己写 MATLAB 程序。下面这套流程是从实测 SPD 到 Ra 和 R1–R14 的常见落地写法也把最容易算错的地方一起说清楚。2. 显色性 CRI 计算的第一步SPD 读入、插值、XYZ 积分2.1 波长轴对不齐后面每一步都会带偏差CRI 计算的理论基础是把测试光源的 SPD 和标准观察者色匹配函数逐点相乘再积分。问题在于仪器导出的 SPD 步长五花八门有的每 5nm 一个点有的每 2nm有的起始波长 360nm有的只测到 760nm。色匹配函数表也常见 1nm 和 5nm 两种版本。只要波长轴不一样trapz积分出来的 X、Y、Z 就没法稳定比较。我在写显色性 CRI 计算程序时第一件事永远是先把所有数据插值到同一套 1nm 网格上。网格范围取 380–780nm这是 CRI 标准要覆盖的可见光区间。插值方法用pchip它能保留窄带荧光灯和 LED 的峰形线性插值会让峰值被拉平样条插值又容易在峰附近振出负值。% 读取仪器导出的 SPD 文件第一列波长第二列功率 raw readtable(spd_measured.csv); wl0 raw{:, 1}; s0 raw{:, 2}; % 标准计算网格380~780 nm步长 1 nm wl (380:780); % 丢掉范围之外的点避免外插 valid (wl0 380) (wl0 780); wl0 wl0(valid); s0 s0(valid); % PCHIP 插值到 1nm范围外补 0 spd interp1(wl0, s0, wl, pchip, 0);这段代码里raw用的是readtable对 Excel、CSV、TXT 都能读关键是列顺序。如果你的文件是波长 nm, 功率 W/sr/nm这种带单位的表头raw{:, 1}和raw{:, 2}仍然能取到数值矩阵。插值之后spd和后面加载的色匹配函数共享同一个wl就不需要再做索引对表。2.2 相对功率也能算 CRI但要先把 XYZ 归一化成 Y100SPD 的绝对辐亮度值对 CRI 没有意义因为显色性评价只涉及光谱形态不涉及光通量大小。一般做法是先把 X、Y、Z 按 Y 值归一到 100再计算色坐标。颜色匹配函数我这里不内置到代码里通常工程上会在项目目录放一份ciexyz_1nm.csv列名分别是wl、xbar、ybar、zbar。读取和积分是这样的cie readtable(ciexyz_1nm.csv); xbar cie.xbar; ybar cie.ybar; zbar cie.zbar; % 梯形法积分得到三刺激值 X trapz(wl, spd .* xbar); Y trapz(wl, spd .* ybar); Z trapz(wl, spd .* zbar); % 相对 SPD 归一化到 Y100 scale 100 / Y; X X * scale; Y 100; Z Z * scale; % 色度坐标 x X / (X Y Z); y Y / (X Y Z);这里用trapz是最常见的做法它按梯形求积对 1nm 网格已经足够。spd .* xbar中spd必须和xbar同一维度所以前面把插值网格严格固定成 380:780 就很重要。归一化那行是关键scale 100 / Y把整个 SPD 的“高度”统一到同一标准上。两组不同功率的 SPD只要形状相同归一化后 X、Y、Z 的比值相同最后算出的 CRI 也一致。量公式在 CRI 计算里的角色X∫ S(λ) x̄(λ) dλ色度计算基础Y∫ S(λ) ȳ(λ) dλ归一化基准Y100Z∫ S(λ) z̄(λ) dλ色度计算基础x, yX/(XYZ), Y/(XYZ)求相关色温 CCT 和参考光源白点做好这一步之后你已经有了一个“洗净”的光源光谱特征。下一章开始才是 CRI 真正区别于普通色度计算的地方怎么找参考光源以及怎么处理人眼的适应性。3. CRI 计算程序的第二步参考光源与色适应3.1 为什么不能直接把标准色卡的 XYZ 拿来比很多人第一次写显色性 CRI 计算程序时拿着标准色卡在 D65 下的 XYZ 表再用测试光源算一遍 XYZ直接求色差。这样得出的数值会非常离谱尤其是暖色温光源色差会大到你不敢相信。问题出在人眼的色适应能力。同一块红布在 3000K 白炽灯下和在 6500K 日光下人眼看过去不会觉得是两种完全不同的颜色因为大脑会把白点“拉回来”。CRI 设计的初衷就是剔除这种色适应带来的偏移只看光源对物体颜色还原的额外偏差。所以必须给每个测试光源配一个“参考光源”这个参考光源的色温要跟测试光源接近标准满色度指数算出来的结果才有可比性。按照 CIE 13.3 的常见做法CCT 在 5000K 以下用黑体辐射 SPD 做参考5000K 以上用 CIE 日光 SPD。实际工程里很多 LED 色温落在 2700K、3000K、4000K所以黑体参考用得相当多。日光区段我一般直接预置 D50、D55、D65、D75 的表格按最接近的色温去取。3.2 在 MATLAB 里生成黑体参考 SPD黑体 SPD 用普朗克公式就能精确生成不需要查表。波长单位要先换算成米普朗克第一、第二辐射常数直接写进函数里。function spd blackbody_spd(T, wl) % T: 黑体色温单位 K % wl: 波长向量单位 nm c1 3.741771e-16; % W·m^2 c2 1.438777e-2; % m·K wl_m wl * 1e-9; spd c1 ./ (wl_m.^5) ./ (exp(c2 ./ (wl_m * T)) - 1); end这个函数出来后可以把spd换成参考光源的spd_ref再用第 2 章同样的积分方法得到Xr, Yr, Zr。需要提醒一点黑体 SPD 是功率密度的相对形状算 XYZ 时同样要归一化到 Y100。参考光源的 XYZ 后面要用它的白点做色适应变换所以这一步的准确性直接决定最终 Ra。3.3 色适应用 von Kries 近似CRI 结果才稳严格来说CIE 13.3 里面有一整套色适应计算流程但它本质上是在锥体响应空间里做缩放。工程计算里用 von Kries 的三通道缩放模型已经能覆盖绝大多数光源尤其对连续光谱的白炽灯和普通荧光粉 LED 效果很好。窄带 LED 的误差主要来自 SPD 峰形而不是色适应模型这一点后面再讲。色适应要做的事把测试光源下某个标准色样的 XYZ变换到“如果把参考光源看成白色”的条件下这个色样应该呈现的 XYZ。先求白点的 LMS 响应再把两个白点的比例作为缩放系数。% 常用 von Kries 变换矩阵近似值 M_vk [ 0.40024 0.70760 -0.08081; -0.22630 1.16532 0.04570; 0.00000 0.00000 0.91822 ]; % 测试光源白点和参考光源白点 test_white [X; Y; Z]; ref_white [Xr; Yr; Zr]; test_lms M_vk * test_white; ref_lms M_vk * ref_white; % 缩放系数 参考白 / 测试白 gain ref_lms ./ test_lms;真正应用时这个gain要对每个标准色样的 XYZ 都乘一遍。由于标准色样算出来的 XYZ 是三维列向量可以先按列乘再还原sample_adapted M_vk \ (gain .* (M_vk * sample_XYZ));我在实际项目里会把这段封装成一个函数输入是测试光源白点、参考光源白点和样本 XYZ输出是适应后的 XYZ。这样做的好处是主程序只需要关心循环和评分色适应矩阵换版本时也只改一个文件。这里要特别强调色适应不是可选项。我见过只算ΔE不加色适应的程序3000K 卤素灯 Ra 只能算到 60 多修正色适应后立刻回到 99 以上。如果你的 CRI 程序算出来白炽灯 Ra 不高十有八九是这一节出了问题。4. MATLAB CRI 主程序标准色样、Ri 循环与结果导出4.1 14 个标准色样怎么参与计算显色性 CRI 计算标准里用到 14 个色样其中前 8 个是中低饱和度的常规颜色对应一般显色指数 Ra后 6 个是红、黄、绿、蓝等高饱和度颜色提供 R9 到 R14 这些特殊显色指数。灯具厂的规格书里最常见的 R9指的就是第 9 个色样深红色。每个色样需要一条光谱反射率曲线。计算思路是用测试光源 SPD 乘以色样反射率再乘以颜色匹配函数积分得到该色样在测试光源下的 XYZ。参考光源那边同理。所以我建议把反射率矩阵单独放一个 CSV 文件行是波长列是 14 个色样。% 导入14个标准色样的反射率曲线 refl readtable(TCS_reflectance.csv); % 假设第一列是波长后面 14 列是 R1~R14 wl_ref refl{:, 1}; refMat refl{:, 2:15};这个refMat里的每一列都可以和之前得到的测试光源spd、xbar、ybar、zbar相乘后trapz积分。循环里同时计算测试光源和参考光源下的 XYZ然后套用上一章的色适应再计算色差。4.2 主循环和 Ri / Ra 的计算下面这段是显色性 CRI 计算程序的主循环框架。它把“积分求 XYZ → 色适应 → 算 ΔE → 折算 Ri”串在一个 for 循环里前后 14 个色样一次跑完。R zeros(14, 1); for i 1:14 % 色样反射率插值到公共波长网格 rho interp1(wl_ref, refMat(:, i), wl, pchip, 0); % 测试光源下该色样的 XYZ Xi trapz(wl, spd .* rho .* xbar); Yi trapz(wl, spd .* rho .* ybar); Zi trapz(wl, spd .* rho .* zbar); % 参考光源下该色样的 XYZ Xi_r trapz(wl, spd_ref .* rho .* xbar); Yi_r trapz(wl, spd_ref .* rho .* ybar); Zi_r trapz(wl, spd_ref .* rho .* zbar); % 用测试光源下的白点做色适应换算成参考光源下的等效XYZ adapted_XYZ apply_vonkries([Xi; Yi; Zi], ... test_white, ref_white); % CIE 1964 W*U*V* 色差计算 dE delta_e_wuv([Xi_r; Yi_r; Zi_r], adapted_XYZ); % 单个色样显色指数 R(i) 100 - 4.6 * dE; R(i) max(R(i), 0); % CRI 标准下限截断为0 end % 一般显色指数取前8个平均值 Ra mean(R(1:8));代码里apply_vonkries和delta_e_wuv是自定义函数。前者做的是第 3 章提到的三通道缩放后者则要先计算 WUV* 坐标。Ri 100 - 4.6 × ΔE是显色指数评分的核心关系式4.6 是 CIE 在制定标准时定下的比例系数。我不建议把它当可调参数处理实际程序里一旦改了这个系数Ra 就不再是标准意义上的 Ra。4.3 输出图表并整理成能进 PPT 的光谱报告计算跑完光给一个 Ra 数字在汇报时不够直观。我一般会生成一张两联图上面是测试光源 SPD 和参考光源 SPD 的对比下面是 14 个 Ri 的柱状图。这样放到 PPT 里只要截个图就能解释清楚“为什么 R9 偏低”。figure; subplot(2,1,1); plot(wl, spd, k, LineWidth, 1.2); hold on; plot(wl, spd_ref, r--, LineWidth, 1.2); xlabel(Wavelength (nm)); ylabel(Relative Power); legend(Test Source, Reference Source); subplot(2,1,2); bar(1:14, R, FaceColor, [0.7 0.2 0.2]); xlabel(Test Color Sample); ylabel(Ri); ylim([0 100]); grid on; exportgraphics(gcf, CRI_report.png, Resolution, 300);exportgraphics是 MATLAB R2020a 以后更推荐的方式比老式的print或saveas更稳不会出现字体错位。如果一定要直接生成可编辑的 PPT也可以用mlreportgen.ppw创建幻灯片但大多数情况下报告里放一张高分辨率 PNG 就够了。把这张图和上面的表一起放进 PPT就是客户能看懂的完整显色性 CRI 计算结论。5. 用黑体回归测试验证 CRI 计算程序边界与误差修正写完显色性 CRI 计算程序后第一件事不是拿实测 LED 数据跑而是先做一个“已知答案”的回归测试。最常见的验证手段是拿一个黑体 SPD 输入程序因为黑体的 Ra 理论上应该非常接近 100一般认为大于 99 就算正常。% 生成 3200K 黑体光谱跑一遍主程序 spd_test blackbody_spd(3200, wl); % 归一化到Y100和主程序接口保持一致 spd_test spd_test / trapz(wl, spd_test .* ybar) * 100; % 参考光源也用3200K黑体 spd_ref blackbody_spd(3200, wl); spd_ref spd_ref / trapz(wl, spd_ref .* ybar) * 100; % 调用主循环得到的Ra应当接近100 Ra_test compute_cri(spd_test, spd_ref);这个测试有一个好处测试光源和参考光源都是同一个黑体公式生成适应矩阵、积分方向若有系统性偏差会在色差上抵消掉一大半因此 Ra 偏低反而说明问题不在光源模型而在样本积分或评分公式。我遇到过的情况是插值函数误用了linear窄色样反射率曲线在转折处被抹平导致 Ra 掉到 97换成pchip后回到 99.9。这就是回归测试的价值。除了黑体回归还有三个边界值得一测。第一个是低色温光源比如 2500K 琥珀色灯CCT 低于黑体轨迹最低可计算区间参考光源选取会不稳定建议在程序里做色温下限检查低于 2000K 直接给出警告。第二个是窄带 LED这类 SPD 在 450nm 蓝光和窄峰红光处极度不连续插值步长太大会把峰值位置移动几纳米直接影响色度坐标这时要对比 1nm 和 0.5nm 重采样的结果差异超过 0.5 就说明仪器原始数据步长不够。第三个是光谱尾部噪声780nm 附近仪器响应弱数据抖动会被积分放大处理办法是取 750–780nm 均值做平滑不要对最后一两个点做尖锐截断。把这些检查项固化成测试脚本每次改代码后自动重跑一遍。最直接的做法是把compute_cri封装成函数写一个test_cri.m依次调用黑体测试、窄带模拟和边界色温测试任何一项不通过终端就会报错。这套验证流程比反复拿实测数据人工核对更省时间也是显色性 CRI 计算程序能从“跑通”走到“可信”的关键一步。本文还有配套的精品资源点击获取