我们将涵盖从基本的数据加载、提取特定波长/强度值，到更高级的寻峰、积分和绘图等操作。
（图片来源网络，侵删）
准备工作：安装必要的库
你需要安装一些核心的科学计算和数据处理库，最常用的是 NumPy, Pandas, Matplotlib 和 SciPy。
pip install numpy pandas matplotlib scipy
- NumPy: 用于高效的数值计算，特别是处理多维数组（如光谱数据）。
- Pandas: 用于处理类似表格的数据，可以方便地加载和操作 CSV 或 Excel 格式的光谱数据。
- Matplotlib: 用于绘制各种图表,包括光谱图。
- SciPy: 提供了许多科学计算工具，
scipy.signal模块在信号处理（如寻峰）中非常有用。
加载光谱数据
光谱数据通常以文本文件（如 .txt, .csv）或 Excel 文件（.xlsx）的形式存储，其中包含两列：波长和强度。
示例数据格式 (spectrum_data.csv)
Wavelength (nm),Intensity (a.u.)
350, 10
351, 12
352, 15
...
700, 5
701, 3
702, 2使用 Pandas 加载（推荐）
Pandas 提供了非常方便的 read_csv 或 read_excel 函数，可以直接将数据读入一个 DataFrame 中,这是一种结构化的表格数据类型。
import pandas as pd
# 从 CSV 文件加载数据
# 假设你的文件名为 'spectrum_data.csv'，并且第一行是标题
df = pd.read_csv('spectrum_data.csv')
# 查看数据的前几行
print("数据预览:")
print(df.head())
# 获取波长和强度的列数据
wavelengths = df['Wavelength (nm)'].values # .values 将其转换为 NumPy 数组
intensities = df['Intensity (a.u.)'].values
print(f"\n波长数组: {wavelengths[:5]}...") # 打印前5个元素作为示例
print(f"强度数组: {intensities[:5]}...")
使用 NumPy 加载
如果你的文件格式非常简单，只有两列数据，没有标题行，NumPy 的 loadtxt 函数更直接。
（图片来源网络，侵删）
import numpy as np
# 从文本文件加载数据
# 假设文件没有标题行，并且是空格或逗号分隔
data = np.loadtxt('spectrum_data.txt', delimiter=',') # 如果是CSV，用delimiter=','
# 将数据分成两列
wavelengths = data[:, 0]
intensities = data[:, 1]
print(f"波长数组: {wavelengths[:5]}...")
print(f"强度数组: {intensities[:5]}...")
核心提取操作
一旦你有了 wavelengths 和 intensities 这两个数组,就可以进行各种提取操作了。
操作 1：提取特定波长范围内的光谱
这是最常见的操作之一，只关心 400nm 到 500nm 的光谱区域。
# 定义你感兴趣的波长范围
min_wavelength = 400
max_wavelength = 500
# 使用布尔索引进行筛选
# (wavelengths >= min_wavelength) 会返回一个布尔数组
# (wavelengths <= max_wavelength) 也是
# 两者相与 (&)，得到一个满足条件的布尔数组
mask = (wavelengths >= min_wavelength) & (wavelengths <= max_wavelength)
# 应用这个掩码来提取对应的数据
extracted_wavelengths = wavelengths[mask]
extracted_intensities = intensities[mask]
print(f"\n在 {min_wavelength}-{max_wavelength} nm 范围内的数据点数: {len(extracted_wavelengths)}")
操作 2：提取特定波长值对应的光谱强度
如果你想查找某个精确波长处的强度值。
target_wavelength = 550.0
# 方法 A: 直接索引（如果你的波长是等间隔的，且数组已排序）
# 这种方法不推荐，因为很难精确匹配浮点数
# index = np.where(wavelengths == target_wavelength)
# intensity_at_target = intensities[index]
# 方法 B: 插值（推荐！）
# 真实世界的数据点往往是离散的，插值可以得到更精确的值
from scipy.interpolate import interp1d
# 创建一个插值函数
# kind='linear' 是线性插值，'cubic' 是三次样条插值，更平滑
f = interp1d(wavelengths, intensities, kind='linear')
# 使用这个函数计算目标波长的强度
intensity_at_target = f(target_wavelength)
print(f"\n在波长 {target_wavelength} nm 处的插值强度为: {intensity_at_target:.2f}")
操作 3：提取峰值信息
识别光谱中的峰位和峰高是分析的关键步骤。
from scipy.signal import find_peaks
# 寻找峰值
# prominence: 峰的显著性，用于过滤掉小的噪声峰
# distance: 峰之间的最小距离，避免识别到同一个峰的多个点
peaks, properties = find_peaks(intensities, prominence=50, distance=10)
# 提取峰位和峰高
peak_wavelengths = wavelengths[peaks]
peak_intensities = intensities[peaks]
print("\n找到的峰值信息:")
for i in range(len(peak_wavelengths)):
print(f" 峰 {i+1}: 波长 = {peak_wavelengths[i]:.2f} nm, 强度 = {peak_intensities[i]:.2f}")
操作 4：提取某个峰的积分面积（峰面积）
峰面积通常比峰高更能代表该物质的含量。
# 假设我们想计算第一个峰的面积
# 我们需要知道这个峰的起始和结束波长
# find_peaks 的 properties 中可以找到 'left_ips' 和 'right_ips' (峰的左右边界点索引)
if len(peaks) > 0:
first_peak_index = peaks[0]
left_ip_index = properties['left_ips'][first_peak_index]
right_ip_index = properties['right_ips'][first_peak_index]
# 提取峰范围内的数据
peak_area_wavelengths = wavelengths[int(left_ip_index):int(right_ip_index)+1]
peak_area_intensities = intensities[int(left_ip_index):int(right_ip_index)+1]
# 使用梯形法则计算面积
from scipy.integrate import trapz
peak_area = trapz(peak_area_intensities, peak_area_wavelengths)
print(f"\n第一个峰的积分面积为: {peak_area:.2f}")
可视化
将提取的结果绘制出来是验证和展示分析结果的最佳方式。
import matplotlib.pyplot as plt
# 创建一个图形和坐标轴
fig, ax = plt.subplots(figsize=(10, 6))
# 1. 绘制原始光谱
ax.plot(wavelengths, intensities, label='原始光谱', color='gray', linestyle='--')
# 2. 高亮显示提取的波长范围
ax.fill_between(extracted_wavelengths, 0, extracted_intensities,
alpha=0.3, color='blue', label='提取范围 (400-500 nm)')
# 3. 标记找到的峰值
ax.plot(peak_wavelengths, peak_intensities, "x", color='red', markersize=10, label='检测到的峰')
# 添加标签和标题
ax.set_xlabel('波长 (nm)')
ax.set_ylabel('强度 (a.u.)')
ax.set_title('光谱数据提取与可视化')
ax.legend()
ax.grid(True)
# 显示图形
plt.show()
完整代码示例
下面是一个将上述所有步骤整合在一起的完整脚本。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.interpolate import interp1d
from scipy.integrate import trapz
# --- 1. 加载数据 ---
# 假设有一个名为 'spectrum_data.csv' 的文件
try:
df = pd.read_csv('spectrum_data.csv')
wavelengths = df['Wavelength (nm)'].values
intensities = df['Intensity (a.u.)'].values
print("数据加载成功！")
except FileNotFoundError:
print("错误：未找到 'spectrum_data.csv' 文件。")
# 创建一些模拟数据以便演示
print("正在生成模拟数据...")
wavelengths = np.linspace(350, 750, 1000)
intensities = np.exp(-((wavelengths - 450)**2) / (2 * 20**2)) * 1000 + \
np.exp(-((wavelengths - 600)**2) / (2 * 30**2)) * 800 + \
np.random.normal(0, 50, len(wavelengths))
intensities = np.maximum(intensities, 0) # 确保强度不为负
# --- 2. 提取操作 ---
# 提取特定范围
min_w, max_w = 400, 500
mask = (wavelengths >= min_w) & (wavelengths <= max_w)
extracted_w = wavelengths[mask]
extracted_i = intensities[mask]
# 插值查找特定波长值
target_w = 550
f_interp = interp1d(wavelengths, intensities, kind='linear')
intensity_at_target = f_interp(target_w)
# 寻找峰值
peaks, _ = find_peaks(intensities, prominence=100, distance=20)
peak_ws = wavelengths[peaks]
peak_is = intensities[peaks]
# 计算第一个峰的面积
if len(peaks) > 0:
# 获取峰的边界索引
from scipy.signal import peak_widths
widths, width_heights, left_ips, right_ips = peak_widths(intensities, peaks, rel_height=0.5)
# 选择第一个峰
left_idx = int(left_ips[0])
right_idx = int(right_ips[0])
peak_area_w = wavelengths[left_idx:right_idx+1]
peak_area_i = intensities[left_idx:right_idx+1]
peak_area = trapz(peak_area_i, peak_area_w)
# --- 3. 可视化 ---
fig, ax = plt.subplots(figsize=(12, 7))
# 绘制原始光谱
ax.plot(wavelengths, intensities, label='原始光谱', color='gray', alpha=0.7)
# 高亮提取范围
ax.fill_between(extracted_w, 0, extracted_i, alpha=0.4, color='blue', label=f'提取范围 ({min_w}-{max_w} nm)')
# 标记目标波长和插值强度
ax.plot(target_w, intensity_at_target, 'go', markersize=8, label=f'插值点 @ {target_w}nm')
# 标记所有峰值
ax.plot(peak_ws, peak_is, 'rx', markersize=10, label='检测到的峰')
# 标记第一个峰的面积
if len(peaks) > 0:
ax.fill_between(peak_area_w, 0, peak_area_i, alpha=0.5, color='red', label=f'第一个峰面积 = {peak_area:.0f}')
# 美化图形
ax.set_xlabel('波长 (nm)')
ax.set_ylabel('强度 (a.u.)')
ax.set_title('Python 光谱数据提取与分析示例')
ax.legend()
ax.grid(True, linestyle='--', alpha=0.6)
plt.tight_layout()
plt.show()
# --- 4. 打印结果 ---
print("\n--- 提取结果 ---")
print(f"在 {min_w}-{max_w} nm 范围内共有 {len(extracted_w)} 个数据点。")
print(f"波长 {target_w} nm 处的插值强度为: {intensity_at_target:.2f}")
print(f"共检测到 {len(peaks)} 个峰值。")
if len(peaks) > 0:
print(f"第一个峰的积分面积为: {peak_area:.2f}")
总结与建议
- 数据格式是关键：首先确认你的光谱数据格式，是 CSV、TXT 还是其他，这决定了你使用
pandas还是numpy来加载。 - Pandas 是你的朋友：对于大多数带有标题和复杂格式的数据，使用 Pandas 的
DataFrame是最方便、最不容易出错的。 - 插值优于直接索引：当你需要查找非数据点处的值时，使用
scipy.interpolate.interp1d是更科学、更准确的方法。 - SciPy 是信号处理的利器：
scipy.signal模块中的find_peaks和相关函数（如peak_widths）是进行光谱峰分析的强大工具。 - 可视化是最好的验证：每一步提取操作后，都建议将结果绘制出来,这样可以直观地检查你的操作是否正确。
希望这个全面的指南能帮助你顺利地使用 Python 进行光谱数据的提取和分析！
