5
$\begingroup$

Here is a 10 seconds-long 440hz sine wave normalized at $0\textrm{ dBFS}$.

When computing the STFT (with the code below) of this audio file, I noticed that max(abs(STFT)) is around 248.33. (more generally, it seems to be approximately fftsize/4 for this particular file).

Even after doing magdB = 20 * math.log10(abs(STFT)) to have decibel values, we get a max of 47.9 dB.

But this $47.9\textrm{ dB}$ doesn't mean anything! We would like to have ~ $0\textrm{ dB}$ instead because this audio only contains 1 sine wave component (1 "harmonic") at $0\textrm{ dB}$ volume.

Audio tools display it properly, here is the spectrum display of this sine wave:

enter image description here

  • Question: how to have an absolute, canonical $\rm dB$ values from a FFT, that makes that a pure $0\textrm{ dBFS}$ sinewave has a peak of $0\textrm{ dB}$ in the spectrum display ?

  • More generally, is there a canonical way to go from FFT values to $\textrm{dB}$ in order to display a spectrum analysis?

NOTE: I know that $\textrm{dB}$ are for ratio between 2 things, so all I mentioned here is referred to $\textrm{dBFS}$ (full-scale)

import scipy, numpy as np
import scipy.io.wavfile as wavfile

def stft(x, fftsize=1024, overlap=4): 
    hop = fftsize / overlap
    w = scipy.hanning(fftsize)
    return np.array([np.fft.rfft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])
Improve this question
$\endgroup$
10
  • $\begingroup$ Ahm, why don't you just rescale your signal, let's say by $FFT_{dB} = 20\log_{10}{\frac{FFT_{orig}}{\textrm{max}\,FFT_{orig}}}$? Yes, dB is a relative measure between two values, one of them being a fixed reference. I assume that you do not have a specially crafted measurement setup that is capable of measuring absolute values. If so, you wold have to use the maximum possible value as the normalization constant in the aforementioned formula. Of course, the FFT-values in this formula are absolute values. $\endgroup$ – M529 Jul 11 '16 at 12:00
  • $\begingroup$ I'm computing and displaying spectrum in realtime, so I don't know the max of future FFT frames. Moreover such a division wouldn't work if signal is constantly zero, etc. $\endgroup$ – Basj Jul 11 '16 at 12:06
  • $\begingroup$ Certainly your spectrum's y-axis will vary if you have signals of different volume, yes. Since you can compute 47.9 dB (in your example), you can therefore also compute the maximum of the absolut FFT signal. If you encounter a higher value after that, update your maximum value. Basically your y-scale will adapt in this situation. Also, there probably always is noise, so division by zero is an unlikely issue - but could also be prevented by a if-statement, if necessary. $\endgroup$ – M529 Jul 11 '16 at 12:17
  • $\begingroup$ According to docs.scipy.org/doc/scipy-0.14.0/reference/generated/…, you get an array of integer samples unless the file itself contains floats. So you need to get the maximum integer of the data type, and take that as your 0dB reference. $\endgroup$ – Sebastian Reichelt Jul 11 '16 at 15:49
  • 2
    $\begingroup$ That's probably due to other problems, but I can only guess: 1. Window normalization, i.e. the window function is defined to compensate how the energy is spread across multiple bins. 2. Incorrect way of determining the amplitude. (It seems the "abs" function works correctly for complex values in Python, so that seems unlikely.) 3. Dropping negative frequencies, resulting in a reduction by 6dB. $\endgroup$ – Sebastian Reichelt Jul 12 '16 at 20:09
8
$\begingroup$

Definitely you will have to calibrate your system. You need to know what is the relationship between dBFS (Decibel Full-Scale) and dB scale you want to measure.

In case of digital microphones, you will find sensitivity given in dBFS. This corresponds to dBFS level, given 94 dB SPL (Sound Pressure Level). For example this microphone for input $94 \;\mathrm{dB SPL}$ will produce signal at $-26 \;\mathrm{dBFS}$. Therefore the calibration factor for spectrum will be equal to $94+26=120\;\mathrm{dB}$.

Secondly, keep in mind scaling of the spectrum while doing windowing and DFT itself. More specifically, given amplitude spectrum (abs(sp)):

  • Total energy of the signal is spread over frequencies below and above Nyquist frequency. Naturally we are interested only in half of the spectrum. That is why, it important to multiply it by 2 and ignore everything above the Nyquist frequency.

  • Whenever doing windowing, it is necessary to compensate for loss of energy due to multiplication by that window. This is defined as division by sum of window samples (sum(win)). In case of rectangular window (or now window), it is as simple as division by N, where N is the DFT length.

Here is some example source code in Python. I am sure that you can take it from here.

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wf

plt.close('all')


def dbfft(x, fs, win=None, ref=32768):
    """
    Calculate spectrum in dB scale
    Args:
        x: input signal
        fs: sampling frequency
        win: vector containing window samples (same length as x).
             If not provided, then rectangular window is used by default.
        ref: reference value used for dBFS scale. 32768 for int16 and 1 for float

    Returns:
        freq: frequency vector
        s_db: spectrum in dB scale
    """

    N = len(x)  # Length of input sequence

    if win is None:
        win = np.ones(1, N)
    if len(x) != len(win):
            raise ValueError('Signal and window must be of the same length')
    x = x * win

    # Calculate real FFT and frequency vector
    sp = np.fft.rfft(x)
    freq = np.arange((N / 2) + 1) / (float(N) / fs)

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

    # Convert to dBFS
    s_dbfs = 20 * np.log10(s_mag/ref)

    return freq, s_dbfs


def main():
    # Load the file
    fs, signal = wf.read('Sine_440hz_0dB_10seconds_44.1khz_16bit_mono.wav')

    # Take slice
    N = 8192
    win = np.hamming(N)
    freq, s_dbfs = dbfft(signal[0:N], fs, win)

    # Scale from dBFS to dB
    K = 120
    s_db = s_dbfs + K

    plt.plot(freq, s_db)
    plt.grid(True)
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')
    plt.show()

if __name__ == "__main__":
    main()

Lastly, there is no better source on this type of spectrum scaling than brilliant publication by G. Heinzel et al. Just keep in mind that if you want to proper RMS scaling, then full-scale sinosuidal signal will be $-3 \;\mathrm{dBFS}$.

You might also find my previous answer being helpful.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Hi, quick question about the relationship between dbSPL and dbFS. If 94 dBSPL gives a response of -26dBFS, shouldn't the scaling from dbFS to dB be a multiplication by (94/-26)? so instead of dB = dbSPL + dbFS it would be dB = dbFS * (dbSPL/dbFS) $\endgroup$ – Richard Bamford Mar 18 '20 at 13:59
  • 1
    $\begingroup$ @RichardBamford It would be a multiplication in linear domain, since decibels are logarithmic, then it becomes addition. $\endgroup$ – jojek♦ Mar 18 '20 at 14:10
  • $\begingroup$ Aha! I see, thanks. $\endgroup$ – Richard Bamford Mar 18 '20 at 15:38
  • $\begingroup$ i have 2 more questions, why is it 20 * log() instead of the usual 10 * log()? additionally, does this equation work in the time domain too? $\endgroup$ – Richard Bamford Mar 19 '20 at 20:28
  • $\begingroup$ @RichardBamford: It's because the equation for SPL is: 10*log10(p^2/p0^2). When you remove the squares from logarithm you end up with 20*log10(p/p0). And yes - you can do exactly the same from time domain (Parseval's identity). For details have a look here. $\endgroup$ – jojek♦ Mar 20 '20 at 9:46
2
$\begingroup$

Input a known magnitude pure sinewave unmodulated reference signal whose frequency is integer periodic in the FFT width, and use the resulting magnitude as a value to rescale all other input. (e.g. calibrate your system).

For other measurements, note that after rescaling you may also have to interpolate (parabolic or Sync kernel) all other magnitude peaks that appear to be of a frequency that is between result bin center frequencies (e.g. not purely integer periodic in aperture).

Improve this answer
$\endgroup$
2
  • $\begingroup$ +1 That's what I thinking but not sure of this "empirical" method. It will work, but isn't there a more canonical way? Isn't there a more standard way to define 0dB in frequency domain / FFT ? $\endgroup$ – Basj Jul 11 '16 at 13:58
  • $\begingroup$ It depends on your particular FFT implementation and window. Different ones require different scale factors. $\endgroup$ – hotpaw2 Jul 11 '16 at 14:02
2
$\begingroup$

Some details to Jojek's post. Shared here for future reference.

The key is to take the window in consideration and taking a reference ref which corresponds to 0dBFS.

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wf

fs, signal = wf.read('Sine_440hz_0dB_10seconds_44.1khz_16bit_mono.wav')  # Load the file
ref = 32768  # 0 dBFS is 32678 with an int16 signal

N = 8192
win = np.hamming(N)                                                       
x = signal[0:N] * win                             # Take a slice and multiply by a window
sp = np.fft.rfft(x)                               # Calculate real FFT
s_mag = np.abs(sp) * 2 / np.sum(win)              # Scale the magnitude of FFT by window and factor of 2,
                                                  # because we are using half of FFT spectrum
s_dbfs = 20 * np.log10(s_mag / ref)               # Convert to dBFS
freq = np.arange((N / 2) + 1) / (float(N) / fs)   # Frequency axis
plt.plot(freq, s_dbfs)
plt.grid(True)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.show()

Then we can get 0dB peak in FFT when working with a 0dB sinusoid :

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ This change was not accepted by me, since it is misleading. You created a "fake" reference value, which is completely different from the one, widely used, in dBFS calculations. Additionally mentioning "tricks", without trying to understand them is not a good idea. $\endgroup$ – jojek♦ Jul 13 '16 at 16:47
  • 1
    $\begingroup$ Why is this a fake reference value? It's the same than yours excepted the 120dB thing. If you're speaking about reference 32768 that I removed in the FFT dB computation, this is totally normal, because I work with signal in float in [-1, 1], by dividing by 32768 in the 6th line. Thus, this division by 32768 later would make no sense. Do you see what I mean? $\endgroup$ – Basj Jul 13 '16 at 16:50
  • $\begingroup$ In any decibel scale, the reference value defines everything. For example in case of acoustic pressure it is $20 \mu \mathrm{Pa}$, for sound intensity level it is $10^{-12} \mathrm{W/m^2}$. In decibel full-scale it is a maximum absolute value: $32768$ for int16 wave format and 1 for float format. My goal was to put variables where they belong, not to shuffle divisions and multiplications in a confusing manner- exact reason why I reverted your change. $\endgroup$ – jojek♦ Jul 13 '16 at 16:53
  • $\begingroup$ @jojek I corrected the code, it's closer to your solution. (Just shared it as an extra answer because it avoids the hardware-dependent 94 dB SPL thing in yours, which is correct, but that adds another complexity when trying to understand the concept). $\endgroup$ – Basj Oct 11 '18 at 20:47

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.