简介：面向MATLAB毕业设计的Voigt模型拟合项目，以简洁的代码与文档展示如何通过自定义Voigt函数结合最小二乘思路对光谱类数据做谱形拟合。Voigt模型融合Lorentzian与Gaussian分布，适用于光谱学、核磁共振和声学等宽峰信号分析场景，项目重点给出拟合主程序与模型函数，并配套演示脚本、说明文档及示例效果图，便于从算法到应用快速理解。压缩包共8个文件，以3个m源码文件为主，辅以md说明、txt文本、png样图及Git配置信息，整体仅31KB，结构轻量；目前已有98人学习，适合需要完成相关课题或补强MATLAB数值拟合能力的本科生作为参考。通过阅读主程序与演示脚本，可掌握Voigt模型建模、参数拟合与结果可视化整个流程，也可迁移至其他线型拟合任务。
1. 从光谱峰形到Voigt拟合：这个MATLAB毕业设计在解决什么
实际光谱分析中，我们测到的峰形往往既不是纯高斯也不是纯洛伦兹，而是两者卷积后的形态。Voigt模型能同时描述多普勒展宽带来的高斯成分和碰撞展宽带来的洛伦兹成分，在红外光谱、核磁共振、X射线衍射等领域非常常用。这个zip包里装的正是一套用MATLAB实现Voigt模型拟合的毕业设计源码：myvoigt.m负责计算Voigt分布曲线，voigtfit.m是拟合主函数，voigtfit_demo.m告诉你如何把两者串起来，最后还会输出一张spectrumfit.png拟合结果图。如果你正在准备“matlab 毕业设计”题目，或者手头有需要分离高斯-洛伦兹贡献的光谱数据，这套代码值得逐行拆解。
2. Voigt模型的数学本质与myvoigt.m的两种实现路径
2.1 Voigt分布是Gaussian与Lorentzian的卷积，不是简单加和
Voigt分布的定义是高斯分布与洛伦兹分布的卷积积分，数学上写成：
V(x) = ∫ G(x') L(x - x') dx'
其中G是高斯函数，L是洛伦兹函数。这个积分没有封闭的解析表达式，只能通过数值积分或近似公式来计算。很多初学者会直接把高斯和洛伦兹曲线做线性叠加，但那只是pseudo-Voigt，不是真正的Voigt。在工程上两种做法都有应用，但原理和误差边界完全不同，需要根据拟合的物理场景来选择。
2.2 数值积分实现：基于quadgk的精确卷积
myvoigt.m中一种直接实现方式是利用MATLAB的integral函数做无穷积分。下面这段代码展示了最朴素的写法：
function y = myvoigt(x, sigma, gamma) % 数值积分实现的Voigt函数 % x: 频率或波长偏移，sigma: 高斯展宽参数，gamma: 洛伦兹展宽参数 y = zeros(size(x)); for i = 1:length(x) f = @(t) exp(-t.^2 / (2 * sigma^2)) ./ ... (1 + ((x(i) - t) ./ gamma).^2); y(i) = (1 / (sqrt(2 * pi) * sigma)) * ... (gamma / pi) * integral(f, -Inf, Inf); end end这个实现的好处是概念上严格对应数学定义，适合用来做理论校验。但实际运行时会发现它非常慢，原因是integral需要多次迭代求积分，而且for循环遍历每个x点，整体复杂度很高。如果光谱数据有上千个点，拟合一次可能要数十秒甚至更久，这在批量拟合场景中不可接受。另外，当sigma或gamma趋近于0时，被积函数变化过于尖锐，数值积分可能不收敛。
2.3 工程近似：pseudo-Voigt加权组合
更常用的做法是pseudo-Voigt模型，它用高斯和洛伦兹的加权和来近似Voigt。权重参数eta控制两者比例，eta=0时是纯高斯，eta=1时是纯洛伦兹。这个近似在拟合精度要求不高时误差通常在2%以内，计算速度却快三个数量级。
function y = myvoigt(x, sigma, gamma, eta) % 伪Voigt近似：Gaussian与Lorentzian线性加权 % eta 为0时是纯Gaussian，为1时是纯Lorentzian % 注意：这里未包含幅值系数，调用时可在外部乘amplitude g = exp(-(x.^2) / (2 * sigma^2)); l = 1 ./ (1 + (x ./ gamma).^2); y = (1 - eta) .* g + eta .* l; end这段代码中的sigma对应高斯宽度，gamma对应洛伦兹半宽，eta做插值。实际拟合时，为了让曲线峰值与数据匹配，通常会在外层乘以一个幅值参数。值得注意的是，pseudo-Voigt的eta与真实Voigt的sigma/gamma之间存在一个经验关系，更严谨的做法是使用Olivero和Longbothum提出的近似公式，但多数光谱拟合场景直接让eta自由拟合，也能得到可接受的曲线形状。
2.4 两种实现的精度与耗时对比表
| 实现方式 | 理论精度 | 单次计算相对耗时 | 适用场景 |
|---|---|---|---|
| 数值积分 | 高 | 约千倍 | 理论验证、生成仿真数据 |
| pseudo-Voigt近似 | 较高 | 极快 | 光谱拟合、批量处理、实时分析 |
从voigtfit_demo.m的实际运行效果看，默认使用的应该是pseudo-Voigt近似，因为spectrumfit.png中拟合曲线平滑且计算很快。若你在自己的数据上发现拟合残差呈现系统性的S形偏差，可以尝试切回数值积分版本，对比两种实现的结果来检查是否近似误差导致的。
3. voigtfit.m 核心拟合算法与MATLAB优化工具箱对接
3.1 函数接口设计与调用约定
voigtfit.m是整套代码的拟合入口。按照MATLAB函数文件惯例，它的调用方式通常类似：
[param, resnorm, residual] = voigtfit(x, y, initParams, lb, ub)参数含义分别为：x是自变量向量，y是观测光谱强度，initParams是模型参数的初始猜测，lb和ub是参数下界和上界。返回的param是拟合得到的最优参数向量，resnorm是残差平方和。如果调用时没有提供lb和ub，则内部默认使用足够宽松的边界，例如所有参数的下界为0，上界为数据最大值的10倍。这样的设计让初级用户可以直接运行demo，但做正式实验时建议显式传入边界。
3.2 用lsqcurvefit做非线性最小二乘拟合
MATLAB的Optimization Toolbox中，lsqcurvefit是求解这类非线性曲线拟合的标准函数。voigtfit.m内部的常见实现方式是把myvoigt包成一个模型函数句柄，然后交给lsqcurvefit迭代优化。下面是一段核心代码框架：
function [param, resnorm] = voigtfit(x, y, init, lb, ub) % 定义带幅值的Voigt模型函数句柄 model = @(p, xdata) myvoigt(xdata, p(1), p(2), p(3)) .* p(4); % 配置优化选项：显示迭代信息，设置高精度和最大迭代次数 opts = optimoptions('lsqcurvefit', ... 'Display', 'final', ... 'TolFun', 1e-10, ... 'TolX', 1e-10, ... 'MaxIter', 2000); % 执行拟合 [param, resnorm] = lsqcurvefit(model, init, x, y, lb, ub, opts); end这里模型句柄中p(1)是sigma，p(2)是gamma，p(3)是eta，p(4)是幅值amplitude。之所以把幅值单独放在myvoigt外面，是因为myvoigt本身体现的是归一化形状，而实际光谱的强度是任意单位的，必须有一个幅值缩放项。注意lsqcurvefit要求x和y都是列向量，如果传入行向量，MATLAB不会报错但会得到错误结果，这是常见坑点之一。
3.3 参数初值、上下界与归一化策略
初始值设置得不好，lsqcurvefit很容易收敛到局部极小值。对于Voigt拟合，合理的初值估计方法是：先找到数据中的最大峰值位置和峰值强度，用峰值强度作为幅值的初值；然后测量半高宽FWHM，由于FWHM与sigma、gamma的关系不是线性的，可以设置sigma=FWHM/2.354，gamma=FWHM/2，eta=0.5。如果数据中噪声很强，可以先用smooth函数做一次平滑再计算FWHM。
在myvoigt的x不包含中心位置时，需要先将数据以峰位为中心进行平移，也就是把x减去峰位后再传入拟合函数。若希望拟合过程中峰位也能自由调整，就需要在模型函数里加入center参数，例如：
model = @(p, xdata) myvoigt(xdata - p(5), p(1), p(2), p(3)) .* p(4);其中p(5)为峰位。此时初值p(5)设为最大峰对应的x坐标，上下界则根据数据范围设定，比如峰位只允许在数据最左侧和最右侧之间变化。边界设置的原则是既不能太紧导致正确解在边界外，也不能太松导致算法搜索过大区间而浪费迭代。
3.4 拟合质量指标：SSE、R²与残差分析
lsqcurvefit返回的resnorm就是SSE，但单看SSE不能判断拟合是否良好。更常用的指标是决定系数R²，计算方式如下：
SSE = resnorm; SST = sum((y - mean(y)).^2); R2 = 1 - SSE / SST;R²越接近1表示拟合越好，但要注意R²高不代表模型正确。如果残差residual = y - model(param, x)存在明显的周期性或形状，说明模型缺少基线项。很多实际光谱数据会有线性漂移，需要在模型里增加一个一次项：model = @(p, xdata) myvoigt(xdata - p(5), p(1), p(2), p(3)) .* p(4) + p(6) + p(7) * xdata;这样扩展后，p(6)是常数基线，p(7)是斜率。添加基线项后，SSE通常会显著下降，但要防止过拟合，可以用赤池信息准则AIC来权衡。
4. voigtfit_demo.m 实战演示：从数据到拟合并保存spectrumfit.png
4.1 生成模拟Voigt谱线数据
voigtfit_demo.m的第一步通常是用已知参数模拟一组光谱数据，用来验证拟合代码的正确性。由于真实数据的真值是未知的，模拟数据可以让我们对比拟合结果与设定值的偏差。常见做法如下：
x = linspace(-8, 8, 600); sigma = 0.8; gamma = 0.6; eta = 0.4; amp = 3.0; y_clean = myvoigt(x, sigma, gamma, eta) * amp; rng(2024); % 固定随机种子，使结果可重现 y_noisy = y_clean + 0.1 * randn(size(x));这里x取[-8,8]共600点，保证峰两侧有足够长的尾部。噪声标准差设为0.1，模拟真实测量中的随机误差。固定随机种子后每次运行结果一致，方便在答辩中演示时保证可复现性。
4.2 调用voigtfit.m完成拟合
有了模拟数据后，直接调用voigtfit.m：
init = [1.0, 0.5, 0.5, max(y_noisy)]; lb = [0.01, 0.01, 0, 0]; ub = [5, 5, 1, 20]; p = voigtfit(x, y_noisy, init, lb, ub);这里初值选择依据是：sigma和gamma分别取1.0和0.5，eta取0.5为中性值，幅值用数据最大值。经过几次迭代后，p应当收敛到接近[0.8, 0.6, 0.4, 3.0]的值。如果发现p与设定值相差超过20%，说明算法可能陷入了局部极小，或者数据点数不够。此时可以尝试将初值改为sigma=0.7, gamma=0.7, eta=0.3，看是否收敛到不同结果。
4.3 绘制并保存结果图
拟合完成后，demo脚本会绘制原始数据点和拟合曲线，并保存为spectrumfit.png。绘图代码通常为：
y_fit = myvoigt(x, p(1), p(2), p(3)) * p(4); figure('Color', 'w'); plot(x, y_noisy, 'ob', 'MarkerSize', 4, 'LineWidth', 1); hold on; plot(x, y_fit, '-r', 'LineWidth', 2); legend('原始数据', 'Voigt拟合', 'Location', 'north'); xlabel('相对位置 (cm^{-1})'); ylabel('强度 (a.u.)'); grid on; saveas(gcf, 'spectrumfit.png');这一段代码涉及MATLAB绘图的基本操作。plot的第一个参数是x，第二个是y，第三个是样式字符串，'ob'表示蓝色圆圈标记，'-r'表示红色实线。legend的Location设为north，意思是图例放在顶部中间。saveas支持png、pdf、eps等格式，写论文时建议输出矢量图：saveas(gcf, 'spectrumfit', 'epsc')。绘图是毕业设计中的加分项，清晰美观的图比参数表格更直观。
4.4 对demo脚本做参数敏感性测试
为了让答辩更有说服力，可以修改demo脚本做一个简单的敏感性分析：把初值在真实值附近上下浮动50%，循环100次调用voigtfit，记录每次的收敛结果。如果100次中绝大多数都收敛到同一组参数，说明拟合的稳定性好；如果发散或收敛到两种不同结果，说明模型在该数据条件下过参数化。这个测试用一个小循环就能完成，也是评审老师比较关心的部分。
5. 毕业设计答辩前必看的边界条件与排错指南
5.1 数据量、采样率对拟合稳定性的影响
如果数据点数太少，比如不足50个，那么Voigt函数的尾部信息不足以约束gamma和eta，拟合结果会非常不稳定。尤其是gamma参数，它对峰尾的拖尾形状非常敏感，采样范围不够时，gamma会倾向于变大或变小以吸收噪声。一般建议采样范围至少覆盖峰中心两侧各3个半高宽，数据点不少于200。如果原始数据不满足这个条件，可以通过插值增密，但插值会引入额外误差，更好的做法是采集时就保证范围足够。
5.2 初始化不当导致的局部极小值与规避
lsqcurvefit是局部优化算法，初值离全局最优太远时容易陷入局部极小。规避方法有三个：第一，先用smooth函数做滑动平均，找到平滑数据中的峰值和半高宽，用平滑结果计算初值；第二，先固定eta为0和1，分别用两参数模型（sigma+amp和gamma+amp）拟合，取残差较小的初始eta；第三，使用全局优化工具箱中的MultiStart或GlobalSearch配合lsqcurvefit，虽然计算量更大，但能显著提高找到全局解的概率。对于毕业设计而言，使用MultiStart可以作为一个展示点，代码量也不大。
5.3 工具箱版本差异与setparam注意事项
不同MATLAB版本中优化选项的设置方式有差异。R2013a之前的版本使用optimset，之后推荐optimoptions。如果你打开老代码，看到optimset('Display','final')，在R2024b上运行通常没问题，但如果是optimset('Algorithm','trust-region-reflective')，则有可能因算法选择受限而报错。另一个容易踩的坑是变量名冲突：如果工作区中存在名为param或model的变量，而函数内部也使用相同变量名，会导致MATLAB提示“Subscript indices must either be real positive integers or logicals”之类的错误。强烈建议运行demo之前先执行clear; clc;清理工作区。
5.4 批量处理多个峰时的策略
如果你的光谱数据有多个峰，不能直接使用现在的voigtfit，因为模型函数只支持单个Voigt。常见做法是用findpeaks函数找到所有局部极大值，然后对每个峰截取一个区间单独拟合，最后合并结果。区间范围可以用相邻峰谷作为分界。这样做的问题在于，重叠峰的尾部会相互干扰，严格的处理方式是把所有峰放在一个模型函数中同时拟合，此时参数数量成倍增加，需要给每个峰分配5个参数（sigma、gamma、eta、amplitude、center）。批量拟合时，不同峰的参数初值可以复用同一个先验值，但边界可能不同，需要动态生成上下界向量。
6. 让这个项目往后延伸的3个实用技巧
6.1 用fittype自定义模型，让fit函数更灵活
如果你的MATLAB安装了Curve Fitting Toolbox，可以用fittype把myvoigt包装成可识别的模型，直接使用fit函数。这个方法适合快速批处理，因为fit会自动计算参数的95%置信区间：
ft = fittype('amp * myvoigt(x - center, sigma, gamma, eta)', ... 'independent', 'x', ... 'problem', 'center'); % 如果不希望center自由拟合，设为problem opt = fitoptions(ft); opt.StartPoint = [1, 0.5, 0.3, 3]; [fresult, gof] = fit(x', y_noisy', ft, opt);这里gof结构体包含sse、rsquare、adjrsquare等统计量，比手算更方便。注意fit要求自变量单调递增，如果数据顺序混乱，先用sort排序。
6.2 把参数结果导出成table，直接用于论文
论文中通常需要表格列出多组拟合参数。把每次拟合的结果存成行向量，最后合并成MATLAB table并使用writetable输出csv或xlsx：
T = table(p(1), p(2), p(3), p(4), SSE, R2, ... 'VariableNames', {'sigma', 'gamma', 'eta', 'amplitude', 'SSE', 'R2'}); writetable(T, 'fit_summary.csv');这样生成的csv文件可以直接导入到Excel或LaTeX表格中，避免手动复制粘贴出错。对于有多个样本的数据，可以在for循环中不断用T_temp拼接。
6.3 扩展为多峰Voigt拟合的思路
把单峰Voigt扩展成多峰的方法很简单：模型函数内部循环累加。下面给出一个支持可变峰数量的模型函数框架：
function y = multi_voigt(x, p, npeaks) y = zeros(size(x)); for k = 1:npeaks idx = (k-1)*4 + 1; sigma = p(idx); gamma = p(idx+1); eta = p(idx+2); amp = p(idx+3); center = p(idx+4); y = y + amp * myvoigt(x - center, sigma, gamma, eta); end end调用时，令初始参数向量长度为npeaks*5，依次排列所有峰的参数。需要注意，峰位置的初值必须通过findpeaks预先确定，否则优化算法很容易让峰位置乱跑。多峰拟合对初值要求更高，建议用MultiStart辅助寻找全局解。如果峰之间重叠严重，参数相关性会增强，拟合结果的可解释性下降，此时可考虑将部分参数固定，比如所有峰的eta使用同一个值，减少自由参数数量。这样一来，这个毕业设计已经从简单的单峰拟合延伸到了真正的分析工具，在答辩时可以明显提升项目复杂度评价。
本文还有配套的精品资源，点击获取
