audio - SPL values from FFT of microphone signal

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

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 frequency domain, you will obtain 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 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 number of samples $N$. For different types of windows this factor that is equal to sum of all window samples.

Additionally we are only interested in one half of a spectrum, therefore amplitude of all samples must be multiplied by $2$ to compensate the loss of energy, except of DC component which appears only once!.

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

plt.close('all')

fs = 5e5
duration = 1
npts = int(fs*duration)
t = np.arange(npts, dtype=float)/fs
f = 1000

ref = 0.004
amp = ref * np.sqrt(2)
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/ref)

# Window signal
win = np.hamming(npts)
signal = signal * win


sp = np.fft.fft(signal)
freq = np.fft.fftfreq(npts, 1.0/fs)

# Scale the magnitude of FFT by window energy and factor of 2, 
# because we are using half of FFT.
sp_mag = np.abs(sp) * 2 / np.sum(win)

# To obtain RMS values, divide by sqrt(2)
sp_rms = sp_mag / np.sqrt(2)


# Shift both vectors to have DC at center
freq   = np.fft.fftshift(freq)
sp_rms = np.fft.fftshift(sp_rms)

# Convert to decibel scale
sp_db = 20 * np.log10( sp_rms/ref ) + 94

plt.semilogx(freq, sp_db)
plt.xlim( (0, fs/2)  )

plt.ylim( (-100, 100))
plt.grid('on')

# Compare the outputs
print dbspl, sp_db.max()

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