6
$\begingroup$

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(???)
$\endgroup$
  • 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
1
$\begingroup$

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!

$\endgroup$
  • $\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.