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(???)
  • 1
    $\begingroup$ Why not use a test signal and calibrate? $\endgroup$ – Bjorn Roche Mar 4 '14 at 17:12
  • $\begingroup$ How would I do this, take the value of the fft at the tone frequency, and use this as 94 dB reference for all? $\endgroup$ – hackyday Mar 4 '14 at 17:42
  • $\begingroup$ Yes. You could try a few frequencies to validate the calibration. $\endgroup$ – Bjorn Roche Mar 5 '14 at 22:25
  • $\begingroup$ Would I not want to use only a frequency for which I receive 4mV peak value as a reference? All other values would be a an unknown dB SPL level (due to speaker roll-off). $\endgroup$ – hackyday Mar 7 '14 at 0:46
  • $\begingroup$ it sounds like you need to calibrate your whole system, including your adc. I would use a tone generator with a known (or measured) output in place of the mic, then validate the mic's sensitivity with a real sound and a dB meter. $\endgroup$ – Bjorn Roche Mar 7 '14 at 4:17

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.xlabel("Frequency [Hz]")
plt.ylabel("SPL [dB]")


enter image description here

Difference in levels: -1.4210854715202004e-14 dB

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

| improve this answer | |
  • $\begingroup$ As I understand your code, this look's like the SPL is only the highest peak of the spectrum, so I don't have to build the average over the whole frequency-range? And what about Leq, is it the average of SPL over a time? $\endgroup$ – user20157 Mar 21 '16 at 23:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.