仓库源文，站点原文

layout: post title: 使用gnuplot进行光谱展宽和拟合 categories:

科 tags:
gnuplot math: true

2022-07-27 22:25:35

在光谱数据的处理方面, 有两类常见任务:

光谱展宽: 将理论计算得到的离散线状谱展宽成类似实验数据的连续谱, 便于比较.
光谱拟合: 将实验得到的连续谱先分峰, 再拟合每个峰的数据.

这两类任务都要涉及光谱中单个峰的数学表达式, 一般称其为线型函数. 基本的线型函数有两类;

洛伦兹线型

$$g_L(ν,γ_L)={γ_L／π}{1／ν^2+γ_L^2}$$

半峰半宽 $γ_L$, 半峰全宽 $2γ_L$

高斯线型, 也称多普勒线型

$$g_G(ν,σ)={1／√{2π}σ}\exp（-{ν^2／2σ^2}）, σ={γ_G／√{2㏑2} }$$

半峰半宽 $γ_G=√{2㏑2}σ$, 半峰全宽 $2γ_G$

线型函数实际上对应于数学上所言的概率密度分布函数, 因此, 也可以将其称为分布. 洛伦兹分布物理上常用, 数学上则称为柯西分布. 高斯分布也就是正态分布. 这两种分布都是t分布的特例, 前者为自由度为1的t分布, 后者为自由度趋向无穷大的t分布.

洛伦兹线型来源于光谱的自然增宽和压力增宽, 表明体系内部具有指数衰减过程. 高斯线型来自多普勒效应, 这也是常称其为多普勒线型的原因. 这种增宽的原因在于温度非零时分子的速度具有麦克斯韦分布(类似高斯分布). 实际条件下测得的光谱是这两种线型函数共同作用的结果, 可写为卷积的形式, 所对应的线型称为Voigt函数. 这个函数为一个积分, 不存在简单的表达式, 计算比较麻烦, 所以除非必要, 一般不会用它进行光谱展宽. 但有时单独用洛伦兹展宽或高斯展宽效果不好, 为此, 又提出了一种赝Voigt函数来近似Voigt函数, 它是洛伦兹和高斯的线性组合, 计算简单, 所以常被用来作为效果稍好些的线型函数.

因为Voigt函数在光谱领域非常重要, 所以有大量论文讨论如何快速精确地计算其值. 各种科学计算程序也大多支持Voigt函数的计算. 目前常用的是libcert函数库, 这也是gnuplot计算Voigt函数所用的方法. Python中也有相应的Voigt实现.

我以前写过简单的高斯, 洛伦兹, 赝Voigt展宽脚本, 但一直没有实现Voigt展宽, 因为自己写代码计算Voigt函数并没那么简单. 最近我发现Windows下的gnuplot 5.4自带了Voigt函数, 所以觉得没有必要再编译libcert或编写Voigt函数了, 可以借助gnuplot实现一个比较通用的光谱展宽方案.

函数测试

先看看gnuplot中VP函数的图像, 并与洛伦兹函数, 高斯函数对比一下, 获得些直观印象.

<table class="highlighttable"><th colspan="2" style="text-align:left">vp-1.png</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height:125%">1 2 3 4 5 6 7 8 9</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%">Gauss(x,σ)=exp(-(x/σ)**2/2)/(sqrt(2*pi)*σ) Lorentz(x,γ)=γ/(pi*(x*x+γ*γ)) VP01(x)=(abs(x-int(x))<20/1E4?VP(x,0,1):1/0) VP10(x)=(abs(x-int(x))<20/1E4?VP(x,1,0):1/0) set xl "x" set yl "g(x)" set tit "Gauss(x,σ)=exp(-½x^2/σ^2)/√(2π)σ=VP(x,σ,0)\nLorentz(x,γ)=γ/(π(x^2+γ^2))=VP(x,0,γ)" p Gauss(x, 1), Lorentz(x,1), VP(x,1,1), VP10(x) w p t"VP(x,1,0)", VP01(x) w p t"VP(x,0,1)"</pre></div> </td></tr></table>

可以看到, 在归一化的情况下, 高斯函数高瘦, 衰减更快, 洛伦兹函数则具有肥尾特征. VP函数在σ=0, γ=0时分别趋近于洛伦兹函数和高斯函数.

在半峰宽相同的条件下, VP函数的表现则有不同, 这是有些误解产生的原因, 具体讨论的可以参考论文A Common Misunderstanding about the Voigt Line Profile以及Comments on “A Common Misunderstanding about the Voigt Line Profile”. 将下面的图与上面的进行对比, 就可以发现其中的区别.

<table class="highlighttable"><th colspan="2" style="text-align:left">vp-2.png</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height:125%">1 2 3 4 5 6 7</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%">Gauss(x,σ)=exp(-(x/σ)**2/2)/(sqrt(2*pi)*σ) Lorentz(x,γ)=γ/(pi*(x*x+γ*γ)) set xl "x" set yl "g(x)" set tit "Profiles with Unit Half Widths" p Gauss(x, 1/sqrt(2*log(2))), Lorentz(x, 1), VP(x, 2/(1+sqrt(5)), 2/(1+sqrt(5)))</pre></div> </td></tr></table>

光谱展宽

光谱展宽的数学比较简单. 假定我们有离散的一系列峰位置和峰高, $c_i, h_i, i=1,⋯,n$, 在每个峰的中心放置一个线型函数$h_i f(ν)$, 然后将所有线型函数的值累加起来作为连续光谱的值:

$$g(ν)=∑_{i=1}^n h_i f(ν-c_i)$$

下面是两个示例. 第一个是苯的红外光谱, 因为分子的对称性较高, 非零强度的频率较少.

另一个是某蛋白部分的红外光谱, 没有对称性, 所以有很多峰.

光谱拟合

有了光谱展宽的经验, 处理光谱拟合也就没有太大问题了. 但光谱拟合的难点其实不在如何拟合, 而是在于分峰. 实验获得的光谱通常比较复杂, 在没有其他信息的条件下, 如何确定分峰的数目是一件很玄学的事情. 虽然也有一些数学上的处理方法, 但也不能作为最后的结论, 还是要结合实际经验和理论推测才成. 网上有很多拟合的示例, 但大都是基于赝Voig函数的:

下面示意一下如何使用gnuplot进行单峰拟合. 如果确定了分峰数目, 多峰拟合的过程大致类似.

可以看到, 洛伦兹函数和Voigt函数的拟合结果区别不大, 这说明测试数据更接近洛伦兹函数.

思考

实际每个峰对应的半峰宽并不一定是常数, 试着改进代码以支持变化的半峰宽
对于不同类型的光谱, 理论计算所给出的强度单位不同, 展宽后峰的高度与实际单位的对应情况也不同, 参考使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图, 绘制具有实际单位的谱图
查找一些实验谱图, 将其数字化, 然后试着拟合, 看不同类型的实验谱图倾向于哪种线型函数
如果对Voigt函数背后的数学感兴趣, 可以看看Voigt线型函数及其最大值的研究, 以及Simple Analytical Expression of the Voigt Proﬁle