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

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() ————————————————