I have a pressure signal and want to do SPL analysis on it

Question

Signal

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;
end
spl=20*log10(abs(Y(1:length(Y)/2)));
semilogx(f,spl,'k','Linewidth',2)

Answer 1

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:

Source code:

#!/usr/bin/env python

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

if __name__ == '__main__':

    plt.close('all')

    # 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.subplot(2,1,1)
    plt.plot(t, s)
    plt.xlabel('Time [s]')
    plt.ylabel('Acoustic pressure [Pa]')
    plt.grid(True)

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