How do I calculate OASPL from SPL ?

I have successfully calculated SPl from my acoustic data but I'm struggling to calculate OASPL now. I understand that there is two ways to do this, one is to calculate OASPL from Prms and the other by integrating over SPL data. I have done both but both give me extremely high OASPL values (150 dB). What am I doing wrong?

1. OASPL from Prms equation

# Constants

P_ref = 20e-6 # Reference pressure in Pascals

# Microphone data

microphones = ['M2', 'M3', 'M5', 'M6']

angles = np.array([63.43, 68.96, 72.64, 90])

Cd = np.array([1.00, 1.25, 1.50, 10.03])

# Function to calculate pressure from raw voltages

def pressure_calc(measured_voltage, sensitivity):

return measured_voltage / sensitivity

# Function to calculate oaspl

def calculate_oaspl(pressure_data):

P_rms = np.sqrt(np.mean(pressure_data**2, axis=0)) # RMS pressure for each microphone

OASPL = 20 * np.log10(P_rms / P_ref) # OASPL in dB

return OASPL

# Function to correct oaspl for directivity using Cd

def correct_oaspl(OASPL, Cd):

return OASPL + 20 * np.log10(Cd)

# Main function to process data from a directory

def DATAPRO(directory):

# Experimental sensitivity for each microphone

experimental_sensitivity = [3.78, 3.44, 3.66, 3.67] # mV ?

OASPL_uncorrected_list = []

OASPL_corrected_list = []

for root, _, files in os.walk(directory):

for filename in files:

if filename.endswith('.csv'):

print(f'Processing: {filename}')

raw_data = pd.read_csv(os.path.join(root, filename))

# renaming columns and drop non-working microphones

columns = ['mic1', 'mic2', 'mic3', 'mic4', 'mic5', 'mic6', 'mic7', 'mic8', 'mic9', 'mic10']

raw_data.columns = columns

raw_data = raw_data.drop(columns=['mic1', 'mic4', 'mic7', 'mic8', 'mic9', 'mic10'])

pressure_data = []

for i, col in enumerate(raw_data.columns):

pressure_d = pressure_calc(raw_data[col], (0.001 * experimental_sensitivity[i]))

pressure_data.append(pressure_d)

pressure_df = pd.DataFrame(pressure_data).transpose()

OASPL_uncorrected = calculate_oaspl(pressure_df.to_numpy())

OASPL_corrected = correct_oaspl(OASPL_uncorrected, Cd[:len(OASPL_uncorrected)]) #corrected oaspl for directivity

OASPL_uncorrected_list.append(OASPL_uncorrected)

OASPL_corrected_list.append(OASPL_corrected)

plot_oaspl_vs_directivity(angles, OASPL_uncorrected, OASPL_corrected, title=f'OASPL vs. Directivity ({filename})')

return OASPL_uncorrected_list, OASPL_corrected_list

# plot oaspl vs directivity

def plot_oaspl_vs_directivity(angles, OASPL_uncorrected, OASPL_corrected, title):

plt.figure()

plt.plot(angles[:len(OASPL_uncorrected)], OASPL_uncorrected, 'o-', label='Uncorrected OASPL')

plt.plot(angles[:len(OASPL_corrected)], OASPL_corrected, 'o-', label='Corrected OASPL')

plt.xlabel('Angle (degrees)')

plt.ylabel('OASPL (dB)')

plt.title(title)

plt.grid(True)

plt.legend()

plt.show()

2. OASPL from integrating over SPL data

import numpy as np

import pandas as pd

import os

from scipy.fftpack import fft, fftfreq

# Constants

P_ref = 20e-6 # Reference pressure in Pascals (20 µPa)

sampling_rate = 139999 # Sampling rate in Hz

# Function to calculate pressure from raw voltages

def pressure_calc(measured_voltage, sensitivity):

"""

Convert raw voltage data to pressure in Pascals.

:param measured_voltage: Raw voltage data (in volts).

:param sensitivity: Microphone sensitivity (in volts per Pascal).

:return: Pressure data in Pascals.

"""

return measured_voltage / sensitivity

# Function to apply Hanning window

def hanning(mic_data):

window = np.hanning(len(mic_data))

return mic_data * window

# Function to perform FFT

def to_fft(data):

columns = data.columns

data_fftd = data.copy()

for col in columns:

X = fft(data[col].to_numpy())

data_fftd[col] = X

data_fftd = data_fftd.iloc[:len(data_fftd)//2]

return data_fftd

# Function to calculate SPL from FFT

def spl(fft_sig, N):

spl_array = []

p_ref = 20e-6

fft_mag = (2/N) * np.abs(fft_sig)

for col in fft_sig.columns:

spl_d = 20 * np.log10(fft_mag[col] / p_ref)

spl_array.append(spl_d)

return spl_array

# Function to get pressure values

def get_pressure_vals(raw_data):

count = 0

pressure_data = []

experimental_sensitivity = [3.78, 3.44, 3.66, 3.67] # Sensitivity for M2, M3, M5, M6 (in mV/Pa)

for col in raw_data.columns:

# Convert sensitivity from mV/Pa to V/Pa

sensitivity_V_per_Pa = experimental_sensitivity[count] * 1e-3

# Calculate pressure (raw data is already in volts)

pressure_d = pressure_calc(raw_data[col], sensitivity_V_per_Pa)

pressure_data.append(pressure_d)

count += 1

pressure_df = pd.DataFrame(pressure_data)

pressure_df = pressure_df.transpose()

return pressure_df

# Function to apply Fourier transform and Hanning window

def apply_fourier_hanning(pressure_df, N):

freqs = fftfreq(N, 1/sampling_rate)[:N//2]

windowed_pressure = pressure_df.apply(hanning)

fftd_pressure = to_fft(windowed_pressure)

spl_fft_data = spl(fftd_pressure, N)

spl_fft_df = pd.DataFrame(spl_fft_data)

spl_fft_df = spl_fft_df.transpose()

return spl_fft_df, freqs

# Function to calculate OASPL from SPL

def calculate_oaspl(spl_df):

"""

Calculate OASPL by integrating SPL over frequency.

:param spl_df: DataFrame containing SPL values for each microphone.

:return: OASPL in dB for each microphone.

"""

spl_values = spl_df.values

oaspl = 10 * np.log10(np.sum(10 ** (spl_values / 10), axis=0)) # Sum over frequency bins

return oaspl

# Main function to process data

def process_directory(directory):

spl_results = []

freqs_results = []

for root, _, files in os.walk(directory):

for filename in files:

if filename.endswith('.csv'):

print(f'Processing: {filename}')

# Load raw data from CSV file

raw_data = pd.read_csv(f'{directory}/{filename}')

# Define column names and drop non-working microphones

columns = ['mic1', 'mic2', 'mic3', 'mic4', 'mic5', 'mic6', 'mic7', 'mic8', 'mic9', 'mic10']

raw_data.columns = columns

raw_data = raw_data.drop(columns=['mic1', 'mic4', 'mic7', 'mic8', 'mic9', 'mic10']) # Keep only working microphones

# Process data

N = len(raw_data)

pressure_df = get_pressure_vals(raw_data)

spl_fft_df, freqs = apply_fourier_hanning(pressure_df, N)

print(f'Finished: {filename}')

# Store results

spl_results.append(spl_fft_df)

freqs_results.append(freqs)

else:

continue

return spl_results, freqs_results

# Example usage

directory = 'acoustic_data/2_0bar'

spl_results, freqs_results = process_directory(directory)

# Calculate OASPL for each file

oaspl_results = [calculate_oaspl(spl_df) for spl_df in spl_results]

# Microphone names

microphones = ['M2', 'M3', 'M5', 'M6']

# Output OASPL for each microphone

for i, oaspl in enumerate(oaspl_results):

print(f'OASPL for File {i+1}:')

for mic, oaspl_value in zip(microphones, oaspl):

print(f'{mic}: {oaspl_value:.2f} dB')

print('-' * 20)

[VBBV]

can you upload the MATLAB code instead of python ?

Answer 1

[VBBV]

pressure_data.append(pressure_d)

this line seems to add the pressure data cumulatively.

you need to use pressure values separately to calculate OASPL

[Nelly]

I see, what would be the correct way to fix this, what function would i need to use here