|
|
import numpy as np
import soundfile as sf
import matplotlib.pyplot as plt
# ============================================================
# 1. 参数设置
# ============================================================
# 你的钢琴 WAV 文件
audio_file = "piano.wav"
# 从第几秒开始分析
start_time = 1.0
# 分析多长时间
#
# 如果 WAV 中包含一个单独的钢琴音符，
# 建议先从 0.1 ~ 0.5 秒开始观察。
analysis_duration = 0.2
# FFT 图横轴最多显示到多少 Hz
max_display_freq = 5000
# ============================================================
# 2. 读取 WAV
# ============================================================
audio, sample_rate = sf.read(audio_file)
print("采样率:", sample_rate)
print("原始数据形状:", audio.shape)
# ============================================================
# 3. 如果是立体声，转换为单声道
# ============================================================
# 立体声通常是：
#
# (采样点数量, 2)
#
# 第 0 列：左声道
# 第 1 列：右声道
#
# 这里把左右声道平均。
if audio.ndim == 2:
audio = np.mean(audio, axis=1)
# ============================================================
# 4. 截取需要分析的声音片段
# ============================================================
start_sample = int(start_time * sample_rate)
end_sample = int(
(start_time + analysis_duration) * sample_rate
)
if start_sample >= len(audio):
raise ValueError("start_time 已经超过 WAV 文件长度")
end_sample = min(end_sample, len(audio))
segment = audio[start_sample:end_sample]
print("分析采样点:", len(segment))
print(
"实际分析时间:",
len(segment) / sample_rate,
"秒"
)
# ============================================================
# 5. 绘制时域波形
# ============================================================
time_axis = np.arange(len(segment)) / sample_rate
plt.figure(figsize=(12, 5))
plt.plot(
time_axis,
segment
)
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title("Piano waveform")
plt.grid(
True,
alpha=0.3
)
plt.tight_layout()
plt.show()
# ============================================================
# 6. 使用 Hann 窗
# ============================================================
# 如果直接把一小段声音截出来做 FFT，
# 相当于在截取位置突然把声音切断。
#
# 这种突然截断会造成频谱泄漏。
#
# Hann 窗可以降低这种影响。
window = np.hanning(len(segment))
windowed_signal = segment * window
# ============================================================
# 7. FFT
# ============================================================
# rfft 专门针对实数信号。
#
# 对于采样率 Fs，
# FFT 的有效频率范围是：
#
# 0 Hz ~ Fs/2
#
# Fs/2 就是 Nyquist 频率。
fft_result = np.fft.rfft(
windowed_signal
)
# ============================================================
# 8. 计算每个 FFT 点对应的频率
# ============================================================
freqs = np.fft.rfftfreq(
len(windowed_signal),
d=1 / sample_rate
)
# ============================================================
# 9. 计算幅度
# ============================================================
magnitude = np.abs(
fft_result
)
# 归一化
magnitude = (
magnitude /
(magnitude.max() + 1e-12)
)
# ============================================================
# 10. 转换成 dB
# ============================================================
# dB 表示：
#
# 20 * log10(幅度)
#
# 这样弱小的泛音也更容易看到。
magnitude_db = 20 * np.log10(
magnitude + 1e-12
)
# ============================================================
# 11. 只显示 0 ~ 5000 Hz
# ============================================================
mask = freqs <= max_display_freq
display_freqs = freqs[mask]
display_magnitude_db = magnitude_db[mask]
# ============================================================
# 12. 绘制 FFT 频谱
# ============================================================
plt.figure(figsize=(12, 6))
plt.plot(
display_freqs,
display_magnitude_db
)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.title(
"Piano FFT Spectrum"
)
plt.xlim(
0,
max_display_freq
)
plt.ylim(
-100,
5
)
plt.grid(
True,
alpha=0.3
)
plt.tight_layout()
plt.show() |
|
