绘制1/3倍频程频谱图

对wav格式文件，绘制1/3倍频程频谱图，频率范围为20-10000Hz。

from scipy.io import wavfile
from scipy.fft import fft
import numpy as np
import matplotlib.pyplot as plt
import librosa
#from scipy.signal import hamming


def one_third_octave(audio_path):
    
    sample_rate, audio_data = wavfile.read(audio_path)
    
    # 确保音频数据是一维的（单声道）
    if len(audio_data.shape) == 2:
        audio_data = audio_data.mean(axis=1)

   # 计算FFT
    n = len(audio_data) # 数据长度为240000
    fft_audio = fft(audio_data, n) # 240000行×1列 复数
    fft_audio = fft_audio[range(n//2)]  # 只取一半，120000行×1列 复数
    
    # 计算功率谱
    power_spectrum = (2.0/n) * np.abs(fft_audio)**2 

    # 计算频率数组
    freqs = float(sample_rate) / n * np.arange(n//2) # 120000行× 1列

    # 计算1/3倍频程带的频率

    fmin = 20  # 最小频率
    fmax = sample_rate / 2  # 最大频率
    n_bins = 29  # 1/3倍频程带的数量
    bins_per_octave = 3  # 每个倍频程的1/3倍频程带数量
    n_octaves = int(np.ceil(np.log2(fmax / fmin)*3))  # 倍频程的数量

    # 计算1/3倍频程带的中心频率
    freqs13 = librosa.cqt_frequencies(n_bins=n_bins, fmin=fmin, bins_per_octave=bins_per_octave)

    # 计算1/3倍频程的边界

    one_third_octaves = freqs13
    # 计算每个1/3倍频程的功率值
    one_third_octave_power = np.zeros(len(one_third_octaves) - 1)
    for i in range(len(one_third_octaves) - 1):
        mask = (freqs >= one_third_octaves[i]) & (freqs <= one_third_octaves[i+1])
        one_third_octave_power[i] = np.trapz(power_spectrum[mask], freqs[mask])

    # 绘制1/3倍频程频谱图
    octaves_13 = np.round(one_third_octaves[:-1],decimals=0) # 取1/3倍频第0-倒数第1个数（0-25），并对元素四舍五入，保留0位小数
    octaves_13 = octaves_13.astype(int) # 转化为整型，去掉小数点
    str_octaves_13 = [str(item) for item in octaves_13] # 将数组转化为字符串列表

    # 绘图
    plt.figure(figsize=(20, 6))
    plt.bar(str_octaves_13, 10 * np.log10(one_third_octave_power),width=0.6)
    plt.title('1/3 Octave Frequency Spectrum')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Power (dB)') #plt.xscale('log') #设置当前轴（通常是x轴）的缩放类型。
    #plt.grid(True, which="both", ls="--")
    plt.show()

# 加载音频文件
audio_path = '2023-03-20-10-34-59.wav'
one_third_octave(audio_path)

结果展示：

备注：librosa.cqt_frequencies生成的1/3倍频程频率结果有偏差，暂未找到原因，欢迎留言指教！