总声压级计算与1/3倍频程声压级计算

    最近在研究声压级计算的算法，说实话网上搜了好多的资料，大都是将某个声压值转换为DB的公式，这公式随便一本书上就有，哪里用得着您细讲嘞？

    废话少说，建议想看懂我这篇博客的先认着找资料研究一下什么是1/3倍频程，以及声压级计算的基本公式。在计算1/3倍频程的声压级时，要搞清楚傅里叶变换的来世今生，目的是使频域的幅值精确对应上。本博客在进行倍频程计算时选择了基数2来确定倍频程带的上下频率。计算声压级最重要的一句话：各个倍频程带内的声压均方值是该频带内频谱各谱线幅值的均方值之和。

    解释一下，均方值是一个统计量，一根谱线的均方值是什么呢？这就要将这根谱线对应到时域信号求其均方值了，也就是这根谱线的幅值除以根号2。而总声压均方值则是各个倍频带内的均方值之和。下面附上我计算时的代码：

import numpy as np
from matplotlib import pyplot as plt
from numpy.fft import fft
import math
 
 
def fft_custom(sig, Fs, zero_padding=False, window=None):
    if zero_padding:
        fft_num = np.power(2, math.ceil(np.log2(len(sig))))
        if window:
            mag = np.fft.fft(sig * np.hanning(len(sig)), fft_num) * 2 / fft_num
            f = np.arange(fft_num) * (Fs / fft_num)
        else:
            mag = np.fft.fft(sig, fft_num) * 2 / fft_num
            f = np.arange(fft_num) * (Fs / fft_num)
    else:
        fft_num = len(sig)
        if window:
            mag = np.fft.fft(sig * np.hanning(len(sig)), fft_num) * 2 / fft_num
            f = np.arange(fft_num) * (Fs / fft_num)
        else:
            mag = np.fft.fft(sig, fft_num) * 2 / fft_num
            f = np.arange(fft_num) * (Fs / fft_num)
    mag[0] /= 2
    return f[:int(fft_num / 2)], abs(mag[:int(fft_num / 2)])
 
 
def get_octave_band_base_2(start=-23, end=14):
    """以2为基数返回1/3倍频程带"""
    bands = list()
    for index, i in enumerate(range(start, end)):
        central_frequency = 1000 * (2 ** (1 / 3)) ** i
        if index == 0:
            bands.append([0, central_frequency * np.power(2, 1 / 6)])
        else:
            bands.append([central_frequency / np.power(2, 1 / 6), central_frequency * np.power(2, 1 / 6)])
    return np.asarray(bands)
 
 
if __name__ == '__main__':
    # 仿真信号
    # Fs = 10000
    # N = 50000
    # n = list(np.arange(N))
    # t = np.asarray([temp / Fs for temp in n])
    # sig = np.sqrt(2) * np.sin(2 * np.pi * 5 * t) + np.sin(2 * np.pi * 18 * t) + np.sin(2 * np.pi * 53 * t) + np.sin(
    #     2 * np.pi * 197 * t)
 
    # 实际信号
    path = r'C:\Users\Administrator\Desktop\yourfilepath.txt'
    Fs = 20833
    sig = np.loadtxt(path)
    N = len(sig)
    n = list(np.arange(N))
    t = np.asarray([temp/Fs for temp in n])
 
    # 1 时域内计算声压级
    temp = [temp * temp for temp in sig]
    p_square_time = np.sum(temp) / len(temp)  # 时域内计算的声压总均方值
    print("p_square_time: ", p_square_time)
 
    f, mag = fft_custom(sig, Fs)
    # 频域内计算声压总均方值
    temp = [(temp / np.sqrt(2)) * (temp / np.sqrt(2)) for temp in mag]
    p_square_frequency = np.sum(temp)
    print("p_square_frequency: ", p_square_frequency)
 
    spl_overall = 10 * np.log10(p_square_time / 0.0000000004)
 
    print("声压级(DB):", spl_overall)
 
    bands = get_octave_band_base_2(start=-21, end=15)  # start 和 end 控制带宽
    spl_of_bands = list()
    f = list(f)
    index_start = 0
    for band in bands:
        index_stop = np.where(f < band[1])[0][-1]
        band_frequencies_mag = mag[index_start:index_stop]
        temp = [(temp / np.sqrt(2)) * (temp / np.sqrt(2)) for temp in band_frequencies_mag]
        p_square_frequency_band = np.sum(temp)
        spl_ = 10 * np.log10(p_square_frequency_band / 0.0000000004)
        if band[0] <= Fs/2:
            spl_of_bands.append([band[0], band[1], spl_])
        else:
            break
        # print("index_start: %d index_stop: %d" % (index_start, index_stop))
        index_start = index_stop
        # print(band)
 
    plt.figure(1)
    plt.subplot(211)
    plt.plot(t, sig)
    plt.subplot(212)
    plt.plot(f, mag)
    # plt.show()
 
    spl_values = list()
    xticks = list()
 
    for temp in spl_of_bands:
        spl_values.append(temp[2])
        xticks.append(str(int(temp[1])))
 
    plt.figure(2)
    plt.bar(range(len(spl_values)), spl_values, facecolor='b', width=0.8)
    plt.title("1/3 octave SPL || Overall SPL is %f " % np.round(spl_overall, 3))
    plt.xticks(list(range(len(spl_values))), xticks,rotation=45)
    plt.show()

    直接给大家看下运行该代码的效果吧：

       大家也可以用我提供的仿真信号来验证该声压级计算算法的准确性，欢迎交流~