I have an acoustic signal from a Ffowcs Williams Hawkings CFD analysis and would like to convert it to the frequency domain and see the SPL and OASPL. I know I need to use fft() but I am unsure why my plot differs so much from the one from FLUENT.

time = data{1};
p = data{2};

Fs = 1/(time(end)-time(end-1));
L = 0.887*Fs+1;
N = 2^9;
Y = fft(p,N);

for k = 1:N/2
   f(k) = (k-1)*Fs/N;
  • $\begingroup$ Welcome to SE.DSP. Could you please expand some acronyms, and provide more details for those not familiar with FLUENT, CFD, etc. $\endgroup$ – Laurent Duval Mar 2 '16 at 19:23

Steps are quite simple:

  • Apply the window to your signal
  • Calculate the DFT of your signal
  • Scale the DFT properly (half of spectra and coherent gain)
  • Convert the result to decibel scale with a reference acoustic pressure of $2\cdot 10^{-5} \mathrm{Pa}$

I don't have MATLAB with me right now so here is some Python code. Considering 1 kHz sinusoid with amplitude of 1 Pascal (94 dB SPL). Resulting plot:

enter image description here

Source code:

#!/usr/bin/env python

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':


    # Generate 94 dB (Peak) SPL sinusoid
    fs = 48000
    duration = 1
    npts = int(fs*duration)
    t = np.linspace(0, 2*np.pi, npts)
    f = 500
    ref = 2.0e-5 # 20 uPa
    amp = 1      # 1 Pa -> 94 dB SPL
    s = amp * np.sin((f*duration) * t)

    # Window signal
    win = np.hamming(npts)
    signal = s * 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.
    # To obtain RMS values, divide by sqrt(2)
    sp_mag = np.abs(sp) * 2 / np.sum(win)

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

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

    plt.plot(t, s)
    plt.xlabel('Time [s]')
    plt.ylabel('Acoustic pressure [Pa]')

    plt.semilogx(freq, sp_db)
    plt.xlim( (0, fs/2)  )
    plt.ylim( (0, 100))
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('SPL [dB]')

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.