Skip to content

Commit

Permalink
pylint
Browse files Browse the repository at this point in the history
  • Loading branch information
@jmcvey3
jmcvey3 committed Sep 10, 2024
1 parent c294b05 commit 7252d9f
Show file tree
Hide file tree
Showing 6 changed files with 238 additions and 166 deletions.
17 changes: 9 additions & 8 deletions examples/acoustics_example.ipynb
View file

Large diffs are not rendered by default.

9 changes: 8 additions & 1 deletion mhkit/acoustics/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,9 @@
from .base import *
"""
The `acoustics` package of the MHKiT (Marine and Hydrokinetic Toolkit) library
provides tools and functionalities for analyzing and visualizing passive
acoustic monitoring data deployed in water bodies. This package reads in raw
wav files and conducts basic acoustics analysis and visualization.
"""

from mhkit.acoustics import io, graphics
from .base import *
103 changes: 67 additions & 36 deletions mhkit/acoustics/base.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
"""
This file includes the main passive acoustics analysis functions. They
are designed to function on top of one another, starting from reading
in wav files from the io submodule.
"""

import warnings
import numpy as np
import xarray as xr
Expand All @@ -6,6 +12,21 @@
from mhkit.dolfyn.time import epoch2dt64, dt642epoch


def _fmax_warning(fn, fmax):
"""
Check that max frequency limit isn't greater than Nyquist frequency
"""

if fmax > fn:
warnings.warn(
"`fmax` = {fmax} is greater than the Nyquist frequency. Setting"
"fmax = {fn}"
)
fmax = fn

return fmax


def minimum_frequency(water_depth, c=1500, c_seabed=1700):
"""
Estimate the shallow water cutoff frequency based on the speed of
Expand All @@ -20,7 +41,7 @@ def minimum_frequency(water_depth, c=1500, c_seabed=1700):
return f_min


def sound_pressure_spectral_density(P, fs, window=1):
def sound_pressure_spectral_density(pressure, fs, window=1):
"""
Calculates the mean square sound pressure spectral density from audio
samples split into FFTs with a specified window_size in seconds and
Expand All @@ -29,7 +50,7 @@ def sound_pressure_spectral_density(P, fs, window=1):
Parameters
----------
P: xarray.DataArray (time)
pressure: xarray.DataArray (time)
Sound pressure in [Pa] or volts [V]
fs: int
Data collection sampling rate [Hz]
Expand All @@ -49,8 +70,8 @@ def sound_pressure_spectral_density(P, fs, window=1):
binner = VelBinner(n_bin=win, fs=fs, n_fft=win)
# Always 50% overlap if numbers reshape perfectly
# Mean square sound pressure
psd = binner.power_spectral_density(P, freq_units="Hz")
samples = binner.reshape(P.values) - binner.mean(P.values)[:, None]
psd = binner.power_spectral_density(pressure, freq_units="Hz")
samples = binner.reshape(pressure.values) - binner.mean(pressure.values)[:, None]
# Power in time domain
t_power = np.sum(samples**2, axis=1) / win
# Power in frequency domain
Expand All @@ -62,7 +83,7 @@ def sound_pressure_spectral_density(P, fs, window=1):
psd_adj,
coords={"time": psd_adj["time"], "freq": psd_adj["freq"]},
attrs={
"units": P.units + "^2/Hz",
"units": pressure.units + "^2/Hz",
"long_name": "Mean Square Sound Pressure Spectral Density",
"fs": fs,
"window": str(window) + "s",
Expand Down Expand Up @@ -104,8 +125,8 @@ def apply_calibration(spsd, sensitivity_curve, fill_value):
calibration = calibration.fillna(fill_value)

# Subtract from sound pressure spectral density
Sf_ratio = 10 ** (calibration / 10) # V^2/uPa^2
spsd /= Sf_ratio # uPa^2/Hz
sensitivity_ratio = 10 ** (calibration / 10) # V^2/uPa^2
spsd /= sensitivity_ratio # uPa^2/Hz
spsd /= 1e12 # Pa^2/Hz
spsd.attrs["units"] = "Pa^2/Hz"

Expand All @@ -130,10 +151,10 @@ def sound_pressure_spectral_density_level(spsd):
"""

# Reference value of sound pressure
P2_ref = 1e-12 # Pa^2 to 1 uPa^2
reference = 1e-12 # Pa^2 to 1 uPa^2

# Sound pressure spectral density level from mean square values
lpf = 10 * np.log10(spsd.values / P2_ref)
lpf = 10 * np.log10(spsd.values / reference)

out = xr.DataArray(
lpf,
Expand Down Expand Up @@ -177,13 +198,9 @@ def band_average(
indexed by time and frequency
"""

# Check fmax
fn = spsdl["freq"].max().values
if fmax > fn:
warnings.warn(
"`fmax` = {fmax} is greater than the Nyquist frequency. Setting"
"fmax = {fn}"
)
fmax = fn
fmax = _fmax_warning(fn, fmax)

bandwidth = 2 ** (1 / octave)
half_bandwidth = 2 ** (1 / (octave * 2))
Expand Down Expand Up @@ -278,24 +295,20 @@ def sound_pressure_level(spsd, fmin=10, fmax=100000):
Sound pressure level [dB re 1 uPa] indexed by time
"""

# Check fmax
fn = spsd.attrs["fs"] // 2
if fmax > fn:
warnings.warn(
"`fmax` = {fmax} is greater than the Nyquist frequency. Setting"
"fmax = {fn}"
)
fmax = fn
fmax = _fmax_warning(fn, fmax)

# Reference value of sound pressure
P2_ref = 1e-12 # Pa^2, = 1 uPa^2
reference = 1e-12 # Pa^2, = 1 uPa^2

# Mean square sound pressure in a specified frequency band from mean square values
P2 = np.trapz(
pressure_squared = np.trapz(
spsd.sel(freq=slice(fmin, fmax)), spsd["freq"].sel(freq=slice(fmin, fmax))
)

# Mean square sound pressure level
mspl = 10 * np.log10(P2 / P2_ref)
mspl = 10 * np.log10(pressure_squared / reference)

out = xr.DataArray(
mspl,
Expand All @@ -312,16 +325,34 @@ def sound_pressure_level(spsd, fmin=10, fmax=100000):


def _band_sound_pressure_level(spsd, bandwidth, half_bandwidth, fmin, fmax):
"""
Calculates band-averaged sound pressure levels
Parameters
----------
spsd: xarray.DataArray (time, freq)
Mean square sound pressure spectral density.
bandwidth: int
Bandwidth to average over
half_bandwidth: int
Half-bandwidth, used to set upper and lower bandwidth limits
fmin: int
Lower frequency band limit (lower limit of the hydrophone). Default: 10 Hz
fmax: int
Upper frequency band limit (Nyquist frequency). Default: 100000 Hz
Returns
-------
out: xarray.DataArray (time, freq_bins)
Sound pressure level [dB re 1 uPa] indexed by time and frequency of specified bandwidth
"""

# Check fmax
fn = spsd.attrs["fs"] // 2
if fmax > fn:
warnings.warn(
"`fmax` = {fmax} is greater than the Nyquist frequency. Setting"
"fmax = {fn}"
)
fmax = fn
fmax = _fmax_warning(fn, fmax)

# Reference value of sound pressure
P2_ref = 1e-12 # Pa^2, = 1 uPa^2
reference = 1e-12 # Pa^2, = 1 uPa^2

center_freq = 10 ** np.arange(
np.log10(fmin),
Expand All @@ -333,20 +364,20 @@ def _band_sound_pressure_level(spsd, bandwidth, half_bandwidth, fmin, fmax):
octave_bins = np.append(lower_limit, upper_limit[-1])

# Manual trapezoidal rule to get Pa^2
P2 = xr.DataArray(
pressure_squared = xr.DataArray(
coords={"time": spsd["time"], "freq_bins": center_freq},
dims=["time", "freq_bins"],
)
for i, key in enumerate(center_freq):
band_min = octave_bins[i]
band_max = octave_bins[i + 1]
P2.loc[dict(freq_bins=key)] = np.trapz(
pressure_squared.loc[{"freq_bins": key}] = np.trapz(
spsd.sel(freq=slice(band_min, band_max)),
spsd["freq"].sel(freq=slice(band_min, band_max)),
)

# Mean square sound pressure level in dB rel 1 uPa
mspl = 10 * np.log10(P2 / P2_ref)
mspl = 10 * np.log10(pressure_squared / reference)

return mspl

Expand All @@ -358,7 +389,7 @@ def third_octave_sound_pressure_level(spsd, fmin=10, fmax=100000):
Parameters
----------
psd: xarray.DataArray (time, freq)
spsd: xarray.DataArray (time, freq)
Mean square sound pressure spectral density.
fmin: int
Lower frequency band limit (lower limit of the hydrophone). Default: 10 Hz
Expand Down Expand Up @@ -391,7 +422,7 @@ def decidecade_sound_pressure_level(spsd, fmin=10, fmax=100000):
Parameters
----------
psd: xarray.DataArray (time, freq)
spsd: xarray.DataArray (time, freq)
Mean square sound pressure spectral density.
fmin: int
Lower frequency band limit (lower limit of the hydrophone). Default: 10 Hz
Expand Down
36 changes: 16 additions & 20 deletions mhkit/acoustics/graphics.py
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,14 @@
import warnings
"""
This submodule includes the main passive acoustics plotting functions. All
functions allow passthrough of matplotlib functionality and commands
to make them fully customizable.
"""

import matplotlib.pyplot as plt
from .base import _fmax_warning


def plot_spectogram(spsdl, fmin=10, fmax=100000, fig=None, ax=None, kwargs={}):
def plot_spectogram(spsdl, fmin=10, fmax=100000, fig=None, ax=None, **kwargs):
"""
Plots the spectogram of the sound pressure spectral density level.
Expand Down Expand Up @@ -32,15 +38,10 @@ def plot_spectogram(spsdl, fmin=10, fmax=100000, fig=None, ax=None, kwargs={}):
# Set dimension names
time = spsdl.dims[0]
freq = spsdl.dims[-1]

# Check fmax
fn = spsdl[freq].max()
if fmax > fn:
warnings.warn(
"`fmax` = {fmax} is greater than the Nyquist frequency. Setting"
"fmax = {fn}"
)
fmax = fn

fmax = _fmax_warning(fn, fmax)
# select frequency range
spsdl = spsdl.sel({freq: slice(fmin, fmax)})

if ax is None:
Expand All @@ -55,7 +56,7 @@ def plot_spectogram(spsdl, fmin=10, fmax=100000, fig=None, ax=None, kwargs={}):
return fig, ax


def plot_spectra(spsdl, fmin=10, fmax=100000, fig=None, ax=None, kwargs={}):
def plot_spectra(spsdl, fmin=10, fmax=100000, fig=None, ax=None, **kwargs):
"""
Plots spectral density. X axis is log-transformed.
Expand Down Expand Up @@ -84,15 +85,10 @@ def plot_spectra(spsdl, fmin=10, fmax=100000, fig=None, ax=None, kwargs={}):

# Set dimension names
freq = spsdl.dims[-1]

# Check fmax
fn = spsdl[freq].max()
if fmax > fn:
warnings.warn(
"`fmax` = {fmax} is greater than the Nyquist frequency. Setting"
"fmax = {fn}"
)
fmax = fn

fmax = _fmax_warning(fn, fmax)
# select frequency range
spsdl = spsdl.sel({freq: slice(fmin, fmax)})

if ax is None:
Expand All @@ -101,6 +97,6 @@ def plot_spectra(spsdl, fmin=10, fmax=100000, fig=None, ax=None, kwargs={}):
left=0.1, right=0.95, top=0.85, bottom=0.2, hspace=0.3, wspace=0.15
)
ax.plot(spsdl[freq], spsdl.T, **kwargs)
ax.set(xlim=(fmin, fmax), xlabel="Frequency [Hz]", ylabel="$dB\ re \ 1 \ uPa^2/Hz$")
ax.set(xlim=(fmin, fmax), xlabel="Frequency [Hz]", ylabel="$dB re 1 uPa^2/Hz$")

return fig, ax
Loading

0 comments on commit 7252d9f

Please sign in to comment.