Skip to content

Commit

Permalink
Use 'bin' instead of 'window' to avoid confusion
Browse files Browse the repository at this point in the history
  • Loading branch information
jmcvey3 committed Oct 9, 2024
1 parent 8fe42de commit 0de1b6a
Show file tree
Hide file tree
Showing 3 changed files with 28 additions and 26 deletions.
28 changes: 14 additions & 14 deletions mhkit/acoustics/analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,22 +103,22 @@ def minimum_frequency(


def sound_pressure_spectral_density(
pressure: xr.DataArray, fs: Union[int, float], window: Union[int, float] = 1
pressure: xr.DataArray, fs: Union[int, float], bin_length: Union[int, float] = 1
) -> xr.DataArray:
"""
Calculates the mean square sound pressure spectral density from audio
samples split into FFTs with a specified window_size in seconds and
at least a 50% overlap. The amplitude of the PSD is adjusted
samples split into FFTs with a specified bin length in seconds, using Hanning
windowing with 50% overlap. The amplitude of the PSD is adjusted
according to Parseval's theorem.
Parameters
----------
pressure: xarray.DataArray (time)
Sound pressure in [Pa] or voltage [V]
fs: int
fs: int or float
Data collection sampling rate [Hz]
window: string (optional)
Length of time in seconds to create FFTs. Default: 1 s.
bin_length: int or float
Length of time in seconds to create FFTs. Default: 1.
Returns
-------
Expand All @@ -131,26 +131,26 @@ def sound_pressure_spectral_density(
raise TypeError("'pressure' must be an xarray.DataArray.")
if not isinstance(fs, (int, float)):
raise TypeError("'fs' must be a numeric type (int or float).")
if not isinstance(window, (int, float)):
raise TypeError("'window' must be a numeric type (int or float).")
if not isinstance(bin_length, (int, float)):
raise TypeError("'bin_length' must be a numeric type (int or float).")

# Ensure that 'pressure' has a 'time' coordinate
if "time" not in pressure.dims:
raise ValueError("'pressure' must be indexed by 'time' dimension.")

# window length of each time series
win = window * fs
nbin = bin_length * fs

# Use dolfyn PSD
binner = VelBinner(n_bin=win, fs=fs, n_fft=win)
binner = VelBinner(n_bin=nbin, fs=fs, n_fft=nbin)
# Always 50% overlap if numbers reshape perfectly
# Mean square sound pressure
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
t_power = np.sum(samples**2, axis=1) / nbin
# Power in frequency domain
f_power = psd.sum("freq") * (fs / win)
f_power = psd.sum("freq") * (fs / nbin)
# Adjust the amplitude of PSD according to Parseval's theorem
psd_adj = psd * t_power[:, None] / f_power

Expand All @@ -161,9 +161,9 @@ def sound_pressure_spectral_density(
"units": pressure.units + "^2/Hz",
"long_name": "Mean Square Sound Pressure Spectral Density",
"fs": fs,
"window": str(window) + "s",
"nbin": str(bin_length) + " s",
"overlap": "50%",
"nfft": win,
"nfft": nbin,
},
)

Expand Down
24 changes: 13 additions & 11 deletions mhkit/tests/acoustics/test_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ class TestAnalysis(unittest.TestCase):
def setUpClass(self):
file_name = join(datadir, "6247.230204150508.wav")
P = acoustics.io.read_soundtrap(file_name, sensitivity=-177)
self.spsd = acoustics.sound_pressure_spectral_density(P, P.fs, window=1)
self.spsd = acoustics.sound_pressure_spectral_density(P, P.fs, bin_length=1)

@classmethod
def tearDownClass(self):
Expand Down Expand Up @@ -59,14 +59,14 @@ def test_averaging(self):
octave = 3
td_spsdl_mean = acoustics.band_average(td_spsdl, octave, fmin=50)

# Time average into 30 s windows
window = 30
td_spsdl_50 = acoustics.time_average(td_spsdl_mean, window, method="median")
# Time average into 30 s bins
lbin = 30
td_spsdl_50 = acoustics.time_average(td_spsdl_mean, lbin, method="median")
td_spsdl_25 = acoustics.time_average(
td_spsdl_mean, window, method={"quantile": 0.25}
td_spsdl_mean, lbin, method={"quantile": 0.25}
)
td_spsdl_75 = acoustics.time_average(
td_spsdl_mean, window, method={"quantile": 0.75}
td_spsdl_mean, lbin, method={"quantile": 0.75}
)

cc = np.array(
Expand Down Expand Up @@ -174,20 +174,22 @@ def test_sound_pressure_spectral_density(self):
data, coords=[time], dims=["time"], attrs={"units": "Pa", "fs": 100}
)

# Adjust window size to get multiple segments
# Adjust bin size to get multiple segments
fs = 100
window = 0.1 # seconds
win_samples = int(window * fs)
bin_length = 0.1 # seconds
win_samples = int(bin_length * fs)

# Run the spectral density function
spsd = acoustics.sound_pressure_spectral_density(pressure, fs=fs, window=window)
spsd = acoustics.sound_pressure_spectral_density(
pressure, fs=fs, bin_length=bin_length
)

# Assert that output is an xarray DataArray with expected dimensions
self.assertIsInstance(spsd, xr.DataArray)
self.assertIn("freq", spsd.dims)
self.assertIn("time", spsd.dims)
self.assertEqual(spsd.attrs["units"], "Pa^2/Hz")
self.assertEqual(spsd.attrs["window"], f"{window}s")
self.assertEqual(spsd.attrs["nbin"], f"{bin_length} s")

# Calculate expected number of segments
overlap = 0.0
Expand Down
2 changes: 1 addition & 1 deletion mhkit/tests/acoustics/test_io.py
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@ def test_calibration(self):
file_name = join(datadir, "6247.230204150508.wav")
td_volt = acoustics.io.read_soundtrap(file_name, sensitivity=None)
td_spsd = acoustics.sound_pressure_spectral_density(
td_volt, td_volt.fs, window=1
td_volt, td_volt.fs, bin_length=1
)

# Run calibration
Expand Down

0 comments on commit 0de1b6a

Please sign in to comment.