Python使用FFT绘制WAV音频频谱与CSV参考结果不符问题求助

现有代码的核心错误 #

普通FFT版本问题 #

1. 未处理单边谱：FFT输出是包含正负频率的双边谱，你直接绘制全长度结果会包含12800Hz以上的镜像冗余数据，同时能量分散在正负频率区间导致幅值偏低
2. 缺少窗函数与修正：无窗函数会导致严重的频谱泄露，加窗后也需要做幅值修正才能得到准确的幅值结果
3. 单位转换逻辑缺失：2e-5是声压（Pa）的参考值，你的FFT输出是信号的量化值，没有转换为实际声压单位就直接计算dB，结果完全没有物理意义
4. 标题与实际使用通道不符：标题写的是First Signal，实际用的是第二通道

Welch版本问题 #

1. 信号缩放完全错误：x = data[1] * (2**15)是将[-1,1]范围的float32信号转为16位整型值的操作，会把信号幅值放大32768倍，对应dB值会偏高约96dB，和你需求完全相反
2. 未使用灵敏度参数做单位转换：和FFT版本问题一致，没有将量化值转为实际声压就计算dB
3. 窗修正缺失：Hann窗会降低信号能量，需要做对应修正才能得到准确幅值
4. 手动加4.5dB没有依据：校准应该通过单位转换逻辑实现，不能随意加偏移值

修正后代码 #

修正版普通FFT代码 #

import numpy as np
import matplotlib.pyplot as plt

sr = 25600
sensitivity = 0.05293 # 52.93mV/Pa 转成V/Pa
full_scale_voltage = 1 # 假设你的float32信号[-1,1]对应采集设备的[-1V,1V]满量程，可根据实际调整

if len(signal.shape) > 1:
    sig = signal[1] # 取第二通道
    n = len(sig)
    # 加Hann窗减少频谱泄露
    window = np.hanning(n)
    sig_windowed = sig * window
    # 计算FFT
    fft_result = np.fft.fft(sig_windowed)
    # 取单边谱（奈奎斯特频率以下部分）
    fft_single = fft_result[:n//2]
    freq = np.fft.fftfreq(n, 1/sr)[:n//2]
    # 幅值修正：单边谱乘2 + 窗函数修正
    magnitude = np.abs(fft_single) * 2 / np.sum(window)
    # 单位转换：量化值 -> 电压 -> 声压
    voltage = magnitude * full_scale_voltage
    pa = voltage / sensitivity
    # 转dB SPL
    spl = 20 * np.log10(pa / 2e-5)

    # 绘图
    plt.figure(figsize=(25, 5))
    plt.semilogx(freq, spl, 'b')
    plt.grid(which='major')
    plt.grid(which='minor', linestyle=':')
    plt.title('Second Channel Full Spectrum')
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("SPL (dB)")
    plt.xlim(20, sr/2) # 只显示有效频率范围
    plt.show()

修正版Welch代码（更推荐用于平滑频谱绘制） #

import numpy as np
from scipy import signal
import librosa
import matplotlib.pyplot as plt

sr = 25600
sensitivity = 0.05293
full_scale_voltage = 1

sig = data[1] # 直接用原始float32信号，不要做缩放
# 窗长调整：频率分辨率=sr/segment_size，2s窗长对应0.5Hz分辨率，可按需调整
segment_size = int(sr * 2)
# Hann窗推荐用75%重叠，平均效果更好
noverlap = int(segment_size * 0.75)

f, Pxx = signal.welch(
    sig,
    fs=sr,
    nperseg=segment_size,
    window='hann',
    nfft=segment_size * 2, # 补零提升频率显示平滑度
    detrend='constant', # 移除直流分量，避免低频偏移
    scaling='spectrum',
    noverlap=noverlap
)
# Hann窗功率谱修正系数
Pxx_corrected = Pxx * 1.633
# 单位转换：功率谱V² -> 声压功率谱Pa²
Pxx_pa = Pxx_corrected * (full_scale_voltage ** 2) / (sensitivity ** 2)
# 转dB SPL，参考声压2e-5Pa对应的功率是(2e-5)^2
ref_pa_sq = (2e-5) ** 2
spl = librosa.power_to_db(Pxx_pa, ref=ref_pa_sq)

# 绘图
plt.figure(figsize=(25, 5))
plt.semilogx(f, spl, 'b')
plt.grid(which='major')
plt.grid(which='minor', linestyle=':')
plt.title('Second Channel Welch Spectrum')
plt.xlabel("Frequency (Hz)")
plt.ylabel("SPL (dB)")
plt.xlim(20, sr/2)
plt.show()

校准调整方法 #

如果修正后和参考CSV仍有整体偏移，调整full_scale_voltage参数即可，不需要手动加dB偏移：

1. 取参考CSV里某一特征频率的dB值作为基准
2. 调整full_scale_voltage数值，直到你绘制的结果中该频率的dB值和参考一致
3. 整体曲线会同步匹配参考结果

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