Should the sample rate be incorporated when calculating sound exposure levels?

Question

A sound exposure level is defined as the integral of squared pressures, over a defined time-period and frequency range, and it is reported with units of dB re: 1 µPa$^2$s (for water-borne acoustics). My understanding is that a sound exposure level is supposed to be an integral, not a sum, with the major difference being that an integral incorporates the distance between each value (e.g. the "dt"), which is dependent on the sample rate.

In my discussions with a few scientists on this topic, I’m finding inconsistent methods in the way that sample rates are (or are not) included in the calculation of sound exposure levels. This is particularly concerning because sound exposure levels are often used in management contexts, like the US National Marine Fisheries Marine Mammal Acoustic Technical Guidance and Southall et al 2019 Marine Mammal Noise Exposure Criteria.

Here's a snippet of R-code which demonstrates how I think the SEL should be calculated

#Calculate a SEL over 100 ms window, from 20 Hz - 1000 Hz 

library(tuneR)
library(seewave)

clip<-readWave(wavfile)
WavClip<-clip@left - mean(clip@left) #Account for DC offset
fs<[email protected]

#Calibration
cal=177 #full system calibration  
cal = 10^(cal/20) #convert from dB to linear
Nbit <- clip@bit
WavBit <- WavClip/(2^(Nbit-1))
WavCal<-WavBit*cal

#Limit frequence range with bandpass filter
WavFilt<- bwfilter(WavCal, f = fs, n =4, from = 20, to=1000, output = "Sample",bandpass = TRUE)

#Limit time range 
winLen_sec<-0.1*fs 
Wav100ms<-WavFilt[1:winLen_sec]

#Integrate squared pressures, divide by sample rate, convert to dB
SEL<-10*log10(sum(Wav100ms^2)/fs) 

Answer 1

My understanding is that a sound exposure level is supposed to be an integral, not a sum, with the major difference being that an integral incorporates the distance between each value (e.g. the "dt"), which is dependent on the sample rate.

Generally speaking An integral (over a continuous variable $t$) can be approximated as a sum (over a discrete variable $k=t/{∆t}$) (see the exact mathematical relashionship) if the infinitesimal variable $dt$ is approximated by a small quantity $∆t$:

$$\int_{t=t_1}^{t=t_2}f(t)dt \simeq\sum_{k=k_1}^{k=k_2}f(t=k\Delta t)\Delta t $$

with $$\Delta t = 1/f_s$$ $$t_1=k_1\Delta t$$ $$t_2=k_2\Delta t$$

In your case $f(k)= $ Wav100ms$^2$ (calibrated squared pressure), $t_1= 0$, $t_2= 0.1$ sec and the shortest time step $∆t$ available is the sample period $1/f_s$. So the factor $1/f_s$ is necessary, otherwise the integral estimate is definitely wrong. Then the SEL estimation $10\log_{10}(\sum$Wav100ms$^2/f_s$) in your R code is correct.

NB: From SEL = 10*log10(sum(*Wav100ms*^2)/fs), it looks like that changing $f_s$ will change SEL. Strictly speaking this is right (the approximation of an integral to a sum of a finite number of elements is dependent of the step length, here the sample period $1/f_s$. However, if your sample rate $f_s$ is standard (i.e. high enough to record the frequencies you are dealing with; Nyquist-Shannon theorem), then the effect of changing $f_s$ the sound exposure level is negligible compared to the variance of the sound exposure level over time and other uncertainties.

Answer 2

It is correct that the sound exposure level (SEL) is an integral and it value is typically increasing with time. The classical approximation is

SEL=10*log10(sum(wav^2) * dt)

As

dt = 1/fs

you snippet is correct

SEL=10*log10(sum(wav^2) /fs)

However, the text is misleading and should read

#Integrate squared pressure and convert to dB

or

#sum squared pressures, divide by sample rate, convert to dB

Answer 3

A common way to calculate the exposure level (LE, also known as sound exposure level, https://www.iso.org/obp/ui/#iso:std:iso:18405:ed-1:v1:en paragraph 3.2.1.5) would be using the "trapz" method (available as in-build function in MatLab/Octave etc.). Indeed this is a common integration method for solving integrals numerically.

Note that in the two examples below the time step is used, and the LE is therefore sample-rate independent.

Something like this in MatLab:

sig_length = 1000; # 1000 samples
sig = rand(sig_length,1); # noise assumed to be in Pascal
fs = 500; # sample rate
PaReference = 1e-6; # reference pressure
time = (0:sig_length-1)/fs; # time in seconds
dT = 1/fs; # time step
LE = 10*log10(trapz(dT,sig.^2)/PaReference^2)
# trapz can also take a vector instead of "dT" if timestep is not constant

# or manually
cumEnergy = zeros(sig_length,1); # preallocate memory
sig_sqr = sig.^2; # squared values for energy calcs
cumEnergy(1) = sig_sqr(1); #  set first sample
for i = 2:sig_length
    # if timestep is not constant use dT = time(i)-time(i-1);
    dPasqr = sig_sqr(i)-sig_sqr(i-1);
    rect = sig_sqr(i-1)*dT; # area of square from last sample to current sample
    tri = (dPasqr*dT)/2; # area of triangle given by dT and dPasqr
    cumEnergy(i) = cumEnergy(i-1) + rect + tri; 
endfor
# above loop can probably be done with vectors only (no loop needed) if you're good

LE_manual = 10*log10(cumEnergy(end)/Parefence^2) 

In the above we could also have "cumEnergy" be a single value that we update, but with this approach we can calculate the LE over time, with LE=10*log10(cumEnergy./PaReference^2)