SPL values from FFT of microphone signal

Question

I am trying to get a frequency response vs dB SPL plot for recorded sounds. I have a microphone with a flat response, a sensitivity of 4mV/Pa, and 1 Pa = 94dB.

To get the dB SPL, I have the following formula:

$dB SPL = 94 + 20 \cdot \log_{10}(AI/0.004)$

where AI in the input voltage from the microphone, and assuming no gain.

This seems to work to get the overall dB SPL level for the signal, using peak voltage as AI.

My question is, how can I apply this to an FFT of the input signal, to get the intensity level per frequency?

I saw this post, and am not sure how to adapt this to input voltage.

Here is some python code to demonstrate what I have thus far, my goal is to apply this to ultrasonic signals that contain many different frequency components.

import numpy as np

fs = 5e5
duration = 1
npts = fs*duration
t = np.arange(npts, dtype=float)/fs
f = 1000
amp = 1.414*0.004   
signal = amp * np.sin((f*duration) * np.linspace(0, 2*np.pi, npts))

rms = np.sqrt(np.mean(signal**2))
dbspl = 94 + 20*np.log10(rms/0.004)

sp = np.fft.fft(signal)
freq = np.arange(npts)/(npts/fs)

sp_db = 94 + 20*np.log10(???)

Answer 1

Let's assume that signal you are analysing is sinusoid with amplitude $A$:

$x=a\sin{2\pi f_{0} t}$

Its RMS value of amplitude is then: $\dfrac{a}{\sqrt{2}}$ as you noticed in your code. Before performing DFT it is good practice to window signal as it decreases leakage, etc.

After transforming this signal into the frequency domain, you will obtain a two-sided spectrum. Then we have to take the magnitude of these complex values. Probably you would notice that all of them are very high. That's because energy is not normalised. Normalisation can be done via dividing values of magnitude by energy of your window. As people tend to take signal "as it is" (which is equal to applying rectangular window) the spectrum is then divided by a number of samples $N$. For different types of windows this factor that is equal to the sum of all window samples.

Additionally, we are only interested in one half of a spectrum, therefore the amplitude of all samples must be multiplied by $2$ to compensate the loss of energy, except DC component (which appears only once) and Nyquist bin (if it's present).

Lastly, as you want to analyse RMS of the signal, you must understand the following. All values of single-sided spectrum are of height $\dfrac{A_{i}^{2}}{2}$, which is equal to $ \left( \dfrac{A_{i}}{\sqrt{2}}\right)^{2}$, where $\dfrac{A_{i}}{\sqrt{2}}$ is the RMS magnitude amplitude for the i-th frequency component. Translating this into Python code you obtain:

import numpy as np
import matplotlib.pyplot as plt

fs = 5e5       # Sampling frequency [Hz]
duration = 10  # Duration of the test signal [s]

N = int(fs * duration)
t = np.arange(N, dtype=float) / fs

f = 1000  # Frequency of the test signal [Hz]
sensitivity = 0.004  # Sensitivity of the microphone [mV/Pa]
p_ref = 2e-5 # Reference acoustic pressure [Pa]
amp = sensitivity * np.sqrt(2)  # Amplitude of sinusoidal signal with RMS of 4 mV (94 dB SPL)
signal = amp * np.sin(2 * np.pi * f * duration * t)  # Signal [V]

# Calculate the level from time domain signal
rms_time = np.sqrt(np.mean(signal**2))
db_time = 20 * np.log10(rms_time / sensitivity / p_ref)

# Apply window to the signal
win = np.hamming(N)
signal = signal * win

# Get the spectrum and shift it so that DC is in the middle
spectrum = np.fft.fftshift( np.fft.fft(signal) )
freq = np.fft.fftshift( np.fft.fftfreq(N, 1 / fs) )

# Take only the positive frequencies
spectrum = spectrum[N//2:]
freq = freq[N//2:]

# Since we just removed the energy in negative frequencies, account for that
spectrum *= 2
# If there is even number of samples, do not normalize the Nyquist bin
if N % 2 == 0:
    spectrum[-1] /= 2

# Scale the magnitude of FFT by window energy
spectrum_mag = np.abs(spectrum) / np.sum(win)

# To obtain RMS values, divide by sqrt(2)
spectrum_rms = spectrum_mag / np.sqrt(2)
# Do not scale the DC component
spectrum_rms[0] *= np.sqrt(2) / 2

# Convert to decibel scale
spectrum_db = 20 * np.log10(spectrum_rms / sensitivity / p_ref)

# Compare the outputs
print("Difference in levels: {} dB".format(db_time - spectrum_db.max()))

plt.semilogx(freq, spectrum_db)
plt.xlim( (1, fs/2)  )
plt.ylim( (0, 120) )
plt.grid('on')
plt.xlabel("Frequency [Hz]")
plt.ylabel("SPL [dB]")

plt.show()

Difference in levels: -1.4210854715202004e-14 dB

As you can see, this is producing consistent results for both time and frequency domain. Good luck!