python---声压级计算

声压级计算代码如下：

# 在spl.py文件里面的计算声压和响度函数
import numpy as np
import math
#np.seterr(divide='ignore', invalid='ignore')
def spl_en(x,fs,flen):
    '''
    :param x:代表输入的语音信号
    :param fs: 采样率
    :param flen: 一帧信号的时间长度，单位是毫秒
    :return: 声压级的值
    '''
    length = len(x)  # 得到输入语音的长度
    #每一帧信号的离散的点数
    m = int(flen * fs/1000)
    if length != m:
        print('length',length)
        print('m',m)
        print('输入信号长度与所定义帧长不等！')

    # 计算声压级根据定义计算有效声压，pa = sqrt((x(1)^2+x(2)^2+...+x(M)^2)/M)
    # 单位是Pa

    pp = 0
    for i in range(0,m-1):
        pp += x[i] ** 2

    pa = np.sqrt(pp / m)

    # 计算声压级，声压级值spl=20*log10(pa/p0),单位为dB
    #基准声压p0，单位为Pa
    p0 = 2 * 1e-5
    spl_value = 20*np.log10(pa / p0)

    return spl_value


# 计算响度函数
def iso226(phon):
    f = np.array([20 ,25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500, 630, 800,
         1000 ,1250 ,1600, 2000, 2500, 3150, 4000, 5000, 6300, 8000, 10000, 12500])

    af = np.array([0.532, 0.506, 0.480, 0.455, 0.432 ,0.409 ,0.387 ,0.367, 0.349, 0.330, 0.315,
          0.301, 0.288, 0.276, 0.267, 0.259, 0.253, 0.250, 0.246, 0.244, 0.243, 0.243,
          0.243, 0.242, 0.242, 0.245, 0.254, 0.271, 0.301])

    Lu =np.array([-31.6 ,- 27.2, - 23.0, - 19.1 ,- 15.9, - 13.0 ,- 10.3 ,- 8.1 ,- 6.2 ,- 4.5, - 3.1,
          - 2.0 ,- 1.1, - 0.4,   0.0 ,  0.3 ,  0.5 ,  0.0 ,- 2.7, - 4.1 ,- 1.0,  1.7,
          2.5 ,  1.2 ,- 2.1 ,- 7.1 ,- 11.2, - 10.7 ,- 3.1])

    Tf = np.array([78.5 , 68.7 , 59.5 , 51.1 , 44.0 , 37.5 , 31.5,  26.5 , 22.1 , 17.9,  14.4,
          11.4 ,  8.6 ,  6.2 ,  4.4,   3.0,   2.2 ,  2.4 ,  3.5 ,  1.7 ,- 1.3 ,- 4.2,
          - 6.0, - 5.4, - 1.5,   6.0 , 12.6 , 13.9,  12.3])

    # 错误处理

    if ((phon < 0) or (phon > 90)):
        print('因素值超出了范围！')
        spl=0
        freq=0
    else:
        Ln = phon
        #从响度级计算声压级
        Af = 4.47 * math.e -3 * (10*(0.025*Ln) - 1.15) + (0.4*10*(((Tf+Lu)/10)-9 ))*af
        Lp = ((10 / af)* np.log10(np.abs(Af))) - Lu + 94  # 取对数这里是不能出现复数的，所以要加绝对值

        spl = Lp
        freq = f

    return [spl,freq]

其他调用函数：

#在spl_jisuan.py文件里面
from spl import spl_en
import numpy as np
from matplotlib import  pyplot as plt
import librosa


data,fs = librosa.load('C2_2_y.wav',sr=None)

length = len(data)

# % 每帧大小为M，当语音长度不是帧长的整数倍时：
# % (1)若剩余长度大于等于帧长的二分之一，则补零至帧长
# % (2)若剩余长度小于帧长的二分之一，则舍弃

framlen = 100
#每帧信号的离散点数
M = int(fs * framlen/1000)   # 一帧有多少点数
m = np.mod(length,M)  # 求余数
if m >= M//2:  #补零
    x = np.pad(data,M-m,mode='constant',constant_values=0)
    length = len(x)

else:
    l = np.floor(length//M)  # 下取整函数,要舍弃
    x = data[0:M*l]  # 选择M*l个点数，刚好把多余的部分截断了
    length = len(x)

N = length // M  #最终的分帧数

s = np.zeros((1,M))  # 一帧存储M个点数
spl_f = np.zeros((1,int(N)))  # 分帧数
for i in range(1,int(N)):
    s = x[(i-1)*M:i*M]  # 一次一帧个点数
    spl_f[0,i] = spl_en(s,fs,framlen)  #算spl,存放到数组里面，每一帧算出一个标量

t = np.arange(0,length)   # 语句的长度，选取时间的长度
SPL = np.zeros((1,length))
for r in range(1,int(N)):  # 把结果存放到对应的SPL数组里面
    SPL[0,(r-1)*M:r*M]=spl_f[0,r]

plt.figure()
plt.subplot(2,1,1)
plt.plot(t/fs,x)  # 算出时间是多少秒，对应的是x的语句的长度，这个是截断以后的长度
plt.xlabel('time/s')
plt.ylabel('wave spec')
plt.subplot(2,1,2)
plt.step(t/fs,SPL[0],'r')  # 时间，算出来的SPL
plt.xlabel('time/s')
plt.ylabel('db')
plt.title('yuyin xinhao de shengyaji db')
plt.show()

#在com_iso.py文件里面
from spl import spl_en,iso226
import numpy as np
from matplotlib import  pyplot as plt
import librosa

phon=50
spl = iso226(phon)
plt.figure()  # 使用plt.subplots()的时候，不需要添加plt.figure()也是可以的，
# 或者fig=plt.figure()
# ax=fig.add_subplots(121)
plt.semilogx(spl[1],spl[0],':','color','k')
# plt.xlim(20,20000)
# plt.ylim(-10,130)
# plt.axis([20,20000,-10,130])
y_ticks = range(-10,200,10)
plt.yticks(y_ticks,y_ticks,rotation=60)
plt.xlabel('freq/hz')
plt.ylabel('shengya  jibie  db')
plt.title('phone=50')
#plt.yticks(np.arange(start=-10,stop=130,step=10),np.arange(start=-10,stop=130,step=10),rotation=60)
plt.grid(b=True,linestyle='--',linewidth=0.5)  # 是否显示网格线
plt.box()  # 添加边框，还有其他的用法也是可以用的。plt.gca()创建这个对象，再条用spine的一系列属性设置值。
plt.show()

结果是：

 

 显示第二张图的时候，是不能用plt.xlim()函数设置x轴的范围的，我也不太明白为什么，如果设置了，还是现在这样子的，但是半对数曲线就显示不出来了，可能是因为范围的问题吧。

用到的参考资料为：

Matplotlib绘图（一）-边框线及坐标轴的设置_帅兄心安否的博客-CSDN博客_matplotlib 边框https://blog.csdn.net/weixin_44560088/article/details/107165510matplotlib之pyplot模块——显示/关闭子图边框（轴脊）（box()）_mighty13的博客-CSDN博客_matplotlib 子图边框https://blog.csdn.net/mighty13/article/details/116571742声压与声压级的计算_jianghuzhijian的博客-CSDN博客_声压级计算https://blog.csdn.net/jianghuzhijian/article/details/103798289