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

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)

• Welcome to SE.DSP. Could you please expand some acronyms, and provide more details for those not familiar with FLUENT, CFD, etc. – 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:

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()