SPL graph over time

[silviaurraa]

I'm using the ASA object with a test folder. To do this, I selected 3 bands (63Hz, 125Hz, 1kHz) and 'rms' as the method_name to use asa.evolution. The problem is I end up getting negative values for the 1kHz frequency, and I don't think it makes sense. I don't know if it's a matter of scale or some other factor I'm not taking into account.

import pyhydrophone as pyhy
import pypam
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
mpl.rcParams['text.usetex'] = False

# Soundtrap
model = "SoundTrap 300 HF"
name = "H1_6338"
serial_number = 6338
soundtrap = pyhy.soundtrap.SoundTrap(
    name=name, model=model,serial_number=serial_number)

# SURVEY PARAMETERS
nfft = 4096 #puntos FFT, resolucion espectral
# para afinar más las bandas nfft = 8192
binsize = 60.0 #en segundos, fragmento temporal para el analisis

#bandas de interes
band_63 = [56.23, 70.79]
band_125 = [11.36, 140.36]
band_1000 = [889.05, 1122.46]
band_list = [band_63, band_125, band_1000]

features = ["rms", "peak", "sel"] #metricas acusticas que se calculan
#rms: energia total, peak: pico maximo, SEL: sound sxposure level

third_octaves = None
dc_subtract = True

include_dirs = False
zipped_files = False

# Run the different parameters
asa = pypam.ASA(
    hydrophone=soundtrap,
    folder_path="//pnetapp/PAM/TFM_PRUEBAS/2.despliegue/6338/",
    binsize=binsize,
    nfft=nfft,
    fft_overlap = 0.5,
    timezone="UTC",
    include_dirs=include_dirs,
    zipped=zipped_files,
    dc_subtract=dc_subtract,
)


#visualizacion RMS () evolucion temporal
rms_evolution = asa.evolution('rms', band_list=band_list)
#RMS (Root Mean Square) Nivel medio de presión acústica durante un intervalo de tiempo
#cuantificar ruido de naturaleza continua (fuentes como trafico marino o ruido ambiente)
#devuelve un array

#visualizar evolution 
tiempo = rms_evolution['datetime'].values
rms = rms_evolution['rms'].values

band_labels = ['63 Hz', '125 Hz', '1 kHz']

plt.figure(figsize=(12,6))

for i,label in enumerate(band_labels):
    plt.plot(tiempo, rms[:,i], label=label)
    
plt.xlabel('Tiempo')
plt.ylabel('SPL (dB)')
plt.legend()
plt.show()

[cparcerisas]

Hi Silvia,

Can you attach the obtained plot?
And could you tell me the version of pypam you're using?
Thanks!

[silviaurraa]

The version is 0.3.2

Thanks!

[cparcerisas]

Hi Silvia, thanks for the extra information!
You have just found a bug in the code, which I will need to fix! Thanks for that, I will work on it so it does not happen in future versions.
The bug only affects when you use multiple bands for the analysis. You can have a work-around at the moment like this:

import pyhydrophone as pyhy
import pypam
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
mpl.rcParams['text.usetex'] = False

# Soundtrap
model = "SoundTrap 300 HF"
name = "H1_6338"
serial_number = 6338
soundtrap = pyhy.soundtrap.SoundTrap(
    name=name, model=model, serial_number=serial_number)

# SURVEY PARAMETERS
nfft = 4096  # puntos FFT, resolucion espectral
# para afinar más las bandas nfft = 8192
binsize = 60.0  # en segundos, fragmento temporal para el analisis

# bandas de interes
band_63 = [56.23, 70.79]
band_125 = [11.36, 140.36]
band_1000 = [889.05, 1122.46]
band_list = [band_63, band_125, band_1000]


dc_subtract = True
include_dirs = False
zipped_files = False

# Run the different parameters
asa = pypam.ASA(
    hydrophone=soundtrap,
    folder_path="//pnetapp/PAM/TFM_PRUEBAS/2.despliegue/6338/",
    binsize=binsize,
    nfft=nfft,
    fft_overlap=0.5,
    timezone="UTC",
    include_dirs=include_dirs,
    zipped=zipped_files,
    dc_subtract=dc_subtract,
)

# visualizacion RMS () evolucion temporal
# calcular un rms por banda
rms_evolution1 = asa.evolution('rms', band_list=[band_63])
rms_evolution2 = asa.evolution('rms', band_list=[band_125])
rms_evolution3 = asa.evolution('rms', band_list=[band_1000])

# Crear diccionary
rms = {'63 Hz': rms_evolution1, '125 Hz': rms_evolution2, '1 kHz': rms_evolution3}

# visualizar evolution
fig, ax = plt.subplots(figsize=(12, 6))

for label, rms_i in rms.items():
    rms_i.rms.plot(ax=ax, label=label)

plt.xlabel('Tiempo')
plt.ylabel('SPL (dB)')
plt.legend()
plt.show()

However, I see you're basically trying to obtain decidecade bands and plot them - there are easier methods in pypam for that.
I give you here two examples (depends if you're interested in using the fft method or the filtering method"

FFT method:

milli_psd = asa.hybrid_millidecade_bands(db=True, method="density", band=[0, 4000], percentiles=None)
decidecade = pypam.utils.hmb_to_decidecade(ds=milli_psd, data_var='millidecade_bands', freq_coord='frequency_bins', fs=milli_psd.fs + 1)
selected_freqs = [63.1 , 125.89, 1000]

fig, ax = plt.subplots()
for freq in selected_freqs:
    rms_freq = decidecade.where(decidecade.frequency_bins==freq, drop=True)
    rms_freq.millidecade_bands.plot(x='datetime', ax=ax, label=f'{freq} Hz')
plt.xlabel('Tiempo')
plt.ylabel('SPL (dB)')
plt.title('Evolucion SPL')
plt.legend()
plt.show()

Filtering method

third_oct = asa.evolution_freq_dom('third_octaves_levels')
selected_freqs = [62.5, 125, 1000] # los numeros tienen que ser estos exactamente, ya que es el centro de la banda de frequencia
fig, ax = plt.subplots()
for freq in selected_freqs:
    rms_freq = third_oct.where(third_oct.frequency==freq, drop=True)
    rms_freq.oct3.plot(x='datetime', ax=ax, label=f'{freq} Hz')
plt.xlabel('Tiempo')
plt.ylabel('SPL (dB)')
plt.title('Evolucion SPL')
plt.legend()
plt.show()

Ya me diras si funciona! Mientras tanto voy a solucionar el bug!

Saludos,

Clea

[silviaurraa]

Hi Clea!

Thank you su much for the options. The first one and the last one work! The second one gives an error that I can't fix:

IndexError: index 224 is out of bounds for axis 1 with size 182

Saludos,

Silvia

[cparcerisas]

Can you tell me the original sampling frequency of your data? And which nfft and band you're using?

[silviaurraa]

milli_psd.fs = 192000, nfft = 4096, band = [0, 4000]

[cparcerisas]

If you're computing hybrid milidecade bands, nfft has to be double of the maximum of the selected band, can you try that?
nfft = 8000

[silviaurraa]

Yes, it works. With the first option I have obtain graph with some spikes every time the audio starts/ends. Could it be due to calibration pitch? With calibration = 5.0 in the ASA instance, it seems to correct itself. I'm not sure if I'm right.

Thank you so much.

[cparcerisas]

Yes indeed, if you have a sountrap with the calibration tones activated, the calibration tones will spike at the beginning of each file. Do you know if they are activated?

[silviaurraa]

Yes they are activated...

[cparcerisas]

Then there you have your answer, the spikes are the calibration tones :)
And you can indeed ignore them by setting calibration=5, which then ignores the first 5 seconds of each file

[silviaurraa]

Okay perfect! thank you so much again! :)