Sound pressure level FFT (Python)
Introduction
What are the characteristics (natural frequency, order sound, etc.) of the measured time signal in the frequency domain? I will share the program code of ** FFT (Fast Fourier Transform) **, which is one of the signal processing methods used when investigating.
Preface
・ The basic concept will be omitted because it was explained by another person. ・ This is my first post, so I would appreciate it if you could watch over with warm eyes. ・ PyCharm 2020.3.2 (Community Edition)
Purpose
In the articles I have read so far, I could not find a graph with frequency [Hz] on the horizontal axis and sound pressure level [dB] or noise level [dBA] on the vertical axis (only insufficient research). So, this time, I will try frequency analysis by [dB].
FFT condition of data to be analyzed
-The analysis target is the tapping sound data when the outer wall tile is vibrated once per second for 10 seconds in a white noise environment.
・ Sampling rate fs: 8192 Hz
・ FFT block size block_size: 8192 points
・ Overlap overlap: 0.75
・ Frequency resolution df: 1.0 Hz
・ Time required for one FFT: 1.0 sec
・ Number of tapping sound measurement microphone channels: 1 ch
imports
import numpy as np #Numerical library
import matplotlib.pyplot as plt #Drawing library
import japanize_matplotlib #Japanese localization of graph axis labels and legends
import scipy.integrate as si #Integral
import scipy.signal as ss
import sympy as sp #Necessary to calculate how many FFTs must be done after setting the FFT conditions
FFT code
FFT.py
def fft(data, fs, block_size, output_ch, overlap):
max_frequency = fs / 2.56 #Upper limit analysis frequency
line_number = block_size / 2.56 #Number of lines
df = max_frequency / line_number #Frequency resolution
dt = 1 / fs
block_size_range = np.arange(0, block_size, 1)
nf = np.int((fs / 2) + 1) #Nyquist frequency
nf_line = (block_size_range / (dt * block_size))[:nf] #Frequency axis
impact_time_data = data[:, 0] #Time waveform of hammer with force sensor attached (impulse input)
output_time_data = data[:, 1] #Response point microphone (microphone dedicated to striking sound measurement)
time_line = np.arange(0, data.shape[0] / fs, dt) #Time axis
# def time_graph():
# plt.rcParams['font.size'] = 30
# plt.figure(figsize=(10, 10))
# plt.plot(time_line, output_time_data, label='Response point microphone', linewidth=1, color='C1')
# plt.ylabel('Sound pressure[Pa]')
# plt.xlabel('time[s]')
# plt.xlim(0, 1)
# plt.ylim(-5, 5)
# plt.yticks(np.arange(-5, 6, 1))
# plt.grid(True)
# plt.legend(loc="upper right", borderaxespad=0, edgecolor='black', framealpha=1, fancybox=False)
# plt.tight_layout()
#
# plt.figure(figsize=(10, 10))
# plt.plot(time_line, impact_time_data, label='Vibration hammer', linewidth=1, color='C2')
# plt.xlabel('time[s]')
# plt.ylabel('Power[N]')
# plt.xlim(0, 1)
# plt.ylim(0, 50)
# plt.grid(True)
# plt.legend(loc="upper right", borderaxespad=0, edgecolor='black', framealpha=1, fancybox=False)
# plt.tight_layout()
# time_graph()
"iteration calculation"
n = sp.Symbol('n')
right = 0 + (n - 1) * block_size * (1 - overlap)
left = data.shape[0] - block_size
zero = right - left
iteration = np.int(sp.solve(zero)[0])
"FFT starting point"
start_point = [
np.int(0 + (i - 1) * block_size * (1 - overlap))
for i in range(1, iteration + 1)
]
output_use = [
np.reshape(output_time_data[start_point[i]:start_point[i] + block_size], (block_size, output_ch))
for i in range(iteration)
]
"Simple integration"
def hanning_integral():
global hanning_simpson
hanning_simpson = si.simps(np.hanning(block_size), dx=1 / fs)
hanning_correction = 1 / hanning_simpson
print('hanning_The window loss rate is{0}'.format(hanning_simpson))
print('So relatively{0}Need to double'.format(hanning_correction))
return hanning_correction
hanning_window_correction = hanning_integral()
"Exact solution"
def exact_solution(x):
hanning_formula = 0.5 - 0.5 * sp.cos(2 * np.pi * x)
return hanning_formula
global hanning_integration
hanning_integration = si.quad(exact_solution, 0, 1)
print('Exact solution{0} '.format(1 / hanning_integration[0]))
# FFT
output_fft = [
[
((np.fft.fft(np.hanning(block_size) * output_use[i][:, j])[:nf]) /
block_size) * (2 * 1 / hanning_integration[0])
for i in range(iteration)
]
for j in range(output_ch)
]
'--------------------Matrix size conversion (FFT repeat count x number of response point channels)--------------------'
output_reshape1 = [
np.reshape(output_fft[j], (iteration, nf))
for j in range(output_ch)
]
output_reshape2 = [
[
output_reshape1[j][:, k]
for j in range(output_ch)
]
for k in range(nf)
]
output_reshape3 = [
np.reshape(output_reshape2[k], (output_ch, iteration)).T
for k in range(nf)
]
'--------------------Sound pressure level--------------------'
output_dB = [
20 * np.log10(np.sqrt(np.mean(np.power(np.abs(output_reshape3[k]), 2), axis=0)) / 2e-5)
for k in range(nf)
]
output_dB = np.reshape(output_dB, (nf, output_ch))
# def fft_graph():
# plt.figure(figsize=(10, 10))
# plt.plot(nf_line, output_dB, color='C0', linewidth=3, label="Response point")
# plt.xlim(0, 2000)
# plt.xticks(np.arange(0, 2200, 200))
# plt.axvspan(1020, 1100, color='green', alpha=0.3)
# plt.ylim(0, 80)
# plt.yticks(np.arange(0, 90, 10))
# plt.xlabel('frequency[Hz]')
# plt.ylabel("Sound pressure level[dB]")
# plt.grid(True)
# plt.legend(loc="upper right", borderaxespad=0, edgecolor='black', framealpha=1, fancybox=False)
# plt.tight_layout()
# pass
#
# fft_graph()
return output_reshape3, nf, nf_line, output_dB, iteration, fs, block_size, output_ch, overlap
Explanation for each part
iteration calculation (FFT iteration count)
If you set the overlap, the FFT will be performed multiple times in a row, but you need to calculate how many times you have to FFT.
The number of sound pressure signal data measured this time is 81920 (sampling rate x measurement time).
75% overlap is the same as 25% off the FFT 1 block from the first FFT start, right? (Figure below)
So, can the starting point for repeated FFT be treated as an arithmetic progression with the first term being 0 (Because it is Python, the Index starts from 0) and the intersection being 2048 (8192 points x 0.25)?
The number of FFT repetitions can be found by finding n at the 73728th time.
How much amplitude is lost in the raw data in the window function?
Here, the calculation of the correction amount is explained.
By applying a hanning window, both ends of the FFT 1 block will always be 0, and leakage errors can be prevented. In response to this merit, the signal was changed by force, so it is necessary to correct it and return it. Here, let's calculate the area when the window function is not applied and when the Hanning window is applied.
The area is the area below the line of the Hanning window. So the area is compared to when you can't write anything in the window
hanning_simpson = si.simps(np.hanning(block_size), dx=1 / fs)
This will be smaller. Therefore, it is necessary to correct and return (in short, multiply by the reciprocal) as it becomes smaller.
Function call
FFT.py
noise_impact_output_ffts, nf, nf_line, noise_impact_output_ffts_dB, ffts_iteration, fs, block_size, output_ch, overlap = fft(data=noise_impact, fs=8192, block_size=8192, output_ch=1, overlap=0.75)
Results obtained
FFT of white noise + tapping sound data gave such a result. Originally, the resonance point of the tapping sound should appear in the figure, but it is considered that the tapping sound component was buried in the noise when the vibration was applied in the white noise environment this time. (It would have been better to analyze the data excited in a quiet environment ...) Also, it can be seen that the signal is a constant signal of about 50 dB at all frequencies without paying attention to the resonance frequency. Therefore, it matches the characteristics of white noise.
** This time, FFT was performed using all the data, but since the tapping sound is instantaneous, the FFT should be performed after trimming 1 second around the moment when the force is input. (I will write it in another article.) **
How was it?
Thank you for reading this far. This time, I tried frequency analysis in dB notation considering human auditory characteristics. I hope this article helps someone.
Corrections
Tag name change (Python3➡Python)