简介：面向量子振荡数据分析的Python工具包，主要服务凝聚态物理、强磁场输运等研究方向的科研人员与研究生。其围绕Shubnikov-de Haas（SdH）振荡的完整数据处理流程而设计，基于SdHDataSet类对单次磁场扫描的原始与处理数据进行统一管理。实现步骤涵盖：数据导入与清洗、磁场反演与样条插值、扣除多项式磁阻背景、FFT频谱峰识别、对SdH及磁断裂轨道进行滤波分离，并对振幅随逆磁场的变化进行理论拟合，从而提取有效质量、g因子、Dingle温度等关键物理参数。压缩包共4个文件，包括两个可直接调用的Python模块、一个Jupyter Notebook演示样例及一份README说明文档，包体仅367KB，结构紧凑。当前已有175人学习浏览，适合需要系统构建SdH分析流程、快速处理实验数据的物理研究者。借助示例Notebook与峰值检测脚本，用户可快速掌握从数据导入到参数提取的每个环节，并易于迁移到自身测量数据中。 做量子振荡测量的人应该都有这种体验：原始电阻或磁化曲线测出来很漂亮，周期性振荡清清楚楚，可真要从中把费米面极值轨道面积、有效质量、散射率这些物理量干净地挖出来，反而要跟数据处理流程较劲很久。量子振荡数据处理流程代码应用这件事，难点不在于某一个算法有多深，而在于整条链路里有太多细节会悄悄影响最终结果——背景扣不干净，频谱里就全是假峰；磁场轴不重采样，频率峰会整体展宽；窗函数选错，两个靠得很近的频率就再也分不开。这篇文章把我处理SdH和dHvA量子振荡数据时沉淀下来的完整流程、踩过的坑、以及一套可以直接复用的Python代码整理出来，适合刚接触量子振荡数据处理、或者想把整个流程固化成自动化脚本的科研工作者参考。
1. 量子振荡数据处理的完整链条：从"漂亮曲线"到"可信参数"
1.1 一条振荡曲线里到底藏着哪些物理量
量子振荡指的是在强磁场下，材料的电阻（Shubnikov-de Haas振荡）或磁化率（de Haas-van Alphen效应）随1/B呈周期性振荡的现象。它背后连着费米面的拓扑信息，是研究拓扑材料、极低载流子浓度体系、超导体正常态性质的重要实验手段。
处理量子振荡数据，最终目标通常集中在三个物理量上：
- 振荡频率F：由Onsager关系F = (ħ / 2πe) S_F给出，S_F就是费米面极值轨道面积。F是直接从频谱里读出来的，也是整套分析中最基础的一步。
- 有效质量m*：振荡振幅随温度升高而衰减，衰减快慢由Lifshitz-Kosevich公式中的温度因子决定，拟合这个衰减就能得到m*。
- Dingle温度T_D：振荡振幅随1/B增大而指数衰减，衰减速率对应载流子散射率，也就是样品质量的一个表征。
换句话说，振荡信号的"频率、温度依赖、场依赖"三部分信息，分别对应费米面面积、有效质量、散射率。数据处理流程的设计思路，就是想办法把这三个维度的信息干净地分离出来。
1.2 数据处理流水线的整体设计
我习惯把整个过程拆成五步：数据预处理→背景扣除→频谱分析→振幅提取→参数拟合。之所以强调"流程"这个词，是因为这些步骤之间有很强的耦合关系——预处理没做好，背景就扣不干净；背景扣不干净，频谱里就会出现伪峰；伪峰一旦出现，后面拟合出的有效质量基本就是错的。
所以我在给组里学生写自动化脚本时，特别强调一个原则：每个中间步骤都要输出一张图，人眼确认后再进入下一步。后面会提到的大部分"翻车"案例，都是因为跳过了中间检查，直接跑全流程导致结果不可信的。
2. 进入频谱之前：预处理、重采样与坏点处理
2.1 磁场轴非均匀间隔：直接做FFT的隐患
很多实验系统（比如超导磁体配合PPMS）在扫场时，默认是按磁场线性速率扫描的。这意味着采集到的原始数据在B轴上是均匀的，但量子振荡是1/B周期的，最终做傅里叶变换时应该对1/B均匀采样。
如果直接把B轴上的振荡信号做FFT，会发生什么？因为振荡信号在1/B空间里是等周期正弦波，而在B空间里周期会随B变化，也就是"瞬时频率"漂移，直接FFT会把能量展宽到一堆频率上，频率峰变得又矮又胖。我在最初处理数据时就吃过这个亏：拿一段SdH数据直接在B轴上跑FFT，结果频谱里找不到明显的峰，还以为是信号本身太弱。
标准做法是先用1/B = invB 构造均匀网格，把原始数据重采样到这个网格上，再做后续处理。这里有一个细节：原始数据在B轴是均匀的，但B_new = 1/invB_new 是递减的，直接用np.interp时要注意方向。我一般先让invB_new单调递增，再映射回B_new，代码里np.interp要求x单调递增，所以我会写成：
invB_new = np.linspace(invB_raw.min(), invB_raw.max(), 4096) B_new = 1.0 / invB_new data_new = np.interp(invB_new, invB_raw[::-1], data_raw[::-1])意思就是把原始数组倒过来，保证插值点的x轴单调递增。
2.2 坏点剔除与曲线对称化
实验数据里偶尔会有尖峰，可能来自接触热电势抖动、磁体剧烈变化时的感应噪声，也可能仅仅是测量表计瞬间跳变。这些尖峰在直接看曲线时很容易忽略，但重采样加FFT后会在整个频谱上叠加白噪声一样的能量，拉高基线，掩盖弱峰。
我的处理方法是先对原始曲线做一次滑动窗口的中值滤波，或者更直接地，计算相邻点差值的绝对值，把偏离局部几十倍以上的点标记出来做插值替换。这一步不用太精细，目的只是别让个别坏点毁掉整条频谱。
另外，SdH测量如果用的是四探针法，测出来的信号里通常会混入霍尔电压的贡献。严格来说纵向磁阻在磁场反号时是对称的，而霍尔信号是反对称的。所以如果条件允许，我会把正负磁场下的曲线都测了，然后做对称化处理：把正负场数据平均或相减，分离出纯粹的纵向磁阻振荡。没有正负场数据时，至少要在论文里说明这个混叠效应的影响，特别是在低场、高迁移率样品中。
3. 背景扣除：多项式阶数、截取区间和边界震荡
3.1 为什么不能直接把FFT用在原始曲线上
原始磁阻曲线可以看作两部分叠加：一个随磁场缓慢变化的本底背景（主要由经典磁电阻、载流子迁移率温度依赖等贡献），叠加一个高频振荡信号。这个背景在频谱上表现为极低频的大幅分量，如果不去掉，FFT得到的频谱在低频区会有一个巨大的峰，同时由于背景两端不连续，泄漏的能量会污染整个频谱，把真实的振荡峰都淹没掉。
3.2 多项式阶数选择与边界效应规避
背景扣除最经典的方法是多项式拟合。但这里有一个最常见的坑：多项式阶数怎么选。
阶数太低，残差里还留着弯曲的背景；阶数太高，多项式会把一部分振荡当成背景吃掉，导致振荡振幅被严重低估。我的经验是，对大多数SdH曲线，用2到4阶多项式对B拟合就够。5阶以上除非有明确的物理理由，否则不推荐。
判断阶数是否合适的一个直观方法是看残差曲线：如果残差依旧呈"大尺度弯曲"，说明阶数不够；如果残差呈现明显的左右对称、上下基本均匀的振荡，说明背景已经被干净去掉了。另一种更有效的办法是直接在1/B域做高通滤波，或者用Euler方法（对B求导）去除背景，但这会改变振幅信息，后续要做振幅分析时比较麻烦。
边界震荡是另一个容易被忽略的问题。多项式在区间端点附近常常拟合得不好，Fit出的背景曲线在两端会出现"翘起"或"下坠"，扣除后残差两端会出现很大的人工振荡。处理办法是：
- 拟合多项式时，只选取中间一段数据作为拟合窗口，两端各留5%到10%不参与拟合。
- 扣除背景后，把两端边界各裁掉一段，只保留中间振荡信号比较干净的区域。
我个人的习惯是：先用全区间数据做一次多项式和FFT，看看频谱里有没有"红移"的趋势，然后手动调整拟合区间，通常保留磁场窗口的80%左右最稳妥。
4. 傅里叶变换与频率标定：加窗、补零、峰值定位
4.1 用1/B作为变换变量的原因
量子振荡的相位是2πF/B，所以振荡信号在1/B坐标下是严格等周期的正弦波。对1/B做傅里叶变换后，频率轴的单位是特斯拉（T），峰值位置就直接给出振荡频率F。这一步看似简单，但很多人会搞混：如果直接对B做FFT，得到的"频率"单位实际上是1/T，和物理上的F对不上，而且频谱展宽严重。
4.2 窗函数和补零的配合使用
对有限长的振荡信号做FFT，本质上是给信号乘了一个矩形窗，矩形窗的频谱是一个sinc函数，旁瓣很高，会产生频谱泄漏。解决方法是乘一个锥形窗函数，比如Hann窗或Hamming窗，把信号两端削平。
Hann窗的主瓣比矩形窗宽，会略微降低频率分辨率，但旁瓣抑制效果好得多。如果两个频率峰相隔很近，比如一个在35T一个在38T，可以尝试Blackman窗或Kaiser窗（可调beta参数），在频率分辨率和旁瓣抑制之间平衡。我通常默认用Hann窗，多数情况下表现稳定。
补零是另一个常用技巧。在加窗之后把序列后面补一段零，再进行FFT。补零不能提高真实的分辨率（分辨率由数据长度决定），但可以细化频率轴采样，让峰值位置的估计更平滑。一般补零一到两倍长度就够，补太多只会增加计算量，不会带来额外信息。
需要特别注意的是加窗对振幅的影响。加窗会使信号总能量减小，不同窗函数衰减系数不同，所以在提取振幅时必须做归一化修正，也就是把FFT结果的振幅除以窗函数的均值：
from scipy.fft import rfft, rfftfreq window = np.hanning(n_points) signal_windowed = signal * window # 补零到2倍长度 signal_padded = np.concatenate([signal_windowed, np.zeros(n_points)]) spectrum = rfft(signal_padded) freqs = rfftfreq(2 * n_points, d=invB_step) amps = np.abs(spectrum) * 2.0 / (n_points * window.mean())这里除以window.mean()是为了把加窗造成的振幅损失修正回来。
4.3 从频率峰得到费米面极值轨道面积
找到频谱峰位之后，用Onsager关系：S_F = (2πe / ħ) F，就可以算出费米面极值轨道面积。实际中研究者通常直接报告频率F的数值（单位T），因为面积和频率一一对应，换算只是一个单位问题。
频率峰位不应该直接取argmax对应的那个离散频率点，因为离散频谱的采样间隔会带来误差。更稳的方式是对峰附近的几个点做高斯拟合，或者用三点抛物线插值估计真正峰位。我一般对峰两侧各取两三个点做高斯拟合，得到的峰位重复性比直接取最大值好很多。
提到峰值定位，还有一个容易被误导的点：频谱里除了真实振荡峰外，低频区经常会出现一个"零频峰"或基波峰，这是背景没扣干净或者Dingle衰减造成的。不能看见最高峰就当成SdH频率。我通常先看一下频谱整体形态，确认峰的位置是否和预期的载流子口袋面积量级吻合，再决定要读哪个峰。
5. 有效质量与Dingle温度：Lifshitz-Kosevich拟合的约束顺序
5.1 温度依赖项提取有效质量
有效质量的提取依赖Lifshitz-Kosevich公式的温度因子：
R_T = (α T) / sinh(α T)
其中α = 2π² k_B m* / (ħ e B)。这里存在一个B的取值问题：严格说这个因子里的B应该是实际轨道上回旋运动的平均磁场，但通常我们用FFT积分区间内的平均1/B对应的B，也就是调和平均磁场B_avg = 1 / mean(1/B) 来近似。这是领域内很常见的做法，不算严格，但数据窗口选得窄时误差很小。
实际操作流程是：在不同温度下测量同一磁场区间的振荡曲线，分别经过预处理、扣背景、FFT后，取出目标频率峰的振幅A(T)。然后对振幅做温度依赖拟合。注意这里不应该直接拟合振幅绝对值，而是拟合振幅比，把最低温的振幅作为基准，这样可以消掉一些与温度无关的常数因子。
5.2 Dingle项和有效质量参数的耦合问题
Dingle因子R_D = exp(-2π² k_B T_D / (ħ e B / m*))，它随磁场的指数衰减，会让FFT窗口内不同磁场处的振荡振幅不一致。这导致FFT峰的宽度和振幅都会受到Dingle效应影响，尤其是在低磁场端。如果样品散射很强（T_D很高），频谱峰会明显展宽，此时提取的振幅就不纯粹。
还有一个更实际的问题是：当同时拟合m和T_D时，这两个参数会互相"打架"。因为温度因子和Dingle因子在公式里都跟m相关，m大一点、T_D小一点，或者反过来，都可能得到差不多好的拟合效果。正确的顺序是先固定磁场条件，从多温度数据拟合m，再把m*代回去，用单温度数据的场依赖振幅提取T_D。这样耦合效应会小很多。
5.3 拟合参数的初值与边界设置
用scipy.optimize.curve_fit拟合时，参数初值别乱给。m的初值可以先粗略观察：温度从2K升到10K，振幅掉到原来一半左右，m大概在0.3~0.6 m_e的量级；如果几乎不掉，可能m很小（比如0.1以下）。T_D的初值从低场端的振幅衰减速度估计。边界设置上，我通常给m[0.01, 5] m_e，T_D [0, 100] K。这样能防止拟合器跑飞。
拟合结束后的判断也很重要：不只看卡方，还要把拟合曲线叠加到数据点上，人工确认温度依赖的形状是否合理。特别是曲线在高温端是否偏离——如果偏差很大，可能是多频率贡献混叠，或者数据区间里的磁场跨度太大，用单一B_avg近似已经不够理想。
6. 完整代码示例：从模拟数据到物理参数的Pipeline
6.1 模拟SdH数据生成
下面给出一段完整可复现的Python脚本，用模拟数据演示整个流程。模拟数据的好处是读者可以立刻运行、对照结果，也方便调试自己的参数。
import numpy as np from scipy.fft import rfft, rfftfreq from scipy.optimize import curve_fit # 物理常数 e = 1.602e-19 hbar = 1.054e-34 k_B = 1.381e-23 m_e = 9.109e-31 # 模拟参数 B = np.linspace(1.5, 15, 1500) invB = 1.0 / B T_base = 2.0 # 基础温度 m_eff = 0.25 * m_e # 有效质量 T_D = 4.0 # Dingle温度 Freqs = [35.0, 87.0] # 两个振荡频率 phases = [0.0, np.pi/3] # 背景：经典磁电阻 background = 2.0 + 0.12 * B + 0.015 * B**2 # 振荡：L-K振幅 osc = np.zeros_like(B) for F, phi in zip(Freqs, phases): beta = 2 * np.pi**2 * k_B * m_eff / (hbar * e * B) R_T = beta * T_base / np.sinh(beta * T_base) R_D = np.exp(-beta * T_D) osc += R_T * R_D * np.cos(2 * np.pi * F / B + phi) data = background + 0.25 * osc + np.random.normal(0, 0.02, size=B.shape)这里注意振荡振幅系数0.25是经验值，为了让振荡在背景之上明显可见且不淹没在噪声里。实际实验数据里振荡幅度可能只有背景的千分之一，这就要靠后面更精细的背景扣除来处理了。
6.2 预处理与背景扣除实现
# 1/B均匀重采样 invB_new = np.linspace(invB.min(), invB.max(), 4096) B_new = 1.0 / invB_new data_new = np.interp(invB_new, invB[::-1], data[::-1]) # 多项式背景扣除（3阶，仅用中间90%区间拟合） mask = (B_new > B_new.min()*1.05) & (B_new < B_new.max()*0.95) p = np.polyfit(B_new[mask], data_new[mask], 3) background_fit = np.polyval(p, B_new) oscillation = data_new - background_fit这里最关键的一行是np.interp(invB_new, invB[::-1], data[::-1])。因为invB是递减的，插值要求x轴单调递增，所以我把原始数组倒过来，让invB[::-1]变成递增序列，再插值到等间距的invB_new上。很多新手在这一步直接用原始数组，插值结果完全是乱的，却不知道问题出在哪。
6.3 FFT与频率提取实现
# FFT参数 n_points = len(oscillation) window = np.hanning(n_points) osc_windowed = oscillation * window osc_padded = np.concatenate([osc_windowed, np.zeros(n_points)]) sp = rfft(osc_padded) freqs = rfftfreq(2 * n_points, d=invB_new[1] - invB_new[0]) amps = np.abs(sp) * 2.0 / (n_points * window.mean()) # 只保留物理合理的频率区段（比如5到200T） mask_f = (freqs > 5) & (freqs < 200) freq_cut = freqs[mask_f] amp_cut = amps[mask_f] # 用高斯拟合精确定位峰位 from scipy.optimize import curve_fit def gaussian(x, A, mu, sigma): return A * np.exp(-(x - mu)**2 / (2 * sigma**2)) peak_indices = np.argsort(amp_cut)[-2:] # 取两个最高峰 for idx in peak_indices: lo = max(0, idx - 3) hi = min(len(freq_cut), idx + 4) try: popt, _ = curve_fit(gaussian, freq_cut[lo:hi], amp_cut[lo:hi], p0=[amp_cut[idx], freq_cut[idx], 0.5]) print(f"峰位: {popt[1]:.3f} T, 振幅: {popt[0]:.4f}") except RuntimeError: print("局部高斯拟合失败，读取离散最大值") print(f"峰位: {freq_cut[idx]:.3f} T, 振幅: {amp_cut[idx]:.4f}")这段代码里，d=invB_new[1] - invB_new[0]是1/B坐标下的采样间隔，所以FFT频率轴单位的物理含义正好是特斯拉，读出来的峰位就是振荡频率F。这是整个流程中最容易概念混淆的一个点，我在代码注释里特意标了。
6.4 有效质量拟合实现
有效质量拟合需要不同温度下的数据，这里用模拟方式生成多温度曲线，然后提取同一个频率峰在不同温度下的振幅，再拟合温度因子。
# 模拟多温度数据，提取目标频率的振幅 temperatures = np.array([2.0, 3.0, 5.0, 8.0, 12.0]) amplitudes = [] for T in temperatures: osc_T = np.zeros_like(B) for F, phi in zip(Freqs, phases): beta = 2 * np.pi**2 * k_B * m_eff / (hbar * e * B) R_T = beta * T / np.sinh(beta * T) R_D = np.exp(-beta * T_D) osc_T += R_T * R_D * np.cos(2 * np.pi * F / B + phi) data_T = background + 0.25 * osc_T + np.random.normal(0, 0.02, size=B.shape) # 复用前面的预处理、扣背景、FFT流程 data_T_new = np.interp(invB_new, invB[::-1], data_T[::-1]) osc_T_new = data_T_new - np.polyval(np.polyfit(B_new[mask], data_T_new[mask], 3), B_new) w = np.hanning(len(osc_T_new)) sp_T = rfft(np.concatenate([osc_T_new*w, np.zeros(len(osc_T_new))])) amp_T = np.abs(sp_T) * 2.0 / (len(osc_T_new) * w.mean()) # 读取35T附近的振幅 mask_freq = (freqs > 30) & (freqs < 40) idx_peak = np.argmax(amp_T[mask_freq]) amplitudes.append(amp_T[mask_freq][idx_peak]) amplitudes = np.array(amplitudes) # 温度依赖项拟合 B_avg = 1.0 / np.mean(invB_new) def R_T_model(T, alpha_eff): return alpha_eff * T / np.sinh(alpha_eff * T) norm = amplitudes[0] popt, _ = curve_fit(R_T_model, temperatures, amplitudes / norm, p0=[0.5], bounds=[0.05, 5.0]) alpha_eff = popt[0] m_fit = alpha_eff * hbar * e * B_avg / (2 * np.pi**2 * k_B) / m_e print(f"拟合有效质量: {m_fit:.3f} m_e")注意这里拟合的是振幅归一化值amplitudes / norm，这能消掉与温度无关的前置因子，拟合结果更稳健。alpha_eff这个中间量包含了B的贡献，代码里又通过B_avg反推出m*，每一步都有明确的物理对应。
我在实际项目里的习惯是，把这套Pipeline封装成一个类，每个步骤做成独立方法，数据文件路径、磁场窗口、频率搜索范围都做成配置文件。这样新样品来了，只需要改配置文件就能直接批量出结果。日常科研中，数据处理的规范性往往比算法本身更能决定一篇论文数据是否经得起推敲。
本文还有配套的精品资源，点击获取
