哲●科●文Jerkwin 类别 标签 链接 关于 订阅

使用gnuplot进行光谱展宽和拟合

类别: 科   标签: gnuplot   阅读次数:   版权: (CC) BY-NC-SA ✕

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

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

洛伦兹线型

\[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函数的图像, 并与洛伦兹函数, 高斯函数对比一下, 获得些直观印象.

vp-1.png
1
2
3
4
5
6
7
8
9
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)"

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

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

vp-2.png
1
2
3
4
5
6
7
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)))

光谱展宽

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

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

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

vp-3.png
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
$ben <<EOD
 411.52    0.0139
 412.28    0.0000
 623.36    0.0000
 623.88    0.0000
 687.17  129.9979
 727.19    0.0000
 869.06    0.0000
 874.66    0.0000
 996.19    0.0196
1001.05    0.0000
1023.49    0.0009
1025.07    0.0000
1025.25    0.0000
1068.62    7.2319
1071.69    7.0255
1179.01    0.0018
1203.21    0.0000
1210.16    0.0000
1340.03    0.0000
1382.68    0.0000
1519.97    7.8011
1524.40    7.9705
1666.82    0.0000
1667.35    0.0000
3185.61    0.0034
3195.88    0.0000
3196.31    0.0000
3212.36   34.2590
3212.85   34.9293
3223.34    0.0000
EOD

$pro <<EOD
   5.5703    0.2248
   8.202     0.0517
  11.1187    0.1945
  12.4887    0.0312
  14.7264    0.2752
  20.8601    0.9193
  26.9303    0.6772
  30.9304    1.9412
  34.0546    1.6479
  34.6674    1.2732
  39.6549    1.9459
  42.5091    5.699
  47.2868    2.9254
  49.6395    6.3822
  50.7001    2.3203
  56.0899    3.7654
  63.9532    6.0518
  70.2859    0.1166
  75.8429    2.1368
  77.092     1.9241
  78.465     5.9527
  80.969     3.1261
  86.7711    2.9239
  94.2281    0.0847
  96.7878   10.0573
 102.5753    2.2876
 103.5022   13.5716
 109.5904    4.1824
 117.5981    1.4245
 120.7254    9.0271
 127.1606    2.6787
 131.1905   11.464
 146.7241   17.3003
 148.9796    5.2004
 151.8352    4.1021
 164.6458    2.014
 168.4281    1.5959
 178.56      1.2548
 179.8289    1.3345
 192.0131   18.7003
 207.8927    5.8878
 211.2551    4.0618
 214.0112    5.7298
 221.5261   27.1325
 225.232     7.3722
 242.6048    9.2672
 266.769     9.7515
 277.6434   13.3218
 284.3466   12.0448
 296.2707   27.5914
 313.9051    2.5382
 326.4214    0.8269
 331.7138    8.9124
 334.2993    3.4369
 340.2758   15.4909
 360.3065    0.3587
 372.9981    1.7383
 375.1391    2.7505
 395.5045   10.3585
 396.1989    7.1265
 400.5555   28.8306
 409.9106   34.5589
 432.7897    5.0126
 455.7656   72.503
 463.0265   95.8982
 469.9515   22.4316
 485.0086   61.404
 486.3107  144.6949
 491.2587   41.7086
 565.8919    9.448
 567.6956   27.495
 601.6105   40.0251
 616.9299   17.9848
 619.2627  105.0766
 630.8975   17.8575
 658.5038   69.1035
 664.4005    4.2228
 669.5269   70.6619
 681.0125   92.2443
 696.3436   11.7463
 700.3316   56.5851
 701.3512   28.4417
 710.7345   35.8959
 714.5113    6.6981
 716.9605   30.2688
 720.5791   47.7707
 732.277    31.6282
 741.0209   56.6437
 747.9012   12.4558
 750.3384    4.3337
 767.4022   18.4507
 782.7754    3.1596
 790.3709   33.9786
 808.49      3.1682
 836.5271   23.3033
 863.5872    7.3729
 864.7526   10.9439
 881.4292   11.6686
 905.2745   81.3366
 909.3038   13.0176
 936.6016   35.3285
 938.4601   64.7013
 944.148    97.6404
 961.2511   34.9923
 961.47     29.2459
1000.0235    3.5754
1006.7099    1.7483
1016.3939    7.4102
1023.1448    1.9533
1025.4377    2.7167
1064.4731   13.3154
1069.9063   12.2305
1071.8959   22.7919
1073.7106   12.824
1077.035     3.0798
1087.8828   24.381
1088.7724    0.5369
1091.3019    7.156
1092.6817    0.9176
1096.0523    3.5179
1106.5754   12.0345
1110.7173   15.6929
1116.1719   25.9462
1127.6926    5.684
1128.2485   10.4408
1135.2589    9.0076
1135.3634   19.894
1170.9492   12.6672
1171.5141    7.1377
1176.381    14.6126
1185.3512   14.374
1192.013    10.5838
1203.5899  140.7856
1204.1205  168.3034
1219.4349   16.6149
1243.2615   72.2718
1244.3306   19.7603
1249.9352   16.8338
1253.7756    4.2009
1257.1176   12.43
1258.2583   10.9626
1274.3982    8.9235
1277.0126   73.3254
1294.9685   25.9155
1297.3158    1.4366
1297.5511   10.0741
1301.119    18.5913
1313.3604   14.4705
1314.5702   25.6698
1327.0457   11.3781
1348.8651   65.0611
1360.6187   26.1801
1365.4324   18.2422
1365.9126   14.7917
1371.9704   21.2766
1372.2487  181.1006
1398.991    84.3235
1408.7975   48.6287
1415.7099  129.5311
1428.2345    6.1065
1429.7462   51.4848
1443.1475   93.2618
1452.4158   84.7856
1459.1361   12.9864
1459.7117   10.1549
1468.7118   58.5782
1489.0287   37.4598
1490.1846   28.3128
1502.2183   24.6394
1518.0248   68.7913
1528.3239   54.046
1529.8767   12.7381
1546.5526    4.6335
1558.1716   32.5112
1558.2001   35.6352
1562.1631   39.1299
1566.4315   88.1398
1592.0727   64.9987
1603.7566   29.2215
1605.9611   12.4568
1606.14     16.5949
1608.7598   17.4199
1608.9714   45.8923
1611.0471   49.0506
1612.5791   29.7291
1622.6447   34.6206
1625.8612   12.5778
1628.0334    1.641
1629.0302  159.3846
1629.2513    3.1434
1630.6931    1.4976
1631.9014   25.1015
1633.3015    6.605
1633.9271   11.1993
1641.0299    9.9203
1641.7938   29.8011
1649.826    19.5962
1650.9025   12.7789
1676.6164  135.773
1690.6269  367.0851
1697.1001   72.4277
1697.3751  518.5653
1700.597   311.6447
1703.253   234.2651
1708.5578    1.8341
1726.4863  264.9281
1731.3218  153.7944
1802.943    78.9193
1811.1822   33.7786
1821.9025   95.8968
1897.9733  507.8488
1908.7425  740.5449
1912.1115   76.4197
1948.6744  298.6536
1952.2258  190.9105
1964.7986  232.6757
3222.6301    3.2267
3227.4149    1.4081
3233.4723    2.4421
3237.1163   41.8467
3237.5387    9.6233
3242.6788   38.4305
3245.535    31.3443
3251.8003    5.6615
3255.9425    3.5297
3260.8563   16.8681
3276.8333    1.83
3283.2845    7.7968
3288.5065   10.6577
3294.0036   12.0576
3296.7166    7.5492
3308.7381    1.9153
3309.4265    2.9439
3312.3182   25.3147
3317.6555   21.2113
3317.895    20.889
3318.5208   12.974
3328.6494   11.0121
3334.3257    0.5727
3338.0091   19.8689
3340.0895    3.3902
3346.1489    9.9576
3349.8259    8.1264
3447.1727  128.4542
3456.8082    1.9686
3458.2356    5.0307
3465.5113    0.466
3468.7609    4.1457
3470.5013    2.5949
3718.79    592.2654
3737.2726  504.0469
3822.376   132.18
3841.1183   73.0049
3856.8681   52.7314
3863.233    65.2451
3871.7072  157.5003
3884.8637   42.4131
3887.9648   52.3414
EOD

	data="$ben" # 选择绘图数据

	stat data u 1 nooutput       # 统计峰的数目
	Npnt=STATS_records
	array c[Npnt]; array h[Npnt] # 定义数组保存峰位置, 高度

	stat data u (c[$0+1]=$1, 1):(h[$0+1]=$2, 2) nooutput # 填充数组

	# 定义展宽函数
	peak(x,c,h,γG,γL)=h*VP(x-c,γG/sqrt(log(4)),γL)

	# 展宽参数
	γG=100
	γL=50
	# γV=VP_fwhm(γG/sqrt(log(4)),γL) # gnuplot 5.4尚不支持此函数

	set xl "cm^{-1}"
	set yl "a.u."
	set tit sprintf("γ_G=%.1f    γ_L=%.1f", γG, γL)

	# 绘制, 使用累加函数
	set boxwidth 0.01    # 使用柱状图画竖线, 5.5才支持Impulses
	p [0:4000] \
		data axis x1y2 w boxes lw .5 t "", \
		sum [i=1:Npnt] peak(x,c[i],h[i],γG,0) t "Gauss", \
		sum [i=1:Npnt] peak(x,c[i],h[i],0,γL) t "Lorentz", \
		sum [i=1:Npnt] peak(x,c[i],h[i],γG,γL) t "Voigt"

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

光谱拟合

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

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

vp-5.png
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
$dat <<EOD
  0   3.3487290833206163
  1   3.441076831745743
  2   7.7932863251851305
  3   7.519064207516034
  4   7.394406511652473
  5  11.251458210206666
  6   4.679476113847004
  7   8.313048016542345
  8   9.348006472917458
  9   6.086336477997078
 10  10.765370342398741
 11  11.402519337778239
 12  11.151689287913552
 13   8.546151698722557
 14   8.323886291540909
 15   7.133249200994414
 16  10.242189407441712
 17   8.887686444395982
 18  10.759444780127321
 19   9.21095463298772
 20  15.693160143294264
 21   9.239683298899614
 22   9.476116297451632
 23  10.128625585058783
 24  10.94392508956097
 25  10.274287987647595
 26   9.552394167463973
 27   9.51931115335406
 28   9.923989117054466
 29   8.646255122559495
 30  12.207746464070603
 31  15.249531807666745
 32   9.820667193850705
 33  11.913964012172858
 34   9.506862412612637
 35  15.858588835799232
 36  14.918486963658015
 37  15.089436171053094
 38  14.38496801289269
 39  14.42394419048644
 40  15.759311758218061
 41  17.063349232010786
 42  12.232863723786215
 43  10.988245956134314
 44  19.109899560493286
 45  18.344353100589824
 46  17.397232553539542
 47  12.372706600456558
 48  13.038720878764792
 49  19.100965014037367
 50  17.094480819566147
 51  20.801679461435484
 52  15.763762333448557
 53  22.302320507719728
 54  23.394129891315963
 55  19.884812694503303
 56  22.09743700979689
 57  16.995815335935077
 58  24.286037929073284
 59  25.214705826961016
 60  25.305223543285013
 61  22.656121668613896
 62  30.185701748800568
 63  28.28382587095781
 64  35.63753811848088
 65  35.59816270398698
 66  35.64529822281625
 67  36.213428394807224
 68  39.56541841125095
 69  46.360702383473075
 70  55.84449512752349
 71  64.50142387788203
 72  77.75090937376423
 73  83.00423387164669
 74 111.98365374689226
 75 121.05211901294848
 76 176.82062069814936
 77 198.46769832454626
 78 210.52624393366017
 79 215.36708238568033
 80 221.58003148955638
 81 209.7551225151964
 82 198.4104196333782
 83 168.13949002992925
 84 126.0081896958841
 85 110.39003569380478
 86  90.88743461485616
 87  60.5443025644061
 88  71.00628698937221
 89  61.616294708485384
 90  45.32803695045095
 91  43.85638472551629
 92  48.863070901568086
 93  44.65252243455522
 94  41.209120125948104
 95  36.63478075990383
 96  36.098369542551325
 97  37.75419965137265
 98  41.102019290969956
 99  26.874409332756752
100  24.63314900554918
101  26.05340465966265
102  26.787053802870535
103  16.51559065528567
104  19.367731289491633
105  17.794958746427422
106  19.52785218727518
107  15.437635249660396
108  21.96712662378481
109  15.311043443598177
110  16.49893493905559
111  16.41202114648668
112  17.904512123179114
113  14.198812322372405
114  15.296623848360126
115  14.39383356078112
116  10.807540004905345
117  17.405310725810278
118  15.309786310492559
119  15.117665282794073
120  15.926377010540376
121  14.000223621497955
122  15.827757539949431
123  19.22355433703294
124  12.278007446886507
125  14.822245428954957
126  13.226674931853903
127  10.551237809932955
128   8.58081654372226
129  10.329123069771072
130  13.709943935412294
131  11.778442391614956
132  14.454930746849122
133  10.023352452542506
134  11.01463585064886
135  10.621062477382623
136   9.29665510291416
137   9.633579419680572
138  11.482703531988037
139   9.819073927883121
140  12.095918617534196
141   9.820590920621864
142   9.620109753045565
143  13.215701804432598
144   8.092085538619543
145   9.828015669152578
146   8.259655585415379
147   9.424189583067022
148  13.149985946123934
149   7.471175119197948
150  10.947567075630904
151  10.777888096711512
152   8.477442195191612
153   9.585429992609711
154   7.032549866566089
155   5.103962051624133
156   9.285999577275545
157   7.421574444036404
158   5.740841317806245
159   2.3672530845679
EOD

Ag=100; cg=80; sg=6;
Al=100; cl=80; gl=6;
A=100;  c=80;  s=6; g=6;
fit Ag*VP(x-cg, sg, 0) $dat u 1:2 via Ag, cg, sg
fit Al*VP(x-cl, 0, gl) $dat u 1:2 via Al, cl, gl
fit A*VP(x-c, s, g)    $dat u 1:2 via A, c, s, g

set xl"x"
set yl"y"
p $dat w p t"data", \
	Ag*VP(x-cg, sg, 0) t "Gauss", \
	Al*VP(x-cl, 0, gl) t "Lorentz", \
	A *VP(x-c,  s, g)  t "Voigt"

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

思考

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆
前一篇:　曾经的海
后一篇:　VMD的图形矩阵及其应用
访问人次(2015年7月 9日起): | 最后更新: 2024-04-16 06:38:20 UTC | 版权所有 © 2008 - 2024 Jerkwin 个