简介：基于MATLAB的拉曼光谱处理示例代码，面向化学、材料、生物等领域需要分析光谱数据的科研人员与工程师。资源针对光谱数据中常见背景噪声与共线性问题，演示了airPLS偏最小二乘变体算法的完整实现流程，涵盖数据导入、预处理、模型构建、验证与可视化等关键环节。包内共4个文件，包含两个m格式脚本（实现airPLS算法及主测试程序）、一个xlsx原始光谱数据表和一个pptx流程演示文稿，压缩包仅260KB，轻量易用。目前已有1554人学习下载，适合想要快速上手拉曼光谱数据处理、深入理解airPLS算法原理并提升MATLAB编程能力的读者。通过实际运行代码，可掌握从原始数据到分析模型的全链路操作，并依据演示文稿理解每一步设计思路，为后续拓展到更复杂的光谱分析任务打下基础。

1. 这组MATLAB代码解决什么问题：拉曼光谱荧光背景拦路，airPLS直接剥干净

拿到拉曼光谱数据，最先看到的往往不是样品峰，而是一条被荧光背景顶起来的丘陵。峰和背景纠缠在一起，无论做峰强统计还是峰位识别，结果都不靠谱。用这组基于MATLAB的拉曼光谱处理代码，配合airPLS自适应迭代加权惩罚最小二乘算法，能把这条丘陵状基线剥离出来，得到干净的校正光谱。这份资源自带原始数据，路径、格式、调用方式都是现成的，适合正在做拉曼光谱预处理和荧光背景扣除的科研或工程开发人员，也适合想快速看懂airPLS算法如何落地成MATLAB函数的人。解压之后跑一遍demo脚本，三秒钟就能看到原始光谱、拟合基线、校正光谱三条曲线。

2. 项目结构与数据组织：动手跑之前先把家底摸清

拿到资源后，我习惯先不急着运行，把目录翻一遍，把每个文件打开扫一眼。这个习惯帮我少踩了不少坑——不了解输入输出格式的算法代码，跑起来容易变成碰运气。这套项目的结构非常直白，没有复杂的工程依赖，也不需要额外的MATLAB工具箱，对新手相当友好。

2.1 文件清单与职责划分

解压之后会看到几个核心文件，职责划分得清清楚楚：

airPLS.m 是整个项目唯一的算法函数，demo_airPLS.m 则把数据取出来喂给它。这种组织方式看得出设计意图：算法函数保持纯函数形态，只接收数值、返回数值，不掺入文件读取和绘图逻辑。数据文件和代码放在一起，好处是拿过来就能直接验证效果，不用满世界找测试数据。对想二次开发的人来说，这种单入口、单函数的组织也是最容易下手的结构。

2.2 原始数据的格式与读取细节

data/sample1.txt 内容非常朴素，每行两个数字，空格分隔。第一列是拉曼位移，单位通常是 cm⁻¹，第二列是信号强度，单位是 a.u. 之类。文件开头几行长这样：

800.00 120.50 802.00 121.30 804.00 118.90 806.00 125.70 808.00 122.80

用 MATLAB 的 load 函数可以直接读进来，因为这种纯数值文件没有表头，load 会把整个文件读取为一个 N 行 2 列的矩阵。读取之后拆开两列，一列当横轴，一列当纵轴。实际代码只需要三行：

data = load('data/sample1.txt'); shift = data(:, 1); intensity = data(:, 2);

这里有一个值得注意的细节：airPLS 处理的是单条光谱的强度向量，不是整个数据矩阵，所以真正传给算法的只有 intensity 这一列。shift 列只用于画图时定位横坐标。如果文件换成了带表头的 CSV，load 就不好使了，我一般改用 readtable 或者 textscan 处理，但这份资源里的文本文件不需要额外折腾。

读取数据后，建议先用 size(data) 确认一下矩阵形状。偶尔会遇到文件末尾多一个空行，导致 data 末尾出现 0 行或 NaN，如果直接拿去跑算法，结果会莫名多出一截曲线，检查数据是最快的排查手段。

2.3 演示脚本的主流程调用链

demo_airPLS.m 做的事情可以概括成四步：读数据、调算法、计算校正结果、画图。核心调用链在一个代码块里就能说清楚：

% demo_airPLS.m 主流程 % 第一步：清理工作区，避免旧变量干扰 clear; close all; clc; % 第二步：加载原始拉曼光谱 data = load('data/sample1.txt'); shift = data(:, 1); intensity = data(:, 2); % 第三步：调用 airPLS，lambda 取 1e5，二阶差分，最多迭代 20 次 [baseline, corrected] = airPLS(intensity, 1e5, 2, 20); % 第四步：可视化对比 figure('Position', [100 100 800 600]); plot(shift, intensity, 'b-', 'LineWidth', 1); hold on; plot(shift, baseline, 'r-', 'LineWidth', 2); plot(shift, corrected, 'g-', 'LineWidth', 1); xlabel('拉曼位移 / cm^{-1}'); ylabel('强度 / a.u.'); legend('原始光谱', 'airPLS基线', '校正后光谱');

这一步里 subtle 的关键点在于 corrected 是函数内部算好再返回来的，不是拿到 baseline 之后在脚本里做减法。也就是说 airPLS.m 返回值里已经包含了校正光谱，调用者少写一行是一行，也避免了“忘了扣基线”这种低级失误。

2.4 为什么把算法函数和演示脚本分开

把算法封装成独立函数，演示脚本单独存在，看起来是微不足道的工程习惯，但对二次开发影响不小。如果所有代码都堆在一个 main 脚本里，换个数据文件就得打开脚本改路径，调参数也得小心翼翼找位置。现在这种结构下，想换数据只需要改 demo_airPLS.m 的 load 路径，想改算法只需要动 airPLS.m，互不干扰。后面做批量处理时，这个习惯的回报会更加明显——直接循环调用 airPLS 函数就行，完全不用碰脚本内部逻辑。

3. airPLS核心算法拆解：惩罚最小二乘与自适应权重的配合

这一章是重点。网上讲 airPLS 原理的帖子不少，但大多停在公式层面，很少有文章把 MATLAB 实现一点一点掰开讲。这里我用项目里的 airPLS.m 作为蓝本，从数学骨架到代码细节逐段展开，看完之后你就能自己上手调整这个算法。

3.1 先理解惩罚最小二乘的基本骨架

airPLS 的全称是 Adaptive Iteratively Reweighted Penalized Least Squares，核心思想是在惩罚最小二乘的框架上加入自适应权重。它的目标函数可以写成这样：

min ||w .* (y - z)||^2 + lambda * ||D * z||^2

其中 y 是原始光谱强度，z 是需要拟合的基线，w 是每个数据点的权重向量，D 是差分矩阵，lambda 是平滑参数。这个式子的直观含义是：左边一项让基线 z 尽可能贴近原始光谱，右边一项让基线的差分不要太大，也就是保证基线足够平滑。两者之间由 lambda 平衡。

D 矩阵的作用是计算基线的 order 阶差分。order 取 2 时，D 生成的就是二阶差分，相当于约束基线二阶导不能太大，这保证了基线是一条光滑的缓变曲线，不会出现锯齿状的剧烈抖动。lambda 越大，右边惩罚越重，基线越平滑，但可能跟不上光谱的缓慢变化；lambda 越小，基线越能贴合细节，但也可能钻进噪声里。

这个优化问题是一个标准的最小二乘问题，求导后得到一个线性方程组，直接求解即可。关键部分在于权重 w 如何取值，这决定了算法能不能把峰区域巧妙避开。

3.2 自适应权重是怎么一回事

airPLS 的巧妙之处在于权重不是一次性算好的，而是在迭代过程中不断更新。刚开始时所有点权重一致，拟合出一条大致形状的基线，然后计算残差 d = y - z。此时如果某个点的光谱强度高于基线，也就是残差为正，这个点很可能是样品峰的位置，权重会被置为 0，意味着这个点不参与下一轮基线拟合。反之，如果残差为负，说明这个点落在基线下方，是背景区域，权重会被指数放大，放大程度与残差绝对值成正比。

这里面的逻辑是：峰区域的点对基线拟合的“拉力”被直接屏蔽，基线就不会被峰顶带跑，而是尽可能逼近除峰区域以外的底部。负残差点权重越大，下一轮基线就越贴近这些点，于是基线在迭代中一点点向光谱底部沉降。迭代越多，基线越贴近真实背景。这种策略让 airPLS 不需要预先标定峰位置，纯粹靠残差正负号来自适应，非常省事。

3.3 airPLS.m 核心代码逐段讲解

项目里的 airPLS.m 我做了简化整理，核心逻辑保持原样，只增加了一些注释和防溢出处理。完整代码如下：

function [z, y_corrected] = airPLS(y, lambda, order, itermax) % airPLS - 自适应迭代加权惩罚最小二乘基线校正 % 输入: % y : 原始光谱强度，N×1 列向量 % lambda : 平滑参数，越大基线越平滑，典型值 1e2 ~ 1e9 % order : 差分阶数，通常取 2 % itermax : 最大迭代次数，通常 10 ~ 20 % 输出: % z : 拟合出的基线 % y_corrected : 校正后的光谱 = y - z y = y(:); N = length(y); % 构建 order 阶差分矩阵 D，DD = D'*D 用于惩罚项 D = diff(speye(N), order); DD = D' * D; % 初始权重所有点都设成 1 w = ones(N, 1); % 记录上一次迭代的残差，用于收敛判断 d_prev = zeros(N, 1); for iter = 1:itermax % 将权重向量转成稀疏对角矩阵 W W = spdiags(w, 0, N, N); % 求解加权惩罚最小二乘： (W + lambda * DD) * z = w .* y z = (W + lambda * DD) \ (w .* y); % 残差：原始信号减去基线 d = y - z; % 如果所有残差都大于等于 0，基线已在数据之下，无需再迭代 if ~any(d < 0) break; end % 收集负残差点，这些是背景区域 neg_idx = d < 0; dn_mean = mean(abs(d(neg_idx))); % 防止均值过小导致指数爆炸 if dn_mean < eps dn_mean = eps; end % 更新权重：负残差点权重按指数放大，峰区域权重保持 0 w_new = zeros(N, 1); exp_arg = iter .* abs(d(neg_idx)) ./ dn_mean; exp_arg = min(exp_arg, 700); % 防止 exp 溢出 w_new(neg_idx) = exp(exp_arg); w = w_new; % 收敛判断：残差变化量相对于当前残差范数足够小 if norm(d - d_prev) / (norm(d) + eps) < 1e-6 break; end d_prev = d; end % 校正光谱 = 原始信号 - 基线 y_corrected = y - z; end

这段代码里，第一步是构造差分矩阵 D，用 diff(speye(N), order) 一行生成稀疏矩阵，然后累乘得到 DD。这里必须用稀疏矩阵，如果 N 是几万点，稠密矩阵直接能把内存吃光。第二步是迭代求解稀疏线性系统，MATLAB 的反斜杠运算符会自动选择算法，不需要手写求解器。权重更新是整个算法的灵魂：正残差点权重为 0，负残差点权重指数增长，迭代次数 iter 乘在指数里，所以越到后期负残差区域拟合越逼近原始信号。末尾的收敛判断比较当前残差和上一次残差的范数比值，小于 1e-6 就认为基线稳定了，提前退出循环节省时间。

代码里的 min(exp_arg, 700) 是工程上的一个防溢出处理，严格按论文公式推的话指数值可能大到让 exp 返回 Inf，MATLAB 里出现 Inf 会导致后面的矩阵求解全部变成 NaN。我把指数上限截断到 700，这是 exp 函数在双精度浮点下的安全边界，不影响实际效果，但彻底规避了数值爆炸的隐患。初次接触这段代码的人可能觉得这是多此一举，真跑一次极端数据就明白它的价值了。

3.4 参数含义和调节方向

airPLS 需要关注的参数不多，但每一个都对结果影响很大：

调参的顺序我一般是先固定 order=2 和 itermax=20，只动 lambda。如果校正后的光谱在峰区出现下凹，说明 lambda 太大，基线弯不进去，把峰底也兜住了；如果基线明显粗糙、校正谱里带毛刺，说明 lambda 太小，基线在拟合噪声。以 sample1.txt 这条数据来说，1e5 属于比较稳妥的中等值，直接跑就能得到不错的基线。

4. 从零跑通：环境准备、运行流程与结果验证

前面把原理讲清楚了，这一章进入实操环节。目标非常简单：拿到这份资源，在MATLAB里把它跑通，并且确认结果不是“看起来好看”而是真的正确。

4.1 运行环境与准备

airPLS.m 的核心操作是稀疏矩阵构建和线性求解，不依赖任何付费工具箱，所以环境要求非常低。只要是 MATLAB 软件本体能正常运行，版本在 R2016a 或更新，就可以直接跑。操作系统方面，Windows、macOS、Linux 下的 MATLAB 都没有问题，因为代码本身不涉及系统调用。

拿到项目后，第一步是把解压后的文件夹放到一个全英文路径下，避免某些系统环境下中文路径导致 load 读取失败。然后把当前工作目录切换到文件夹位置，可以直接用 cd 命令，或者在编辑器中打开 demo_airPLS.m，右键运行，MATLAB 会自动切换目录。顺手在命令行窗口执行一下 which airPLS，确认函数能被搜索到，这一步能避免后续“找不函数”的低级报错。

4.2 完整运行流程

运行 demo_airPLS.m 是唯一需要手动做的事。脚本执行完毕，工作区会出现几个核心变量：shift、intensity、baseline、corrected。此时可以用 whos 看变量尺寸，如果数据是几百到几千个点的量级，整个过程应该在一两秒内完成。从加载数据到绘图，完整流程在一个代码块里约束得明明白白：

% 4.2 另一种手动分步执行法，方便观察中间结果 clear; close all; clc; % 读取第二份样本来验证算法的泛化能力 data = load('data/sample2.txt'); shift2 = data(:, 1); intensity2 = data(:, 2); % 调用 airPLS [baseline2, corrected2] = airPLS(intensity2, 1e5, 2, 20); % 画出结果，并把校正谱向上平移一段距离便于观察 figure; plot(shift2, intensity2, 'b-'); hold on; plot(shift2, baseline2, 'r-', 'LineWidth', 2); plot(shift2, corrected2 + 500, 'g-'); xlabel('拉曼位移 / cm^{-1}'); ylabel('强度 / a.u.'); legend('原始光谱', 'airPLS基线', '校正后光谱+500'); grid on;

这段代码和 demo_airPLS.m 功能一致，只是换成了第二份数据，并且在校正光谱上加了 500 的偏移。加偏移纯粹是为了视觉上不跟原始光谱叠在一起，不是数据处理步骤。如果你发现校正后的曲线起起伏伏像噪声，先检查一下是不是忘了这个偏移量。

或者，在命令行窗口直接用 readmatrix 之类的其他读取函数替换 load 也可以，但 load 在这里是最省事的。执行完后，图形窗口会弹出三条曲线，蓝线是原始数据，红线是拟合基线，绿线是校正后的光谱。

4.3 结果验证的几个硬指标

很多人跑完图形窗口弹出，看一眼“哇，基线挺顺的”就完事了，很容易被视觉误差骗过。我一般会用三个维度验证结果是不是靠谱。

第一个维度是基线形态。拟合出的红线应该是一条平滑的缓变曲线，它会贴着所有谱峰的底部穿过，但不能有任何峰形结构。如果红色基线上鼓出一个包，说明有一个峰被基线吃掉了，lambda 需要调大或调小，具体看鼓包形态。第二个维度是校正后光谱的零点位置。校正谱在无峰区间应该贴着零线附近摆动，幅度和噪声水平一致，不应出现大段的负基线偏移。第三个维度是峰形完整性。校正后的峰不应该出现明显变形，尤其是峰尾不能出现下凹的“塌陷”，这种塌陷通常是基线在峰区过度弯曲造成的。

如果三个维度都通过，说明这次校正基本可信。如果只追求“看起来平”，忽略了峰形是否真实，后面做定量分析时会埋下大雷。把校正结果用 dlmwrite 存成文本，方便后续导入其他软件：

% 保存校正结果到当前文件夹 dlmwrite('sample1_corrected.txt', [shift, corrected], 'delimiter', '\t', 'precision', 6);

这样导出的是两列制表符分隔文本，第一列是位移，第二列是校正后强度，跟原始数据的格式保持一致，后续做峰拟合或者数据共享都不用再转换格式。

5. 避坑与常见问题排查：airPLS实践里的四个翻车现场

airPLS 在实际使用中翻车概率不算低，几乎每个人都会遇到至少一次莫名其妙的基线。我把常见的坑整理成一个速查表，然后挑四个最典型的场景展开分析。

5.1 高频问题速查表

速查表适合应急，下面几个案例能帮你理解问题背后的机制，下次再遇到不用翻表格也能直接判断。

5.2 翻车一：强峰区基线“上浮”，校正谱出现凹陷

现象是校正后的光谱在一个较强峰的位置不是凸起的峰形，反而凹下去一块，峰顶被削掉一层。我第一次遇到时以为是算法程序有 bug，折腾了很久。

原因在 lambda 上。lambda 取太大时，约束项惩罚过重，基线变得非常僵硬，在峰区域弯曲不下来，只能从峰的半腰处斜穿过去。此时在峰顶附近，原始强度高于基线，但峰的中低区域被基线兜住了，校正之后两边的肩部被扣除过度，就形成了凹陷。

解决办法是先把 lambda 从 1e5 降到 1e3 或者 1e4，观察凹陷是否消失。如果消失，说明基线柔性够了；如果基线开始出现锯齿，说明降过头了，再用二分法找出中间值。我现在的习惯是先看一眼峰宽度和整体信噪比，峰宽度大的数据用偏大的 lambda，窄峰数据用偏小的 lambda，比固定值盲调高效得多。

5.3 翻车二：校正谱整体出现负基线偏移

现象是校正后的光谱所有峰都立起来了，但无峰区间的基线不在零附近，而是整体趴在负几百的位置。这会让后续的峰面积统计失真。

原因是收敛判据过松或者负残差均值被极端值干扰。负残差区域如果混入一个很大的离群值，dn_mean 会被拉大，导致权重普遍偏小，基线下降速度变慢。迭代没到真正收敛就提前退出，返回的基线比真实背景高出一截，校正谱自然整体偏负。

这种情况先检查 itermax，把它从 20 调到 50 再跑一次。如果问题消失，说明之前确实是提前跳出循环。如果还在，检查数据尾部是否有一段异常高的谱段，比如激光残影导致强度突然升高，这会干扰均值统计。把异常段截掉再处理，通常能解决问题。

5.4 翻车三：长光谱处理时内存溢出或运算缓慢

现象是光谱数据有几万甚至十几万个点，跑 airPLS 时 MATLAB 直接报错内存不足，或者运行特别慢，连画图都卡。

原因是构建差分矩阵时，diff(speye(N), order) 会生成一个 N 阶的稀疏矩阵，当 N 超过几万时，DD = D'*D 的存储量急剧增加，即便稀疏也扛不住。而且每次迭代都要重新构建对角矩阵和求解线性系统，时间复杂度接近 O(N^3)，数据越长越吃力。

解决思路有几个。最直接的是下采样，把光谱从 100000 点抽稀到 20000 点，跑完基线再插值回原长度，峰位不会偏。也可以用分段处理，把光谱切成几段分别跑 airPLS，段与段之间保证一定重叠，再把基线拼接起来。如果频繁处理大数据，建议把下采样方法封装进预处理函数里，一劳永逸。

5.5 翻车四：换了一个数据文件，同样的参数结果完全乱套

现象是 sample1 调好的参数，搬到 sample2 上，基线要么飞到天上，要么压到地板下，完全没有可用性。

原因是不同样本强度量级不同。sample1 的强度可能在几十到几百之间，sample2 如果动辄几千，同一 lambda 值对这两个量级的“平滑力度”完全不是一个概念。lambda 的绝对值没有跨样本的普适性，它本质上是相对于信号幅度的相对量。

解决办法是先做一次归一化，把信号强度压缩到 0 到 1 之间再跑 airPLS，这样 lambda 就有了跨样本的稳定含义。归一化可以简单写成 intensity / max(intensity)，跑完之后 baseline 和 corrected 都在原尺度上，画图时再乘回去即可。

6. 二次开发：把airPLS变成自己的光谱预处理管道

理解了算法原理之后，二次开发的重点就不再是单条光谱，而是如何把它接进自己的处理流程里。

6.1 批量处理整个文件夹的光谱

实际项目中很少只处理一条光谱，通常一个文件夹里有几十条。批量循环调用 airPLS 就行了，把数据读取、校正、保存封装成一个循环：

% batch_airPLS.m 批量处理 data 目录下所有 txt 文件 files = dir('data/*.txt'); for i = 1:length(files) data = load(fullfile(files(i).folder, files(i).name)); shift = data(:, 1); intensity = data(:, 2); % 归一化到 0~1，保证 lambda 对所有样本语义一致 norm_intensity = intensity / max(intensity); [~, corrected] = airPLS(norm_intensity, 1e5, 2, 20); corrected = corrected * max(intensity); % 输出文件加 _corrected 后缀 out_file = fullfile(files(i).folder, ... [files(i).name(1:end-4) '_corrected.txt']); dlmwrite(out_file, [shift, corrected], 'delimiter', '\t'); end

这段代码里我对每条谱都先归一化再调 airPLS，跑完再乘回原尺度。这样做的好处是 lambda 的设置对所有样本都通用，不会因为样本间强度差异导致参数失效。输出文件和原始文件一一对应，后续随便拿哪个出来都能直接画图。

6.2 让lambda自动适应数据

固定 lambda 总归是偷懒的做法，碰到形态差异大的数据集还是得手动调。一个简单的启发式办法是用光谱峰面积的大致估计来缩放 lambda，比如先算一下信号绝对值的和，再让 lambda 包一层比例系数。这种自动缩放不求精确，只求能把参数限定在一个不需要反复试的数量级范围内。

我个人的教训是：以前总想找一组万能参数，一个 lambda 打天下，结果换了几批实验数据后频频翻车。从那以后，每组数据我都强制自己先跑一遍归一化，再根据峰宽量级选一个 lambda 的起点值，看一次基线拟合效果再往下微调。这个流程多花一分钟，但省下来的调试时间远不止十分钟。希望这个习惯也能帮到你。

本文还有配套的精品资源，点击获取