首先Overall Level 计算的是RMS,但它就算的是频率谱的均方根(RMS)值,计算步骤如下:
对噪音或振动进行fft变换,结果如下图:
计算公式为:
注意事项:进行计算的时候需要是单边谱,如果存在直流分量时A0分量不需要除以2;
因为阶次分析实际上也是部分频率带宽内的能量计算,上述公式同样适用于阶次提取。频率带宽的定义有下图三种模式,按需定义即可。
对于Overall Level 是不是与手持声级计测量的值是否相同?答案是否定的。
根据 IEC 61672-1定义标准声级计的工作原理如下:
-
标准声级计被称为指数平均声级计,因为它通过一个均方根(RMS)电路将麦克风的交流信号转换为直流信号。
-
为了进行这种转换,声级计需要一个积分的时间常数,称为时间加权。这影响了声级计对声音信号变化的响应速度。
-
标准化的时间加权:
S(Slow):时间加权为1秒,适合测量较慢变化的声音。
F(Fast):时间加权为125毫秒,适合测量快速变化的声音。
I(Impulse):时间加权为35毫秒,最初用于测量冲击噪声。
如何理解指数平均和时间加权是本文的重点内容:
假设声音信号为y(t),时间加权的公式如下:
式中
τ 时间加权的时间常数,单位为秒s;
p0 参考声压级 2*10^-5 pa
h(y2) 实际上是将声音信号进行平方运算,然后用脉冲响应为
的滤波器进行卷积计算。(卷积公式:
)利用拉普拉斯变换的该滤波器s域函数
, 这是一个模拟滤波器,通过双线性变换可得Z域函数:
,便获得了数字滤波器。实际上还是通过滤波器实现了时间加权。不同的时间常数只是影响滤波特性罢了。 下图是用AMEsim 和代码对不同时间常数的时间加权滤波器的仿真,时间常数越大滤波效果越明显:
为什么说声级计是指数平均声级计?
因为Z域函数:
的差分方式为 y[n] = (1 - a) * y[n-1] + a * x[n],
a=
观察上式:时间计权等效滤波器的差分方程就是指数平均方程,指数平均和时间加权其实是一个事儿,只是观察角度不同。
实践部分:计算 A 计权的 Fast 声压级 LAF,流程如下:
声压信号
LAF
代码 :
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import zpk2tf, freqs, bilinear_zpk, zpk2sos, lfilter
import math
# 参考声压,单位为 Pa
P_ref = 20e-6
# 模拟输入声压信号 (单位: Pa)
# 这里假设声压信号为 1 秒钟内的采样数据
fs = 1000 # 采样率为1000 Hz
t = np.linspace(0, 1, fs) # 1秒钟内的时间点
pressure_signal = 0.05 * np.sin(2 * np.pi * 100 * t) # 模拟声压信号 (100Hz)
# 时间加权函数 (指数平滑滤波器)
def time_weighting(signal, fs, tau):
alpha = 1 - np.exp(-1 / (fs * tau))
weighted_signal = np.zeros_like(signal)
weighted_signal[0] = signal[0]
for i in range(1, len(signal)):
weighted_signal[i] = alpha * signal[i] + (1 - alpha) * weighted_signal[i - 1]
return weighted_signal
# 时间加权常数,单位为秒
# Slow: 1s, Fast: 125ms, Impulse: 35ms
tau = 0.125
# 频率计权
def ABC_weighting(curve='A'):
if curve not in 'ABC':
raise ValueError('Curve type not understood')
z = [0, 0]
p = [-2 * math.pi * 20.598997057568145,
-2 * math.pi * 20.598997057568145,
-2 * math.pi * 12194.21714799801,
-2 * math.pi * 12194.21714799801]
k = 1
if curve == 'A':
p.append(-2 * math.pi * 107.65264864304628)
p.append(-2 * math.pi * 737.8622307362899)
z.append(0)
z.append(0)
elif curve == 'B':
p.append(-2 * math.pi * 10 ** 2.2) # exact
z.append(0)
b, a = zpk2tf(z, p, k)
k /= abs(freqs(b, a, [2 * math.pi * 1000])[1][0])
return np.array(z), np.array(p), k
def A_weighting(fs, output='ba'):
z, p, k = ABC_weighting('A')
z_d, p_d, k_d = bilinear_zpk(z, p, k, fs)
if output == 'zpk':
return z_d, p_d, k_d
elif output in {'ba', 'tf'}:
return zpk2tf(z_d, p_d, k_d)
elif output == 'sos':
return zpk2sos(z_d, p_d, k_d)
else:
raise ValueError(f"'{output}' is not a valid output form.")
# 使用 A 频率计权处理信号
b, a = A_weighting(fs, output='ba')
A_signal = lfilter(b, a, pressure_signal)
# 信号平方
A_signal_squared = A_signal ** 2
# 时间加权处理信号
_weighted_signal = time_weighting(A_signal_squared, fs, tau)
# 转换为dB
_weighted_signal_dB = 10 * np.log10(_weighted_signal / P_ref)
# # 可视化加权信号
plt.figure(figsize=(10, 6))
plt.plot(t, _weighted_signal_dB, label='Weighted Signal')
plt.xlabel('Time (s)')
plt.ylabel('Pressure (dB(A))')
plt.legend()
plt.title('L_AF Pressure Signals')
plt.grid(True)
plt.show()
# # 可视化加权信号
plt.figure(figsize=(10, 6))
plt.plot(t, pressure_signal, label='Weighted Signal')
plt.xlabel('Time (s)')
plt.ylabel('Pressure (pa)')
plt.legend()
plt.title('L_AF Pressure Signals')
plt.grid(True)
plt.show()
A计权滤波实现参考本号文章:《A计权的原理、算法及python 实现(有码)》
感谢小伙伴的关注,后续将陆续为大家奉上数字数据处理原理及代码实现等内容....
如果您有场景、有算法...
何不考虑用我们的低成本采集器和免费软件进行二次开发创造自己的产品?
关注本号,联系我们,创造就开始了
千元机,8通道IEPE、ICP
可以快速部署算法及业务的免费软件平台