如何使用Obspy在Python中绘制频谱图？（附SAC波形代码）

首先修正你原有波形绘制代码的顺序——plt.show()应该放在st.plot()之后，否则波形图可能无法正常显示：

import matplotlib.pyplot as plt
st = read('SAC file HERE', header=0, index_col=0, parse_dates=True, squeeze=True)
st.plot()
plt.show()

接下来基于这段代码扩展，实现频谱图绘制，核心是对时间序列做傅里叶变换，具体步骤如下：

前提：确认采样频率 #

频谱计算需要知道数据的采样频率（fs，单位Hz）。如果你的read是用pandas读取的格式化数据，可以通过时间索引间隔计算；如果是用地震数据工具（如obspy）读取SAC原生文件，可直接从文件头提取。

方案1：基于pandas+numpy实现 #

假设你用pandas读取的SAC数据已转为一维时间序列：

import matplotlib.pyplot as plt
import numpy as np

# 保持你的数据读取逻辑
st = read('SAC file HERE', header=0, index_col=0, parse_dates=True, squeeze=True)

# 提取波形数据和计算采样频率
data = st.values
time_interval = st.index[1] - st.index[0]
fs = 1 / time_interval.total_seconds()

# 执行FFT计算频谱
n = len(data)
fft_result = np.fft.fft(data)
freq_axis = np.fft.fftfreq(n, 1/fs)

# 过滤正频率并归一化幅度谱
positive_freq_mask = freq_axis >= 0
freq_pos = freq_axis[positive_freq_mask]
amp_pos = np.abs(fft_result[positive_freq_mask]) / n  # 归一化到原始数据幅度范围

# 绘制频谱图
plt.figure(figsize=(10, 6))
plt.plot(freq_pos, amp_pos)
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.title('SAC数据频谱图')
plt.grid(True)
plt.show()

方案2：用obspy处理原生SAC文件（更专业） #

如果你的read是obspy库的读取函数（SAC是地震数据标准格式，obspy对其支持更完善），可以直接从文件头获取采样率，代码更简洁：

from obspy import read
import matplotlib.pyplot as plt
import numpy as np

# 读取SAC文件（obspy原生支持）
st = read('SAC file HERE')
trace = st[0]  # 提取第一个地震道数据

# 绘制波形（修正顺序）
trace.plot()
plt.show()

# 计算并绘制频谱图
data = trace.data
fs = trace.stats.sampling_rate  # 直接从SAC头读取采样率

n = len(data)
fft_result = np.fft.fft(data)
freq_axis = np.fft.fftfreq(n, 1/fs)

positive_freq_mask = freq_axis >= 0
freq_pos = freq_axis[positive_freq_mask]
amp_pos = np.abs(fft_result[positive_freq_mask]) / n

plt.figure(figsize=(10, 6))
plt.plot(freq_pos, amp_pos)
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.title('SAC数据频谱图')
plt.grid(True)
plt.show()

可选优化 #

• 若要绘制功率谱，将amp_pos替换为np.square(amp_pos)即可。
• 为减少频谱泄漏，可先对原始波形加窗处理，例如：data = data * np.hanning(len(data))，再执行FFT。

内容的提问来源于stack exchange，提问作者abdualaziz836