Hybrid Millidecade bands

This script is an example to compute hybrid millidecade bands, save them as netCDF files and use them to plot

Create the hydrophone object

import pyhydrophone as pyhy
import pypam

# Soundtrap
model = 'ST300HF'
name = 'SoundTrap'
serial_number = 67416073
soundtrap = pyhy.soundtrap.SoundTrap(name=name, model=model, sensitivity=-172.8, serial_number=serial_number)

Set the study parameters

# First, decide band to study. The top frequency should not be higher than the nyquist frequency (sampling rate/2)
band = [0, 4000]

# Then, set the nfft to double the sampling rate. If you want to pass None to your band, that is also an option, but
# then you need to know the sampling frequency to choose the nfft.
nfft = band[1] * 2  # or nfft = 8000

# Set the band to 1 minute
binsize = 60.0

Declare the Acoustic Survey

include_dirs = False
zipped_files = False
dc_subtract = True
asa = pypam.ASA(hydrophone=soundtrap, folder_path='./../tests/test_data', binsize=binsize, nfft=nfft,
                timezone='UTC', include_dirs=include_dirs, zipped=zipped_files, dc_subtract=dc_subtract)

# Compute the hybrid millidecade bands
# You can choose 'density' or 'spectrum' as a method
milli_psd = asa.hybrid_millidecade_bands(db=True, method='density',
                                         band=band,
                                         percentiles=None)
print(milli_psd['millidecade_bands'])

  0%|          | 0/3 [00:00<?, ?it/s]../tests/test_data/67416073.210610033655.wav

  0%|          | 0/5 [00:00<?, ?it/s]
 20%|##        | 1/5 [00:00<00:03,  1.09it/s]
100%|##########| 5/5 [00:01<00:00,  6.19it/s]


 33%|###3      | 1/3 [00:01<00:02,  1.07s/it]../tests/test_data/67416073.210610034155.wav

  0%|          | 0/5 [00:00<?, ?it/s]
 80%|########  | 4/5 [00:00<00:00, 36.37it/s]


 67%|######6   | 2/3 [00:01<00:00,  1.84it/s]../tests/test_data/67416073.210610034655.wav

  0%|          | 0/5 [00:00<?, ?it/s]
 80%|########  | 4/5 [00:00<00:00, 36.17it/s]


100%|##########| 3/3 [00:01<00:00,  2.66it/s]
100%|##########| 3/3 [00:01<00:00,  2.11it/s]
<xarray.DataArray 'millidecade_bands' (id: 18, frequency_bins: 1400)>
array([[107.12148692, 106.35207795,  99.1527314 , ..., -44.37420226,
        -50.02876127, -53.10621923],
       [107.36452089, 106.40475592,  97.8538451 , ..., -44.85179299,
        -49.91661407, -52.66584546],
       [107.36884989, 106.82555761,  97.49484291, ..., -45.10001757,
        -50.91589651, -52.62115783],
       ...,
       [107.39769473, 106.33787223,  93.79044492, ..., -44.63627933,
        -50.29965203, -52.77922535],
       [108.31286521, 106.93937232,  94.29078787, ..., -44.07367095,
        -49.58126669, -52.59687466],
       [ 56.61729663,  53.61831716,  12.93788447, ..., -26.27747371,
        -26.28720262, -28.90494788]])
Coordinates:
  * id               (id) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
    datetime         (id) datetime64[ns] 2021-06-10T03:36:55 ... 2021-06-10T0...
    start_sample     (id) int64 0 480000 960000 ... 1440000 1920000 2400000
    end_sample       (id) int64 480000 960000 1440000 ... 2400000 2880000
    file_path        (id) <U44 '../tests/test_data/67416073.210610033655.wav'...
    _start_frame     (id) int64 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  * frequency_bins   (frequency_bins) float64 0.0 1.0 2.0 ... 3.995e+03 4e+03
    lower_frequency  (frequency_bins) float64 -0.5 0.5 ... 3.99e+03 3.999e+03
    upper_frequency  (frequency_bins) float64 0.5 1.5 2.5 ... 3.999e+03 4e+03
Attributes:
    units:          uPa^2/Hz
    standard_name:  PSD

Now, with the obtained hybrid millidecade bands, make some plots We will first load some pre-computed data

import xarray

milli_psd_day = xarray.open_dataset('./../tests/test_data/test_day.nc')
milli_psd_day = milli_psd_day.where(milli_psd_day.frequency_bins > 10, drop=True)

# Plot the spectrum mean with the standard deviation
pypam.plots.plot_spectrum_mean(milli_psd_day, data_var='millidecade_bands')

# Plot the SPD with percentiles
percentiles = [10, 50, 90]
spd = pypam.utils.compute_spd(milli_psd_day, data_var='millidecade_bands', percentiles=percentiles)
pypam.plots.plot_spd(spd)

# Plot a summary
pypam.plots.plot_summary_dataset(milli_psd_day, data_var='millidecade_bands',
                                 percentiles=percentiles, freq_coord='frequency_bins')

• 
• 
•