目录
https://github.com/Gleland/SpectralAnalysis引言
继续分享一个有关光谱分析的Python代码。
朗伯-比尔定律
朗伯-比尔定律，这是一个描述溶质在溶液中吸收光线强度与溶液浓度和光程长度成正比关系的基本定律。朗伯-比尔定律的数学表达式为：
其中A为吸光度；k为吸收系数；c为溶液浓度；L为紫外线光程；I0为入射光强度；I为透射光强度。
代码目的
该程序将通过透射光谱的比尔-朗伯定律计算吸收系数。
运行代码步骤
该分析项目包含了两个CSV文件作为示例数据。RawBackground.csv和RawData.csv是在92开尔文下收集的空腔室和纯甲烷 (CH4) 液体的光谱数据。运行代码步骤如下所示：
- 需要通过终端运行代码
- 导入用于分析的数据。程序首先会要求输入原始数据集，然后是背景光谱。这些数据文件不需要与位于同一目录中，但如果是这种情况，请确保提供完整路径名。链接中包含了下载中可以用于测试该程序的数据。
- 导入这些数据文件后，程序将保存您正在分析的混合物的温度和浓度。程序将询问这两个信息，并将每个答案分别存储在目录中的一个单独的txt文件中。
- 数据导入后，程序将绘制两个数据集供用户检查（以确保数据正常），并询问用户希望使用哪种滤波方法。选项包括使用 Savitzky-Golay滤波器（本质上是一种低通滤波器）或快速傅里叶变换方法（FFT）手动去除不需要的频率。此滤波步骤用于去除可能来自实验方法的任何信号。在此项目中，有一个信号来自实验室设置中蓝宝石窗的反射。有关更多信息，请参见 Protopapa 2015和Grundy2002。建议使用FFT 方法，因为这允许用户对处理的数据有更多控制。
- 用户确定滤波方法后，程序将询问用户选择的进一步输入。如果选择了Savitzky-Golay滤波器，用户需要定义一个窗口框大小和算法的多项式。请参阅维基百科文章以了解其工作原理。如果用户选择了FFT选项，将对数据进行FFT，并且用户手动选择要去除的频率。程序监听在显示的图表上的鼠标点击，左键点击会添加一条垂直线，中键点击会移除最近的线，右键点击表示用户对选择满意并希望继续。在提供的gif中可以看到FFT是对称的。程序将选定的频率（在gif中约为0.7cm）及其在原点另一侧的对应频率设为零。
- 滤波过程完成后，下一步是拟合一个连续谱。这个连续谱定义了没有吸收发生的地方。程序将要求用户放大希望拟合连续谱的区域，然后询问要使用的多项式（线性、二次等）。在下方的gif中，使用了三次（立方）拟合，但用户可以选择任何拟合。当用户在图表上点击时，程序会对添加到图表上的所有数据点进行多项式拟合。左键点击添加一个点，中键点击移除最近的点，右键点击结束采集。
- 拟合连续谱后，程序将要求用户输入在光谱采集期间使用的池的厚度。这是因为光的吸收取决于通过比尔-朗伯定律的池腔厚度。我的实验中使用的池腔厚度为2cm，因此将输入该值。最后，程序将生成吸收系数并显示给用户检查。注意，并非所有的光谱都可以依赖于合法数据，终点（小于90000 cm-1和大于15000 cm-1）不被认为是显著的，因此请确保不在计算中包括这些值（2中数据）。
- 查看图表结束后，用户输入，程序关闭。程序运行期间生成的所有数据和图表将保存到开始时输入的目录中。
运行过程展示
上面提到的GIF是做作者给出的程序运行过程填写数据和信息的GIF图片，这里我自己运行了一遍。
结果图
这里附上光谱分析的代码：
import os
import csv # 用于处理CSV文件的读写
import time # 用于处理时间相关的功能
import warnings # 用于处理警告信息
import numpy as np # 用于科学计算
import matplotlib.pyplot as plt # 用于绘图
from scipy.optimize import curve_fit as cf # 用于曲线拟合
from scipy.fftpack import fft, fftfreq, ifft # 用于快速傅里叶变换
from scipy.signal import savgol_filter as sgf # 用于Savitzky-Golay滤波
from scipy.integrate import trapz # 用于数值积分
def main():
# 选择保存数据的文件夹
folder_to_save = choose_dir()
# 选择要分析的文件
raw_x, raw_y, raw_xbg, raw_ybg = choose_files(folder_to_save)
print("绘制导入的数据...")
# 绘制原始数据图
plotting_data_for_inspection(raw_x, raw_y, 'Raw Data', 'Wavenumber (cm-1)', '% Transmittance', 'rawspectrum.pdf', folder_to_save, False)
# 绘制原始背景数据图
plotting_data_for_inspection(raw_xbg, raw_ybg, 'Raw Background', 'Wavenumber (cm-1)', '% Transmittance', 'rawbackground.pdf', folder_to_save, False)
# 用户在检查图后选择方法
user_method = str(input('按 "s" 使用Savitzky-Golay滤波，按 "f" 使用FFT滤波\n:'))
choosing = True
while choosing:
if user_method.lower() == 's':
# 选择Savitzky-Golay滤波方法
choosing = False
args_list = [folder_to_save, raw_y, raw_ybg, raw_x]
raw_x, norm_smooth = sgf_calc(args_list)
plot_data(raw_x, norm_smooth, folder_to_save)
elif user_method.lower() == 'f':
# 选择FFT滤波方法
choosing = False
frq_x, frq_xbg, fft_y, fft_ybg = fft_calculation(raw_x, raw_y, raw_xbg, raw_ybg, folder_to_save)
plot_figure, plot_axis = plotting_data_for_inspection(frq_x, np.log(abs(fft_ybg)), 'FFT of raw bg', 'Cycles/Wavenumber (cm)', 'Log(Power/Frequency)', 'fft_background.pdf', folder_to_save, False)
filt_y = fft_y.copy()
filt_ybg = fft_ybg.copy()
input('放大到合适的范围，然后按回车键开始')
print('左键添加，中央键移除最近的，右键完成')
# 连接鼠标点击事件
vert_lines = []
frq_cid = plot_figure.canvas.mpl_connect('button_press_event', lambda event: freq_click(event, [frq_x, fft_ybg, plot_figure, plot_axis, vert_lines, filt_y, filt_ybg, folder_to_save, raw_x]))
plt.show()
plot_figure.canvas.mpl_disconnect(frq_cid)
# 断开鼠标点击事件
def save_as_csv(folder_to_save, title, column1_title, column2_title, column1_data, column2_data):
# 切换到保存数据的目录
os.chdir(folder_to_save)
with open(title, "w") as f:
writer = csv.writer(f)
writer.writerow([column1_title, column2_title])
writer.writerows(list(zip(column1_data, column2_data)))
# 切换回上一级目录
os.chdir('..')
def fft_calculation(raw_x, raw_y, raw_xbg, raw_ybg, folder_to_save):
"""
计算数据的FFT以去除不需要的频率
"""
# 计算y数据的FFT
fft_y = fft(raw_y)
fft_ybg = fft(raw_ybg)
# 获取数据FFT的频率和采样间隔
frq_x = fftfreq(len(fft_y), ((max(raw_x) - min(raw_x)) / len(fft_y)))
frq_xbg = fftfreq(len(fft_ybg), ((max(raw_xbg) - min(raw_xbg)) / len(fft_ybg)))
save_as_csv(folder_to_save, "FFT_Raw_bg_data.csv", "frq_x", "log(abs(fft_bg))", frq_x, np.log(abs(fft_ybg)))
return frq_x, frq_xbg, fft_y, fft_ybg
def choose_dir():
"""
用户选择保存所有工作的目录，并生成时间戳以供日后参考
"""
# 所有后续工作保存的目录
folder_to_save = input('输入要保存所有数据的目录名称\n:')
# 创建并切换到用户命名的目录
os.mkdir(folder_to_save)
os.chdir(folder_to_save)
# 记录程序运行的日期和时间，并将其保存到文件中
with open("time_created.txt", "w") as text_file:
text_file.write("程序运行时间: {} \n".format(time.strftime("%Y-%m-%d %H:%M")))
# 切换回上一级目录
os.chdir('..')
return folder_to_save
def plotting_data_for_inspection(xdata, ydata, plot_title, plot_xlabel, plot_ylabel, filename_for_saving, folder_to_save, block_boolean):
"""
在程序内绘制数据供用户查看
参数
----------
xdata, ydata: 要绘制的x和y数据
plot_xlabel, plot_ylabel: 图表x轴和y轴的标签
file_name_for_saving: 保存文件时的文件名字符串
block_boolean: 布尔值，指示程序是否等待图表关闭
"""
plot_figure, plot_axis = plt.subplots()
plt.plot(xdata, ydata, color='blue')
plt.xlabel(plot_xlabel)
plt.ylabel(plot_ylabel)
plt.suptitle(plot_title)
plt.show(block=block_boolean)
os.chdir(folder_to_save)
plt.savefig(filename_for_saving)
os.chdir('..')
return plot_figure, plot_axis
def choose_files(folder_to_save):
"""
让用户确定要导入进行分析的文件，并保存这些偏好以供以后参考
"""
raw_import = str(input('输入要分析的原始数据集\n:'))
print("\n好的！现在导入...\n")
raw_x, raw_y = import_data(raw_import)
bg_import = str(input('输入要分析的原始背景数据\n:'))
print("\n好的！现在导入...\n")
raw_xbg, raw_ybg = import_data(bg_import)
os.chdir(folder_to_save)
with open("data_files_used.txt", "w") as text_file:
text_file.write("使用的原始数据文件: {} \n".format(raw_import))
text_file.write("使用的原始背景数据文件: {}".format(bg_import))
concentration = str(input('输入混合物的浓度\n:'))
# 将浓度保存到文本文件以供日后使用
with open("concentration.txt", "w") as f:
f.write(concentration)
temperature = str(input('输入混合物的温度\n:'))
# 将温度保存到文本文件以供日后使用
with open("temperature.txt", "w") as f:
f.write(temperature)
os.chdir('..')
return raw_x, raw_y, raw_xbg, raw_ybg
# 假设数据存储在CSV文件中，因为冰实验室的所有数据都以CSV格式存储
def import_data(filename):
# 从CSV文件中加载数据
raw_data = np.loadtxt(open(filename, "rb"), delimiter=",")
# 提取第一列数据作为x轴数据
xdat = raw_data[:, 0]
# 提取第二列数据作为y轴数据
ydat = raw_data[:, 1]
# 返回x轴和y轴数据
return xdat, ydat
def freq_click(event, args_list):
# 如果button_click是左键：添加左线
# 如果button_click是中键：移除最近的线
# 如果button_click是右键：完成
# 将点击的数据点添加到列表中
frq_x, fft_ybg, plot_figure, plot_axis, vert_lines, filt_y, filt_ybg, folder_to_save, raw_x = args_list
# 设置x轴和y轴的范围
plt.xlim(plt.gca().get_xlim())
plt.ylim(plt.gca().get_ylim())
if event.button == 1:
# 左键点击，添加垂直线
vert_lines.append(event.xdata)
plot_axis.plot(frq_x, np.log(np.abs(fft_ybg)), color='blue')
for val in vert_lines:
plt.axvline(x=val, color='black')
plt.xlabel('Cycles/Wavenumber')
plt.ylabel('Relative Intensity')
plt.draw()
if event.button == 2:
# 中键点击，移除最近的垂直线
print('pop!')
xlims = plt.gca().get_xlim()
ylims = plt.gca().get_ylim()
plot_axis.cla()
plot_axis.plot(frq_x, np.log(np.abs(fft_ybg)), color='blue')
plt.xlim(xlims)
plt.ylim(ylims)
plt.xlabel('Cycles/Wavenumber')
plt.ylabel('Relative Intensity')
xindx = np.abs(vert_lines - event.xdata).argmin()
del vert_lines[xindx]
for line in vert_lines:
plt.axvline(x=line, color='black')
plt.draw()
if event.button == 3:
# 右键点击，结束点击检测
os.chdir(folder_to_save)
plt.savefig('FFT_filter.pdf')
with open("freq_window.csv", "w") as f:
writer = csv.writer(f)
writer.writerow(["Xposition of vert. line"])
writer.writerows(list(zip(vert_lines)))
os.chdir('..')
args_dict = {"vert_lines": vert_lines, "frq_x": frq_x, "filt_y": filt_y, "filt_ybg": filt_ybg}
plt.close("all")
argslist = [vert_lines, frq_x, filt_y, filt_ybg]
filt_y, filt_ybg = window_filter(argslist)
fft_calc(filt_y, filt_ybg, raw_x, folder_to_save)
def fft_calc(filt_y, filt_ybg, raw_x, folder_to_save):
# 将过滤后的y数据与过滤后的背景数据相除
norm_fft = ifft(filt_y) / ifft(filt_ybg)
save_as_csv(folder_to_save, "fft_data.csv", "raw_x", "fft_filt", raw_x, norm_fft.real)
plot_data(raw_x, norm_fft.real, folder_to_save)
def sgf_calc(args_list):
folder_to_save, raw_y, raw_ybg, raw_x = args_list
# 使用sgf选项时的警告
warnings.filterwarnings(action="ignore", module="scipy", message="^internal gelsd")
window_param = int(input('输入窗口大小（必须是奇数）\n:'))
poly_param = int(input('输入平滑多项式的阶数\n:'))
# 保存选择的参数以供将来检查
os.chdir(folder_to_save)
with open("sgf_params.txt", "w") as sgf_file:
sgf_file.write("使用的窗口参数: {} \n".format(window_param))
sgf_file.write("使用的多项式参数: {}".format(poly_param))
smoothed_y = sgf(raw_y, window_param, poly_param, delta=(abs(raw_y)[1] - raw_y)[0])
smoothed_ybg = sgf(raw_ybg, window_param, poly_param, delta=(abs(raw_ybg)[1] - raw_ybg)[0])
norm_smooth = smoothed_y / smoothed_ybg
rows = list(zip(raw_x, norm_smooth))
with open("sgf_data.csv", "w") as f:
writer = csv.writer(f)
writer.writerow(["window", "polynomail order"])
writer.writerow([window_param, poly_param])
writer.writerow(["raw_x", "sgf_filt"])
writer.writerows(rows)
os.chdir('..')
return raw_x, norm_smooth
# 过滤频率范围
def window_filter(args_list):
vert_lines, frq_x, filt_y, filt_ybg = args_list
window_min, window_max = vert_lines[-2], vert_lines[-1]
for i in range(len(frq_x)):
if (frq_x[i] >= window_min and frq_x[i] <= window_max) or (frq_x[i] > -1 * window_max and frq_x[i] < -1 * window_min):
filt_y[i] = 0
filt_ybg[i] = 0
return filt_y, filt_ybg
def plot_data(x, y, folder_to_save):
plot_figure, plot_axis = plotting_data_for_inspection(x, y, "Divide and Filtered Spectrum", "Wavenumber cm-1", "Relative Intensity", "dv_filt_spectrum.pdf", folder_to_save, False)
order = int(input('调整视图到满意后输入多项式阶数以进行连续拟合\n:'))
xcoords, ycoords = [], []
global cid
cid = plot_figure.canvas.mpl_connect('button_press_event', lambda event: onclick(event, [xcoords, ycoords, plot_figure, plot_axis, order, folder_to_save, x, y]))
print('左键添加，中键删除最近点，右键完成')
plt.show()
def onclick(event, argslist):
xcoords, ycoords, plot_figure, plot_axis, order, folder_to_save, x, y = argslist
global pvals
if event.button == 1:
# 左键点击
plt.xlim(plt.gca().get_xlim())
plt.ylim(plt.gca().get_ylim())
try:
if len(plot_axis.lines) != 1: plot_axis.lines.remove(plot_axis.lines[-1])
except: UnboundLocalError
xcoords.append(event.xdata)
ycoords.append(event.ydata)
plot_axis.scatter(xcoords, ycoords, color='black')
plt.xlabel('Wavenumber cm-1')
plt.ylabel('Relative Intensity')
plt.draw()
xvals = np.array(xcoords)
yvals = np.array(ycoords)
warnings.simplefilter('ignore', np.RankWarning)
p_fit = np.polyfit(xvals, yvals, order)
pvals = np.poly1d(p_fit)
plot_axis.plot(x, pvals(x), color='black')
plt.draw()
if event.button == 2:
# 中键点击，移除最近的点
print('pop!')
xlims = plt.gca().get_xlim()
ylims = plt.gca().get_ylim()
plot_axis.cla()
plot_axis.plot(x, y)
plt.xlim(xlims)
plt.ylim(ylims)
plt.xlabel('Wavenumber cm-1')
plt.ylabel('Relative Intensity')
xindx = np.abs(xcoords - event.xdata).argmin()
del xcoords[xindx]
yindx = np.abs(ycoords - event.ydata).argmin()
del ycoords[yindx]
plot_axis.scatter(xcoords, ycoords, color='black')
plt.draw()
xvals = np.array(xcoords)
yvals = np.array(ycoords)
warnings.simplefilter('ignore', np.RankWarning)
p_fit = np.polyfit(xvals, yvals, order)
pvals = np.poly1d(p_fit)
plot_axis.plot(x, pvals(x), color='black')
plt.draw()
if event.button == 3:
# 右键点击，结束点击检测
plot_figure.canvas.mpl_disconnect(cid)
os.chdir(folder_to_save)
plt.savefig('continuum_chosen.pdf')
with open("continuum_polynomial.txt", "w") as save_file:
save_file.write("%s *x^ %d " % (pvals[0], 0))
for i in (range(len(pvals))):
save_file.write("+ %s *x^ %d " % (pvals[i + 1], i + 1))
os.chdir('..')
calc_coeffs(pvals, x, y, folder_to_save)
def calc_coeffs(pvals, x, y, folder_to_save):
fit_y = pvals(x)
new_continuum = y / fit_y
thickness = int(input('\n输入样品厚度（单位：厘米）\n:'))
err_settings = np.seterr(invalid='ignore')
alpha_coeffs = -np.log(new_continuum) / thickness
plotting_data_for_inspection(x, alpha_coeffs, "Alpha Coefficients", "Wavenumber cm-1", "Absorption cm-1", "alpha_coeffs.pdf", folder_to_save, False)
save_as_csv(folder_to_save, "alpha_coeffs.csv", "x", "alpha", x, alpha_coeffs)
x_mask1 = x[(x > 10000) & (x < 10500)]
x_mask2 = x[(x > 11200) & (x < 12000)]
y_mask1 = alpha_coeffs[(x > 10000) & (x < 10500)]
y_mask2 = alpha_coeffs[(x > 11200) & (x < 12000)]
save_as_csv(folder_to_save, "10000_peak.csv", "x", "y", x_mask1, y_mask1)
save_as_csv(folder_to_save, "11200_peak.csv", "x", "y", x_mask2, y_mask2)
area10000 = trapz(y_mask1, x_mask1)
area11200 = trapz(y_mask2, x_mask2)
os.chdir(folder_to_save)
with open("10000area.txt", "w") as f:
f.write(str(area10000))
with open("11200area.txt", "w") as f:
f.write(str(area11200))
os.chdir('..')
finish_prog = input("完成后请按'y'\n:")
check = True
while check:
if finish_prog == "y": check = False
plt.close('all')
print("完成!")
quit() # 程序结束
if __name__ == '__main__':
main()
扫一扫
2909
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
