WAV局部频域绘制：自动截段、分频段平滑及dB轴转换问题

背景 #

我需要绘制音频样本的频域图，对信号处理一无所知，Python基础薄弱，参考示例后用fft和fftfreq实现初步效果。所有音频时长约30秒，含两个相同极短脉冲声，需在浏览器绘图，仅保留相关片段少量数据点。当前已实现的完整代码如下：

def doFFT(dataw: ndarray, samplerate: int): 
    dataw = dataw[(round(24.98*samplerate)):(round(24.98*samplerate)+32767)]
    datafft = rfft(dataw)  
    # 获取实部和虚部的绝对值
    fftabs = abs(datafft)     
    freqs = rfftfreq(dataw.shape[0], 1/(samplerate))  
    smoothed = savgol_filter(fftabs,  10,  3)  
 
    return [ freqs[:int(freqs.size/2)].tolist() , smoothed[:int(freqs.size/2)].tolist() ]

将结果转为JSON后用Plotly在浏览器绘图，当前效果如图：

问题与解决方案 #

问题1：手动选取时间段方法粗糙，如何自动定位脉冲片段？ #

当前手动截取代码：

samplerate, dataw= wavfile.read('x.wav')
dataw = dataw[(round(24.98*samplerate)):(round(24.98*samplerate)+32767)]

解决方案：
通过信号能量自动定位脉冲，再截取有效片段：

import numpy as np
from scipy.io import wavfile

samplerate, dataw = wavfile.read('x.wav')
# 多声道转单声道
if len(dataw.shape) > 1:
    dataw = np.mean(dataw, axis=1)
# 计算信号能量
energy = np.abs(dataw)
# 设定阈值（可根据音频实际情况调整，比如取能量均值的5倍）
threshold = 5 * np.mean(energy)
# 找出超过阈值的脉冲索引
pulse_indices = np.where(energy > threshold)[0]

if len(pulse_indices) > 0:
    # 扩展截取范围，保留脉冲前后足够数据
    start_idx = max(0, pulse_indices[0] - 1024)
    end_idx = min(len(dataw), pulse_indices[-1] + 1024)
    dataw = dataw[start_idx:end_idx]

这种方法无需硬编码时间点，能自动适配不同音频的脉冲位置。

问题2：如何实现分区间不同平滑因子的频谱平滑？ #

朋友的Matlab分区间平滑代码：

%% 划分频谱三个区间，用不同平滑因子提升可读性
smooth_spect_1 = fastsmooth(Spectrum_dB(1:70), 5, 1, 1); 
smooth_spect_2 = fastsmooth(Spectrum_dB(71:700), 7, 1, 1);         
smooth_spect_3 = fastsmooth(Spectrum_dB(701:end), 41 , 1, 1);        
total_smooth = [smooth_spect_1; smooth_spect_2; smooth_spect_3];
plot(f, total_smooth);

当前用单一savgol_filter处理全频谱，高频平滑不足。

解决方案：
在Python中对频谱按频率区间拆分，分别应用不同平滑参数后拼接：

from scipy.signal import savgol_filter
import numpy as np

# 假设已得到freqs（频率数组）和fftabs（频谱幅值）
# 先转dB（提前处理问题3）
spect_dB = 20 * np.log10(fftabs / np.max(fftabs))
spect_dB = np.where(fftabs == 0, -np.inf, spect_dB)

# 定义频率区间与对应平滑窗口大小
freq_ranges = [
    (0, 700),    # 低频区间
    (700, 7000), # 中频区间
    (7000, None) # 高频区间
]
window_sizes = [5, 7, 41]

smoothed_parts = []
filtered_freqs = []
for (low, high), win_size in zip(freq_ranges, window_sizes):
    # 筛选对应频率区间的索引
    if high is None:
        mask = freqs >= low
    else:
        mask = (freqs >= low) & (freqs < high)
    # 平滑对应区间
    smoothed_part = savgol_filter(spect_dB[mask], win_size, 1)
    smoothed_parts.append(smoothed_part)
    filtered_freqs.append(freqs[mask])

# 拼接结果
total_smooth = np.concatenate(smoothed_parts)
final_freqs = np.concatenate(filtered_freqs)

注：频率区间数值需根据采样率和实际频谱范围调整，若偏好按索引划分，直接替换为索引切片即可。

问题3：如何将Y轴转换为dB单位？ #

解决方案：
采用音频频谱转dB的标准公式：dB = 20 * log10(幅值 / 参考幅值)，示例代码：

import numpy as np

# 方法1：归一化到频谱峰值（0dB为最大值）
spect_dB = 20 * np.log10(fftabs / np.max(fftabs))
# 处理0值避免log报错
spect_dB = np.where(fftabs == 0, -np.inf, spect_dB)

# 方法2：基于音频采样位数（如16位音频参考值为32767）
# reference = 32767
# spect_dB = 20 * np.log10(fftabs / reference)

将转换后的spect_dB作为Y轴数据传入Plotly，即可实现dB刻度显示。

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