diff --git a/examples/acoustics_example.html b/examples/acoustics_example.html new file mode 100644 index 000000000..5bd113b14 --- /dev/null +++ b/examples/acoustics_example.html @@ -0,0 +1,622 @@ + +Analyzing Passive Acoustic Data with MHKiT

Analyzing Passive Acoustic Data with MHKiT

The following example illustrates how to read and analyze some basic parameters for passive acoustics data. Functionality to analyze .wav files recorded using hydrophones has been integrated into MHKiT to support analysis based on the IEC-TS 62600-40 standard.
The standard workflow for passive acoustics analysis is as follows:
  1. Import .wav file
  2. Calibrate data
  3. Calculate spectral density
  4. Calculate other parameters
  5. Create plots

Read in Hydrophone Measurements

All hydrophone .wav files can be read in MHKiT using a base function called read_hydrophone from the acoustics.io submodule. Because the sampling frequency is so fast, measurements are stored in the lowest memory format possible and need to be scaled and transformed to return the measurements in units of voltage or pressure.
The read_hydrophone function scales and transforms raw measurements given a few input parameters. Most parameters needed to convert the raw data are stored in the native .wav format header blocks, but two, the "peak_voltage" of the sensor's analog-to-digital converter (ADC) and file "start_time" (usually stored in the filename) are required.
Two other inputs, the hydrophone "sensitivity" and an amplifier "gain" can also be input. If a sensitivity value is provided, the function will convert voltage to pressure; otherwise the sensitivity(ies) can be provided later using a calibration curve. Gain should be provided if the instrument utilizes an amplifier gain (typically for custom hydrophone builds), which is then added to the sensitivity. For this example we will be using a different file, but
filename = "./data/acoustics/RBW_6661_20240601_053114.wav";
peak_voltage = 3; % V
sensitivity = -177; % db
gain = 0;
start_time= "2024-06-01T05:31:14";
P = read_hydrophone(filename,peak_voltage,sensitivity,gain,start_time)
P = struct with fields:
time: [01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 … ] (1×30720000 datetime) + data: [30720000×1 double] + units: 'Pa' + sensitivity: 1.4125e-09 + valid_min: -2.1238e+03 + valid_max: 2.1238e+03 + fs: 512000 + filename: "RBW_6661_20240601_053114.wav" +
“Smart” hydrophones are those where the hydrophone element, pre-amplifier board, ADC, motherboard and memory card are sold in a single package. Companies that sell these often store metadata in the .wav file header, and MHKiT has a couple of wrapper functions for these hydrophones.

Smart Hydrophone Examples

OceanSonics icListen and OceanInstruments Soundtrap are two common smart hydrophone models, with examples as follows.
For icListen datafiles, only the filename is necessary to provide to return file contents in units of pressure. The stored sensitivity calibration value can be overridden by setting the "sensitivity" input, and to return measurements in units of voltage, set sensitivity to empty and use_metadata to false.
P = read_iclisten(filename) % get sound pressure
P = struct with fields:
time: [01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 … ] (1×30720000 datetime) + data: [30720000×1 double] + units: 'Pa' + sensitivity: 1.4125e-09 + valid_min: -2.1238e+03 + valid_max: 2.1238e+03 + fs: 512000 + filename: "RBW_6661_20240601_053114.wav" + serial_number: 'icListen HF #6661' + bits_per_sample: 24 + peak_voltage: 3 + humidity: '24.0 % RH' + temperature: '8.6 deg C' + accelerometer: 'Acc(-980,-18,141)' + magnetometer: 'Mag(3603,3223,-598)' + count_at_peak_voltage: '8388608 = Max Count' + sequence_num: '2589798400000 = Seq #' +
V = read_iclisten(filename, [], false) % get raw voltage only
V = struct with fields:
time: [01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 01-Jun-2024 05:31:14 … ] (1×30720000 datetime) + data: [30720000×1 double] + units: 'V' + valid_min: -3 + valid_max: 3 + fs: 512000 + filename: "RBW_6661_20240601_053114.wav" + serial_number: 'icListen HF #6661' + bits_per_sample: 24 + peak_voltage: 3 + humidity: '24.0 % RH' + temperature: '8.6 deg C' + accelerometer: 'Acc(-980,-18,141)' + magnetometer: 'Mag(3603,3223,-598)' + count_at_peak_voltage: '8388608 = Max Count' + sequence_num: '2589798400000 = Seq #' +

SoundTrap Example

For Ocean Instruments Soundtrap datafiles, the filename and sensitivity should be provided to return the measurements in units of pressure. If the hydrophone has been calibrated, set the sensitivity to empty to return the measurements in units of voltage.
sfile = "./data/acoustics/6247.230204150508.wav";
P = read_soundtrap(sfile, -177)
P = struct with fields:
time: [04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 … ] (1×28800672 datetime) + data: [28800672×1 double] + units: 'Pa' + sensitivity: 1.4125e-09 + valid_min: -707.9458 + valid_max: 707.9458 + fs: 96000 + filename: "6247.230204150508.wav" +
V = read_soundtrap(sfile, [])

Mean Square Sound Pressure Spectral Density

After the .wav file is read in, either in units of pressure or voltage, we calculate the mean square sound pressure spectral density (SPSD) of the time-series using sound_pressure_spectral_density. This splits the time-series into windows and uses fast Fourier transforms to convert the raw measurements into the frequency domain, with units of Pa^2/Hz or V^2/Hz depending on the input. The function takes the original datafile, the hydrophone’s sampling rate (“fs”), which is stored as an attribute of the measurement time-series, and a window size (“bin_length”) in seconds as input.
The IEC-40 considers an acoustic sample to have a length of 1 second, so we’ll set the bin length as such here.
% create mean square spectral densities using 1 second bins
spsd = sound_pressure_spectral_density(V, V.fs, 1)
spsd = struct with fields:
data: [48000×599 double] + time: [04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:09 04-Feb-2023 15:05:09 04-Feb-2023 15:05:10 04-Feb-2023 15:05:10 04-Feb-2023 15:05:11 04-Feb-2023 15:05:11 04-Feb-2023 15:05:12 … ] (1×599 datetime) + freq: [48000×1 double] + name: "Mean Square Sound Pressure Spectral Density" + units: "V^2/Hz" + fs: 96000 + nbin: "1 s" + overlap: "50%" + nfft: 96000 +

Applying Calibration Curves

For conducting scientific-grade analysis, it is critical to use calibration curves to correct the SPSD calculations. Hydrophones should be calibrated (i.e., a sensitivity calibration curve should be generated for a hydrophone) every few years. The IEC-40 asks that a hydrophone be calibrated both before and after the test deployment.
A calibration curve consists of the hydrophone's sensitivity (in units of dB rel 1V^2/uPa^2) vs frequency and should be applied to the spectral density we just calculated.
The easiest way to apply a sensitivity calibration curve in MHKiT is to first copy the calibration data into a CSV file, where the left column contains the calibrated frequencies and the right column contains the sensitivity values. Here we use the function in the following section to read in a CSV file created with the column headers "Frequency" and "AnalogSensitivity".
% read in csv table
sense_table = readtable("./data/acoustics/6247_calibration.csv")
Warning: Column headers from the file were modified to make them valid MATLAB identifiers before creating variable names for the table. The original column headers are saved in the VariableDescriptions property.
Set 'VariableNamingRule' to 'preserve' to use the original column headers as table variable names.
sense_table = 40×2 table
 FrequencyAnalogSensitivity
11-223.4900
21.1830-220.8000
31.3990-218.1300
41.6550-215.4100
51.9580-212.6800
62.3160-209.9100
72.7400-207.1200
83.2410-204.2900
93.8340-201.4500
104.5350-198.5800
115.3640-195.6900
126.3450-192.7900
137.5060-189.8500
148.8790-186.9000
1510.5000-183.9300
1612.4200-180.9300
1714.7000-177.9200
1817.3800-174.9700
1920.5600-172.1600
2024.3300-169.6100
2128.7800-167.6900
2234.0400-166.5200
2340.2600-165.9600
2447.6300-165.8100
2556.3400-165.8500
2666.6500-165.9500
2778.8400-166.0500
2893.2600-166.1300
29110.3000-166.2000
30130.5000-166.2500
31154.4000-166.2800
32182.6000-166.2900
33216-166.2900
34255.5000-166.2700
35302.2000-166.2400
36357.5000-166.1700
37422.9000-166.0300
38500.3000-165.7900
39591.8000-165.4700
40700-164.8700

Apply Calibration Matrix

Now, lets simplify the table into a matrix of two columns where the first is frequency and the second holds the calibration values. Once we have the calibration data in a matrix, we can apply that to the SPSD using the apply_calibration function. Calibration curves typically do not cover the entire range of the hydrophone, so this function will linearly interpolate the missing values. A fill_value can be provided to extrapolate outside of the calibrated frequencies.
sense_curve = table2array(sense_table);
fill_Sf = sense_curve(end, 2); % use last value in curve for higher frequencies
spsd_cal = apply_calibration(spsd, sense_curve, fill_Sf)
spsd_cal = struct with fields:
data: [48000×599 double] + time: [04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:09 04-Feb-2023 15:05:09 04-Feb-2023 15:05:10 04-Feb-2023 15:05:10 04-Feb-2023 15:05:11 04-Feb-2023 15:05:11 04-Feb-2023 15:05:12 … ] (1×599 datetime) + freq: [48000×1 double] + name: "Calibrated Sound Pressure Spectral Density" + units: "Pa^2/Hz" + fs: 96000 + nbin: "1 s" + overlap: "50%" + nfft: 96000 +

Mean Square Sound Pressure Spectral Density Level

We can use the function sound_pressure_spectral_density_level to calculate the mean square sound pressure spectral density levels (SPSDLs) from the calibrated SPSD. This function converts absolute pressure into relative pressure in log-space, the traditional means with which we measure sound, in units of decibels relative to 1 uPa [dB rel 1 uPa], the standard for underwater sound.
Sidenote: Sound in air is measured in decibels relative to 20 uPa, the minimum sound pressure humans can hear. To convert between [dB rel 1 uPa] and [dB rel 20 uPa], one simply needs to subtract 26 dB from the [dB rel 1 uPa] value.
spsdl = sound_pressure_spectral_density_level(spsd_cal)
spsdl = struct with fields:
data: [48000×599 double] + time: [04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:09 04-Feb-2023 15:05:09 04-Feb-2023 15:05:10 04-Feb-2023 15:05:10 04-Feb-2023 15:05:11 04-Feb-2023 15:05:11 04-Feb-2023 15:05:12 … ] (1×599 datetime) + freq: [48000×1 double] + name: "Sound Pressure Spectral Density Level" + units: "dB re 1 uPa^2/Hz" + fs: 96000 + nbin: "1 s" + overlap: "50%" + nfft: 96000 +

Creating Spectrograms

Now that the SPSDL is calculated, we can create spectrograms, or waterfall plots, using the plot_spectrogram function in the graphics submodule. While spectrograms aren't required by the IEC-40, they are useful to do quality control so we can avoid using contaminated soundbytes in further analysis (i.e., we'd remove the boat noise shown here from further analysis of a marine energy device).
To do this, we'll give the function the minimum and maximum frequencies to plot, define the colormap, and provide min and max ranges of the colorbar. For these measurements, we're setting fmin = 10 Hz, the minimum specified by the IEC-40, and fmax = 48,000 Hz, the Nyquist frequency for these data.
Note, the IEC-40 requires a maximum frequency of 100,000 Hz, so a hydrophone capable of sampling faster than 200,000 Hz should be used for IEC testing.
% Specify spectrogram range
fmin = 10;
fmax = 48000;
 
plot_spectrogram(spsdl, ...
'fmin', fmin, ...
'fmax', fmax, ...
'cm', 'hot', ...
'use_smoothing', true, ...
'smoothing_stride', 3, ...
'plot_method', 'imagesc')
If you see something interesting in the spectrogram, the next step you should do is listen to the .wav file. This can tell you a lot about what you're looking at. If you listen to this file, you'll hear the boat cruising by around 3 minutes in. To do this, we will utilize MATLAB's built-in functions.
player = audioplayer(P.data, P.fs); % define audioplayer with data
% play(player) % uncomment to listen
% stop(player) % uncomment to stop
% export_audio('mhkit_acoustics_output.wav', P, 1) % uncomment to export audio

IEC-40 Stats

The IEC-40 requires a few aggregate statistics for characterizing the sound of marine energy devices. For the first, the IEC-40 asks for plots showing the 25%, 50%, and 75% quantiles of the SPSDL during specific marine energy device states. For current energy devices, the IEC-40 requires 10 SPSDL samples at a series of turbine states (braked, freewheel, 25% power, 50% power, 75% power, 100% power). For wave energy devices, the spec requires 30 SPSDL samples in each wave height and period bin observed.
For this example notebook we'll keep it simple and use a random set of 30 samples and collate them together. Typically, one will pick and choose which 1 second samples to collate together and analyze. Then we can find the median and quantiles of those 30 samples.
spsdl_median = spsdl;
samples = randperm(length(spsdl.data(1,:)),30);
temp = spsdl.data(:,samples);
 
spsdl_median.data = median(temp,2);
spsdlq25 = quantile(temp, 0.25, 2);
spsdlq75 = quantile(temp, 0.75, 2);
 
plot_spectra(spsdl_median,fmin,fmax)
patch([spsdl.freq' fliplr(spsdl.freq')], [spsdlq25' fliplr(spsdlq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3)
legend('Median','Quantiles')

Window Aggregating

If desired, one can group a series of measurements into blocks of time, though this isn't required by the IEC-40. In the following block, we'll take our 5 minutes of measurements, time_aggregate them into 30 second intervals, and find the median, 25% and 75% quantiles of each interval. We then plot the stats of the median parameters.
window = 30; %seconds
method = "median";
spsdl_time = time_aggregate(spsdl, window, method)
spsdl_time = struct with fields:
data: [48000×10 double] + time: [04-Feb-2023 15:05:23 04-Feb-2023 15:05:53 04-Feb-2023 15:06:23 04-Feb-2023 15:06:53 04-Feb-2023 15:07:23 04-Feb-2023 15:07:53 04-Feb-2023 15:08:23 04-Feb-2023 15:08:53 04-Feb-2023 15:09:23 04-Feb-2023 15:09:53] + freq: [48000×1 double] + units: "dB re 1 uPa^2/Hz" + name: "Sound Pressure Spectral Density Level" +
We can then use the plot_spectra function in the graphics submodule to plot the median and quantiles of the first 30s window.
spsdl_med = spsdl_time;
spsdl_med.data = median(spsdl_time.data, 2);
spsdlq25 = quantile(spsdl_time.data, 0.25, 2);
spsdlq75 = quantile(spsdl_time.data, 0.75, 2);
 
plot_spectra(spsdl_med, fmin, fmax)
patch([spsdl.freq' fliplr(spsdl.freq')], [spsdlq25' fliplr(spsdlq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3)
legend(method,'Quantiles')

Frequency Band Analysis

Frequency band analysis can be completed by grouping the data into specific frequency bands, called "band aggregating" here. In other words, instead of aggregating by the time dimension, we aggregate by the frequency dimension. The band_aggregate function operates by taking the SPSDL and grouping it based on a specified octave and octave base.
In the following section, we use the same plotting function as above, but do so by creating the decidecade frequency bands (10th octave, octave base 10 => ). We'll calculate the median and quantiles of each sample, then take the median of each of those parameters. For the IEC-40, one will calculate these parameters for each device operating state (current energy converters) or wave state matrix (wave energy converters).
octave = [10, 10]; % [octave, octave base]
method = "median";
spsdl10 = band_aggregate(spsdl,octave,fmin,fmax,method)
spsdl10 = struct with fields:
data: [38×599 double] + freq: [10 12.5893 15.8489 19.9526 25.1189 31.6228 39.8107 50.1187 63.0957 79.4328 100 125.8925 158.4893 199.5262 251.1886 316.2278 398.1072 501.1872 630.9573 794.3282 1.0000e+03 1.2589e+03 1.5849e+03 1.9953e+03 2.5119e+03 3.1623e+03 … ] (1×38 double) + time: [04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:09 04-Feb-2023 15:05:09 04-Feb-2023 15:05:10 04-Feb-2023 15:05:10 04-Feb-2023 15:05:11 04-Feb-2023 15:05:11 04-Feb-2023 15:05:12 … ] (1×599 datetime) + name: "Sound Pressure Spectral Density Level" + units: "dB re 1 uPa^2/Hz" +
spsdl10m = spsdl10;
spsdl10m.data = median(spsdl10.data, 2);
spsdlq25 = quantile(spsdl10.data, 0.25, 2);
spsdlq75 = quantile(spsdl10.data, 0.75, 2);
 
plot_spectra(spsdl10m, fmin, fmax)
patch([spsdl10.freq fliplr(spsdl10.freq)], [spsdlq25' fliplr(spsdlq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3)
legend(method,'Quantiles')
The plot above shows significant spread in sound measurements at higher frequency due to the vessel noise, but less so at lower frequency. This mirrors what we can see in the sound pressure density level figure we first plotted.
You can further build on the results with more complicated statistical methods. In the following block of code, we find the empirical quantile function (the empirical version of the cumulative distribution function, CDF) of each decidecade frequency band and plot the decidecade band centered nearest to 160 Hz.
f = 160; % target frequency
[val, idxf] = min(abs(spsdl10.freq-f));
n = numel(spsdl10.data(1,:));
p = 1.0 - (0:(n-1))' / (n + 1);
x = sort(spsdl10.data(idxf,:));
plot(x,p)
xlabel('SPSDL [dB re 1 uPa^2/Hz]'); ylabel("Exceedance Probability");
title(strcat(num2str(spsdl10.freq(idxf))," Hz"))
grid on;

Sound Pressure Level

The IEC-40 has two requirements considering calculations of sound pressure level (SPL). We'll first calculate the SPL over the full frequency range of the turbine and/or hydrophone using the function sound_pressure_level. First, however, note that the IEC-40 asks that the range be set from 10 to 100,000 Hz. The lower limit can be increased due to flow noise or low frequency signal loss due to shallow water. Shallow water cutoff frequency Low frequency sound is absorbed into the seabed in shallow water depths. We can use the function minimum_frequency to get an approximation of what our minimum frequency should be. This approximation uses the water depth, estimates of the in-water sound speed, and sea/riverbed sound speed to determine what the cutoff frequency will be. The difficult part with this approximation is figuring out the speed of sound in the bed material, which generally ranges from 1450-1800 m/s.
This function should only be used as a rough approximation and sanity check if significant attenuation is seen at various low frequencies and harmonics.
depths = [1, 5, 10, 20, 40, 80];
c = 1500; % speed of sound in water
c_seabed = 1700; % speed of sound in seabed
fmin = minimum_frequency(depths, c, c_seabed);
% print out min freq for corresponding depths
for i = 1:length(depths)
fprintf('%d m = %.1f Hz\n', depths(i), fmin(i));
end
1 m = 796.9 Hz +5 m = 159.4 Hz +10 m = 79.7 Hz +20 m = 39.8 Hz +40 m = 19.9 Hz +80 m = 10.0 Hz
Though the IEC-40 says we should default to a minimum frequency of 10 Hz, as you can see above, unless we're measuring from a depth of around 80 +/- 10 m, our minimum frequency should be higher. One can play around with the bed soundspeed to see how these change with varying bed densities/compositions.

Flow noise

Flow noise, or psuedo-sound, is the other reason to increase the minimum frequency of our SPL measurements. Flow noise is caused by one of three things: turbulence advected past the hydrophone element, turbulence caused by the hydrophone element, and the sensitivity of the hydrophone element to temperature inhomogeneities in the advected flow. Flow noise is most noticeably apparent when flow speeds increase above 0.5 m/s, seen in spectrograms as a logarithmic increase in pressure with decreasing frequency.
The particular data shown here was measured in around 8-10 m of water, and a mix of mild flow noise below 20 Hz and low frequency attenuation below ~50 Hz can be seen in the spectrogram. We'll again use the Nyquist frequency of 48,000 Hz.

Cumulative SPL

Running the code block below, we can see our cumulative SPL start out at 86 dB and then peak at 125 dB as the boat drives by. If you haven't listened to the audio track, this peak SPL of 125 dB rel 1 uPa (underwater) is equivalent to 99 dB rel 20 uPa (air). For reference, the OSHA time limit for workers experiencing 100 dB rel 20 uPa of sound is 2 hours. Vessel traffic can be quite loud and is one of the largest contributors to noise in the marine environment.
% sound pressure level
fmin = 50;
fmax = 48000;
spl = sound_pressure_level(spsd_cal, fmin, fmax);
% make figure
plot(spl.time,spl.data)
xlabel('Time'); ylabel({spl.name; spl.units}); grid on;

Decidecade Sound Pressure Levels

The last stat that IEC-40 requests are the decidecade SPLs, where the SPL is calculated for each decidecade frequency band. Note that the IEC-40 labels these as synonymous the third-octave SPLs; while for some applications this may be true, mathematically they have different definitions.
To explain, a true octave is a frequency band where the upper frequency is double (i.e., base 2) that of the lower frequency. Third octaves are often measured because mammals to have evolved to interpret sound at this bandwidth. The true one-third octave is a frequency band where the upper frequency is 2^(1/3) = 1.25992 times the lower frequency. The decidecade band referenced by the IEC-40 refers to the one-tenth octave of base 10, where the upper frequency is 10 times that of the lower frequency. Mathematically this means the decidecade band has a bandwidth of 10^(1/10) = 1.25892. So, when reporting frequency analysis, it is important to note both the octave and its bandwidth.
We can calculate the SPL in each decidecade band using the function decidecade_sound_pressure_level. This function uses the same calculation as sound_pressure_level above and runs it on each tenth octave band. It returns 1 SPL in each frequency band every timestamp. Boats are loud.
spl10 = decidecade_sound_pressure_level(spsd_cal, fmin, fmax)
spl10 = struct with fields:
data: [31×599 double] + freq: [50.0000 62.9463 79.2447 99.7631 125.5943 158.1139 199.0536 250.5936 315.4787 397.1641 500.0000 629.4627 792.4466 997.6312 1.2559e+03 1.5811e+03 1.9905e+03 2.5059e+03 3.1548e+03 3.9716e+03 5.0000e+03 6.2946e+03 7.9245e+03 … ] (1×31 double) + time: [04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:09 04-Feb-2023 15:05:09 04-Feb-2023 15:05:10 04-Feb-2023 15:05:10 04-Feb-2023 15:05:11 04-Feb-2023 15:05:11 04-Feb-2023 15:05:12 … ] (1×599 datetime) + units: 'dB re 1 uPa' + name: 'Decidecade Sound Pressure Level' +
spl10m = spl10;
spl10m.data = median(spl10.data, 2);
splq25 = quantile(spl10.data, 0.25, 2);
splq75 = quantile(spl10.data, 0.75, 2);
plot_spectra(spl10m, fmin, fmax)
patch([spl10.freq fliplr(spl10.freq)], [splq25' fliplr(splq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3)
legend('Median','Quantiles')

Third Octave Sound Pressure Level

One can also calculate the true 1/3 octave SPLs using third_octave_sound_pressure_level if desired. Note the results are quite similar to decidecade_sound_pressure_level.
spl3 = decidecade_sound_pressure_level(spsd_cal, fmin, fmax)
spl3 = struct with fields:
data: [31×599 double] + freq: [50.0000 62.9463 79.2447 99.7631 125.5943 158.1139 199.0536 250.5936 315.4787 397.1641 500.0000 629.4627 792.4466 997.6312 1.2559e+03 1.5811e+03 1.9905e+03 2.5059e+03 3.1548e+03 3.9716e+03 5.0000e+03 6.2946e+03 7.9245e+03 … ] (1×31 double) + time: [04-Feb-2023 15:05:08 04-Feb-2023 15:05:08 04-Feb-2023 15:05:09 04-Feb-2023 15:05:09 04-Feb-2023 15:05:10 04-Feb-2023 15:05:10 04-Feb-2023 15:05:11 04-Feb-2023 15:05:11 04-Feb-2023 15:05:12 … ] (1×599 datetime) + units: 'dB re 1 uPa' + name: 'Decidecade Sound Pressure Level' +
spl3m = spl3;
spl3m.data = median(spl3.data, 2);
splq25 = quantile(spl3.data, 0.25, 2);
splq75 = quantile(spl3.data, 0.75, 2);
plot_spectra(spl3m, fmin, fmax)
patch([spl3.freq fliplr(spl3.freq)], [splq25' fliplr(splq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3)
legend('Median','Quantiles')

Sound Exposure Level and Marine Mammal Weightings

While the IEC-40 does not have recommendations for sound exposure level (SEL), it is a better parameter for quantifying auditory stress/injury to organisms. Sound exposure level is a metric that quantifies sound pressure level relative to exposure duration, based on the idea that that the effects on an organism's hearing are the same for equivalent totals sound energy received - even if that energy was received at different rates. Numerically, an SPL of 190 dB re 1 uPa over a 1 second interval and an SPL 160 dB re 1 uPa for 1000 seconds both have a SEL of 190 dB re 1 uPa^2 s.
The National Marine Fisheries Service (NMFS) publishes auditory weighting functions for five groups of marine mammals. These weighting functions are designed to try to emulate the auditory sensitivity of each group to better predict auditory injury thresholds. They are mathematically equivalent to band-pass filters. A link to the latest recommendations from NMFS is https://www.fisheries.noaa.gov/national/marine-mammal-protection/marine-mammal-acoustic-technical-guidance-other-acoustic-tools
The SEL function sound_exposure_level has a few inputs. One is the pressure spectral density (PSD) in Pa^2/Hz, found from sound_pressure_spectral_density. The 'n_bin' input to that function should be the same length of time as you want to calculate SEL for. If you would like to calculate the 24 hr SEL, you need 24 hours worth of data - the best way to do this and save computer RAM would be to calculate the PSD from each 1-5 minute individual file and then concatenate the PSD from each file together.
Set "group=None" to calculate the standard SEL. Set "group" to "LF" for 'low frequency' cetaceans; "HF" for 'high frequency' cetaceans; "VHF" for 'very high frequency' cetaceans; "PW" for phocid pinnepeds; and "OW" for otariid pinnepeds.
% five minute SEL
spsd_300s = sound_pressure_spectral_density(V, V.fs, 300);
% calibrate PSD
fill_Sf = sense_curve(end, 2); % use last value in curve for higher frequencies
spsd_300s = apply_calibration(spsd_300s, sense_curve, fill_Sf);
 
group = ["LF", "HF", "VHF", "PW", "OW"];
sel = sound_exposure_level(spsd_300s, [], fmin, fmax);
fprintf('Group-Standard SEL = %.1f\n', sel.data)
Group-Standard SEL = 136.9
for g = 1:length(group)
sel = sound_exposure_level(spsd_300s, group(g), fmin, fmax);
fprintf('Group-%s SEL = %.1f\n', group(g), sel.data)
end
Group-LF SEL = 136.1 +Group-HF SEL = 133.6 +Group-VHF SEL = 128.5 +Group-PW SEL = 135.1 +Group-OW SEL = 132.7

Marine Mammal Auditory Weighting Functions

We can look at the specific marine mammal auditory weighting and noise exposure functions as well using nmfs_auditory_weighting, given a frequency vector and one of the mammal groups. It outputs the weighting function and the exposure function (the inverse of the former) in units of dB. To convert back to a unitless magnitude, use 10.^(func / 10). The exposure function shows the SEL in dB at and above which temporary or permanent hearing damage can occur to an individual in the specified group. The minimum value in the exposure function is the known or estimated injury level for a given group.
% Calculate weighting and exposure functions
[weight_func, exp_func] = nmfs_auditory_weighting(spsd_300s.freq, 'LF');
 
% Visualize weighting and exposure
figure('Position', [100, 100, 800, 400]);
 
% Weighting Function subplot
subplot(1, 2, 1);
semilogx(spsd_300s.freq, weight_func, 'LineWidth', 2, 'DisplayName', 'Weighting Function');
hold on;
yline(0, 'k--', 'LineWidth', 1.5, 'DisplayName', 'Highest Sensitivity');
xlabel('Frequency [Hz]');
ylabel('Transmission [dB]');
ylim([-50, 20]);
xlim([10, 48000]);
legend('Location', 'northeast');
grid on;
title('Auditory Weighting Function');
 
% Exposure Function subplot
subplot(1, 2, 2);
semilogx(spsd_300s.freq, exp_func, 'LineWidth', 2, 'DisplayName', 'Exposure Function');
hold on;
yline(min(exp_func), 'k--', 'LineWidth', 1.5, 'DisplayName', 'Highest Sensitivity');
 
% Fill injury region using fill()
% Note: For large datasets, this downsamples for visualization performance
if length(spsd_300s.freq) > 10000
% Downsample for plotting efficiency while preserving shape
downsample_target_size = 5000; % target size for plotting
downsample_factor = ceil(length(spsd_300s.freq) / downsample_target_size);
freq_plot = spsd_300s.freq(1:downsample_factor:end);
exp_plot = exp_func(1:downsample_factor:end);
else
freq_plot = spsd_300s.freq;
exp_plot = exp_func;
end
 
% Create filled region FROM exp_func UP TO level fill_max_y_level
fill_max_y_level = 300;
fill([freq_plot; flipud(freq_plot)], ...
[exp_plot; repmat(fill_max_y_level, size(exp_plot))], ...
'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'DisplayName', 'Injury Region');
 
xlabel('Frequency [Hz]');
ylabel('Exposure Level [dB]');
ylim([min(exp_func)-20, min(exp_func)+50]);
xlim([10, 48000]);
legend('Location', 'southeast');
grid on;
title('Auditory Exposure Function');
 
set(gcf, 'Color', 'white');
sgtitle('NMFS Marine Mammal Auditory Functions (LF Group)');
+
+ +
\ No newline at end of file diff --git a/examples/acoustics_example.m b/examples/acoustics_example.m new file mode 100644 index 000000000..1399b34ff --- /dev/null +++ b/examples/acoustics_example.m @@ -0,0 +1,495 @@ +%% Analyzing Passive Acoustic Data with MHKiT +% +% The following example illustrates how to read and analyze some basic parameters +% for passive acoustics data. Functionality to analyze .wav files recorded using +% hydrophones has been integrated into MHKiT to support analysis based on the +% IEC-TS 62600-40 standard. +% +% The standard workflow for passive acoustics analysis is as follows: +%% +% # Import .wav file +% # Calibrate data +% # Calculate spectral density +% # Calculate other parameters +% # Create plots + +%% Read in Hydrophone Measurements +% All hydrophone .wav files can be read in MHKiT using a base function called +% |read_hydrophone| from the |acoustics.io| submodule. Because the sampling frequency +% is so fast, measurements are stored in the lowest memory format possible and +% need to be scaled and transformed to return the measurements in units of voltage +% or pressure. +% +% The |read_hydrophone| function scales and transforms raw measurements given +% a few input parameters. Most parameters needed to convert the raw data are stored +% in the native .wav format header blocks, but two, the "peak_voltage" of the +% sensor's analog-to-digital converter (ADC) and file "start_time" (usually stored +% in the filename) are required. +% +% Two other inputs, the hydrophone "sensitivity" and an amplifier "gain" can +% also be input. If a sensitivity value is provided, the function will convert +% voltage to pressure; otherwise the sensitivity(ies) can be provided later using +% a calibration curve. Gain should be provided if the instrument utilizes an amplifier +% gain (typically for custom hydrophone builds), which is then added to the sensitivity. +% For this example we will be using a different file, but + +filename = "./data/acoustics/RBW_6661_20240601_053114.wav"; +peak_voltage = 3; % V +sensitivity = -177; % db +gain = 0; +start_time= "2024-06-01T05:31:14"; +P = read_hydrophone(filename,peak_voltage,sensitivity,gain,start_time) + +%% +% “Smart” hydrophones are those where the hydrophone element, pre-amplifier +% board, ADC, motherboard and memory card are sold in a single package. Companies +% that sell these often store metadata in the .wav file header, and MHKiT has +% a couple of wrapper functions for these hydrophones. +% +%% Smart Hydrophone Examples +% +% OceanSonics icListen and OceanInstruments Soundtrap are two common smart hydrophone +% models, with examples as follows. +% +% For icListen datafiles, only the filename is necessary to provide to return +% file contents in units of pressure. The stored sensitivity calibration value +% can be overridden by setting the "sensitivity" input, and to return measurements +% in units of voltage, set |sensitivity| to empty and |use_metadata| to false. + +P = read_iclisten(filename) % get sound pressure +V = read_iclisten(filename, [], false) % get raw voltage only + +%% SoundTrap Example +% +% For Ocean Instruments Soundtrap datafiles, the filename and sensitivity should +% be provided to return the measurements in units of pressure. If the hydrophone +% has been calibrated, set the sensitivity to empty to return the measurements +% in units of voltage. + +sfile = "./data/acoustics/6247.230204150508.wav"; +P = read_soundtrap(sfile, -177) +V = read_soundtrap(sfile, []) + +%% Mean Square Sound Pressure Spectral Density +% After the .wav file is read in, either in units of pressure or voltage, we +% calculate the mean square sound pressure spectral density (SPSD) of the time-series +% using |sound_pressure_spectral_density|. This splits the time-series into windows +% and uses fast Fourier transforms to convert the raw measurements into the frequency +% domain, with units of Pa^2/Hz or V^2/Hz depending on the input. The function +% takes the original datafile, the hydrophone’s sampling rate (“fs”), which is +% stored as an attribute of the measurement time-series, and a window size (“bin_length”) +% in seconds as input. +% +% The IEC-40 considers an acoustic sample to have a length of 1 second, so we’ll +% set the bin length as such here. + +% create mean square spectral densities using 1 second bins +spsd = sound_pressure_spectral_density(V, V.fs, 1) + +%% Applying Calibration Curves +% +% For conducting scientific-grade analysis, it is critical to use calibration +% curves to correct the SPSD calculations. Hydrophones should be calibrated (i.e., +% a sensitivity calibration curve should be generated for a hydrophone) every +% few years. The IEC-40 asks that a hydrophone be calibrated both before and after +% the test deployment. +% +% A calibration curve consists of the hydrophone's sensitivity (in units of +% dB rel 1V^2/uPa^2) vs frequency and should be applied to the spectral density +% we just calculated. +% +% The easiest way to apply a sensitivity calibration curve in MHKiT is to first +% copy the calibration data into a CSV file, where the left column contains the +% calibrated frequencies and the right column contains the sensitivity values. +% Here we use the function in the following section to read in a CSV file created +% with the column headers "Frequency" and "AnalogSensitivity". + +% read in csv table +sense_table = readtable("./data/acoustics/6247_calibration.csv") + +%% Apply Calibration Matrix +% +% Now, lets simplify the table into a matrix of two columns where the first +% is frequency and the second holds the calibration values. Once we have the calibration +% data in a matrix, we can apply that to the SPSD using the |apply_calibration| +% function. Calibration curves typically do not cover the entire range of the +% hydrophone, so this function will linearly interpolate the missing values. A +% fill_value can be provided to extrapolate outside of the calibrated frequencies. + +sense_curve = table2array(sense_table); +fill_Sf = sense_curve(end, 2); % use last value in curve for higher frequencies +spsd_cal = apply_calibration(spsd, sense_curve, fill_Sf) + +%% Mean Square Sound Pressure Spectral Density Level +% +% We can use the function |sound_pressure_spectral_density_level| to calculate +% the mean square sound pressure spectral density levels (SPSDLs) from the calibrated +% SPSD. This function converts absolute pressure into relative pressure in log-space, +% the traditional means with which we measure sound, in units of decibels relative +% to 1 uPa [dB rel 1 uPa], the standard for underwater sound. +% +% Sidenote: Sound in air is measured in decibels relative to 20 uPa, the minimum +% sound pressure humans can hear. To convert between [dB rel 1 uPa] and [dB rel +% 20 uPa], one simply needs to subtract 26 dB from the [dB rel 1 uPa] value. + +spsdl = sound_pressure_spectral_density_level(spsd_cal) + +%% Creating Spectrograms +% +% Now that the SPSDL is calculated, we can create spectrograms, or waterfall +% plots, using the |plot_spectrogram| function in the graphics submodule. While +% spectrograms aren't required by the IEC-40, they are useful to do quality control +% so we can avoid using contaminated soundbytes in further analysis (i.e., we'd +% remove the boat noise shown here from further analysis of a marine energy device). +% +% To do this, we'll give the function the minimum and maximum frequencies to +% plot, define the colormap, and provide min and max ranges of the colorbar. For +% these measurements, we're setting fmin = 10 Hz, the minimum specified by the +% IEC-40, and fmax = 48,000 Hz, the Nyquist frequency for these data. +% +% Note, the IEC-40 requires a maximum frequency of 100,000 Hz, so a hydrophone +% capable of sampling faster than 200,000 Hz should be used for IEC testing. + +% Specify spectrogram range +fmin = 10; +fmax = 48000; + +plot_spectrogram(spsdl, ... + 'fmin', fmin, ... + 'fmax', fmax, ... + 'cm', 'hot', ... + 'use_smoothing', true, ... + 'smoothing_stride', 3, ... + 'plot_method', 'imagesc') + +%% +% If you see something interesting in the spectrogram, the next step you should +% do is listen to the .wav file. This can tell you a lot about what you're looking +% at. If you listen to this file, you'll hear the boat cruising by around 3 minutes +% in. To do this, we will utilize MATLAB's built-in functions. + +player = audioplayer(P.data, P.fs); % define audioplayer with data +% play(player) % uncomment to listen +% stop(player) % uncomment to stop +% export_audio('mhkit_acoustics_output.wav', P, 1) % uncomment to export audio + +%% IEC-40 Stats +% +% The IEC-40 requires a few aggregate statistics for characterizing the sound +% of marine energy devices. For the first, the IEC-40 asks for plots showing the +% 25%, 50%, and 75% quantiles of the SPSDL during specific marine energy device +% states. For current energy devices, the IEC-40 requires 10 SPSDL samples at +% a series of turbine states (braked, freewheel, 25% power, 50% power, 75% power, +% 100% power). For wave energy devices, the spec requires 30 SPSDL samples in +% each wave height and period bin observed. +% +% For this example notebook we'll keep it simple and use a random set of 30 +% samples and collate them together. Typically, one will pick and choose which +% 1 second samples to collate together and analyze. Then we can find the median +% and quantiles of those 30 samples. + +spsdl_median = spsdl; +samples = randperm(length(spsdl.data(1,:)),30); +temp = spsdl.data(:,samples); + +spsdl_median.data = median(temp,2); +spsdlq25 = quantile(temp, 0.25, 2); +spsdlq75 = quantile(temp, 0.75, 2); + +plot_spectra(spsdl_median,fmin,fmax) +patch([spsdl.freq' fliplr(spsdl.freq')], [spsdlq25' fliplr(spsdlq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3) +legend('Median','Quantiles') + +%% Window Aggregating +% +% If desired, one can group a series of measurements into blocks of time, though +% this isn't required by the IEC-40. In the following block, we'll take our 5 +% minutes of measurements, |time_aggregate| them into 30 second intervals, and +% find the median, 25% and 75% quantiles of each interval. We then plot the stats +% of the median parameters. + +window = 30; %seconds +method = "median"; +spsdl_time = time_aggregate(spsdl, window, method) + +%% +% We can then use the |plot_spectra| function in the graphics submodule to plot +% the median and quantiles of the first 30s window. + +spsdl_med = spsdl_time; +spsdl_med.data = median(spsdl_time.data, 2); +spsdlq25 = quantile(spsdl_time.data, 0.25, 2); +spsdlq75 = quantile(spsdl_time.data, 0.75, 2); + +plot_spectra(spsdl_med, fmin, fmax) +patch([spsdl.freq' fliplr(spsdl.freq')], [spsdlq25' fliplr(spsdlq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3) +legend(method,'Quantiles') + +%% Frequency Band Analysis +% +% Frequency band analysis can be completed by grouping the data into specific +% frequency bands, called "band aggregating" here. In other words, instead of +% aggregating by the time dimension, we aggregate by the frequency dimension. +% The |band_aggregate| function operates by taking the SPSDL and grouping it based +% on a specified octave and octave base. +% +% In the following section, we use the same plotting function as above, but +% do so by creating the decidecade frequency bands (10th octave, octave base 10 +% => $10^{1/10}$). We'll calculate the median and quantiles of each sample, then +% take the median of each of those parameters. For the IEC-40, one will calculate +% these parameters for each device operating state (current energy converters) +% or wave state matrix (wave energy converters). + +octave = [10, 10]; % [octave, octave base] +method = "median"; +spsdl10 = band_aggregate(spsdl,octave,fmin,fmax,method) +spsdl10m = spsdl10; +spsdl10m.data = median(spsdl10.data, 2); +spsdlq25 = quantile(spsdl10.data, 0.25, 2); +spsdlq75 = quantile(spsdl10.data, 0.75, 2); + +plot_spectra(spsdl10m, fmin, fmax) +patch([spsdl10.freq fliplr(spsdl10.freq)], [spsdlq25' fliplr(spsdlq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3) +legend(method,'Quantiles') + +%% +% The plot above shows significant spread in sound measurements at higher frequency +% due to the vessel noise, but less so at lower frequency. This mirrors what we +% can see in the sound pressure density level figure we first plotted. +% +% You can further build on the results with more complicated statistical methods. +% In the following block of code, we find the empirical quantile function (the +% empirical version of the cumulative distribution function, CDF) of each decidecade +% frequency band and plot the decidecade band centered nearest to 160 Hz. + +f = 160; % target frequency +[val, idxf] = min(abs(spsdl10.freq-f)); +n = numel(spsdl10.data(1,:)); +p = 1.0 - (0:(n-1))' / (n + 1); +x = sort(spsdl10.data(idxf,:)); +plot(x,p) +xlabel('SPSDL [dB re 1 uPa^2/Hz]'); ylabel("Exceedance Probability"); +title(strcat(num2str(spsdl10.freq(idxf))," Hz")) +grid on; + +%% Sound Pressure Level +% +% The IEC-40 has two requirements considering calculations of sound pressure +% level (SPL). We'll first calculate the SPL over the full frequency range of +% the turbine and/or hydrophone using the function |sound_pressure_level|. First, +% however, note that the IEC-40 asks that the range be set from 10 to 100,000 +% Hz. The lower limit can be increased due to flow noise or low frequency signal +% loss due to shallow water. +% Shallow water cutoff frequency +% Low frequency sound is absorbed into the seabed in shallow water depths. We +% can use the function |minimum_frequency| to get an approximation of what our +% minimum frequency should be. This approximation uses the water depth, estimates +% of the in-water sound speed, and sea/riverbed sound speed to determine what +% the cutoff frequency will be. The difficult part with this approximation is +% figuring out the speed of sound in the bed material, which generally ranges +% from 1450-1800 m/s. +% +% This function should only be used as a rough approximation and sanity check +% if significant attenuation is seen at various low frequencies and harmonics. + +depths = [1, 5, 10, 20, 40, 80]; +c = 1500; % speed of sound in water +c_seabed = 1700; % speed of sound in seabed +fmin = minimum_frequency(depths, c, c_seabed); +% print out min freq for corresponding depths +for i = 1:length(depths) + fprintf('%d m = %.1f Hz\n', depths(i), fmin(i)); +end + +%% +% Though the IEC-40 says we should default to a minimum frequency of 10 Hz, +% as you can see above, unless we're measuring from a depth of around 80 +/- 10 +% m, our minimum frequency should be higher. One can play around with the bed +% soundspeed to see how these change with varying bed densities/compositions. +% +%% Flow noise +% +% Flow noise, or psuedo-sound, is the other reason to increase the minimum frequency +% of our SPL measurements. Flow noise is caused by one of three things: turbulence +% advected past the hydrophone element, turbulence caused by the hydrophone element, +% and the sensitivity of the hydrophone element to temperature inhomogeneities +% in the advected flow. Flow noise is most noticeably apparent when flow speeds +% increase above 0.5 m/s, seen in spectrograms as a logarithmic increase in pressure +% with decreasing frequency. +% +% The particular data shown here was measured in around 8-10 m of water, and +% a mix of mild flow noise below 20 Hz and low frequency attenuation below ~50 +% Hz can be seen in the spectrogram. We'll again use the Nyquist frequency of +% 48,000 Hz. +% +%% Cumulative SPL +% +% Running the code block below, we can see our cumulative SPL start out at 86 +% dB and then peak at 125 dB as the boat drives by. If you haven't listened to +% the audio track, this peak SPL of 125 dB rel 1 uPa (underwater) is equivalent +% to 99 dB rel 20 uPa (air). For reference, the OSHA time limit for workers experiencing +% 100 dB rel 20 uPa of sound is 2 hours. Vessel traffic can be quite loud and +% is one of the largest contributors to noise in the marine environment. + +% sound pressure level +fmin = 50; +fmax = 48000; +spl = sound_pressure_level(spsd_cal, fmin, fmax); +% make figure +plot(spl.time,spl.data) +xlabel('Time'); ylabel({spl.name; spl.units}); grid on; + +%% Decidecade Sound Pressure Levels +% +% The last stat that IEC-40 requests are the decidecade SPLs, where the SPL +% is calculated for each decidecade frequency band. Note that the IEC-40 labels +% these as synonymous the third-octave SPLs; while for some applications this +% may be true, mathematically they have different definitions. +% +% To explain, a true octave is a frequency band where the upper frequency is +% double (i.e., base 2) that of the lower frequency. Third octaves are often measured +% because mammals to have evolved to interpret sound at this bandwidth. The true +% one-third octave is a frequency band where the upper frequency is 2^(1/3) = +% 1.25992 times the lower frequency. The decidecade band referenced by the IEC-40 +% refers to the one-tenth octave of base 10, where the upper frequency is 10 times +% that of the lower frequency. Mathematically this means the decidecade band has +% a bandwidth of 10^(1/10) = 1.25892. So, when reporting frequency analysis, it +% is important to note both the octave and its bandwidth. +% +% We can calculate the SPL in each decidecade band using the function |decidecade_sound_pressure_level|. +% This function uses the same calculation as |sound_pressure_level| above and +% runs it on each tenth octave band. It returns 1 SPL in each frequency band every +% timestamp. Boats are loud. + +spl10 = decidecade_sound_pressure_level(spsd_cal, fmin, fmax) +spl10m = spl10; +spl10m.data = median(spl10.data, 2); +splq25 = quantile(spl10.data, 0.25, 2); +splq75 = quantile(spl10.data, 0.75, 2); +plot_spectra(spl10m, fmin, fmax) +patch([spl10.freq fliplr(spl10.freq)], [splq25' fliplr(splq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3) +legend('Median','Quantiles') + +%% Third Octave Sound Pressure Level +% +% One can also calculate the true 1/3 octave SPLs using |third_octave_sound_pressure_level| +% if desired. Note the results are quite similar to |decidecade_sound_pressure_level|. + +spl3 = decidecade_sound_pressure_level(spsd_cal, fmin, fmax) +spl3m = spl3; +spl3m.data = median(spl3.data, 2); +splq25 = quantile(spl3.data, 0.25, 2); +splq75 = quantile(spl3.data, 0.75, 2); +plot_spectra(spl3m, fmin, fmax) +patch([spl3.freq fliplr(spl3.freq)], [splq25' fliplr(splq75')],'b','EdgeColor', 'none', 'FaceAlpha', 0.3) +legend('Median','Quantiles') + +%% Sound Exposure Level and Marine Mammal Weightings +% +% While the IEC-40 does not have recommendations for sound exposure level (SEL), +% it is a better parameter for quantifying auditory stress/injury to organisms. +% Sound exposure level is a metric that quantifies sound pressure level relative +% to exposure duration, based on the idea that that the effects on an organism's +% hearing are the same for equivalent totals sound energy received - even if that +% energy was received at different rates. Numerically, an SPL of 190 dB re 1 uPa +% over a 1 second interval and an SPL 160 dB re 1 uPa for 1000 seconds both have +% a SEL of 190 dB re 1 uPa^2 s. +% +% The National Marine Fisheries Service (NMFS) publishes auditory weighting +% functions for five groups of marine mammals. These weighting functions are designed +% to try to emulate the auditory sensitivity of each group to better predict auditory +% injury thresholds. They are mathematically equivalent to band-pass filters. +% A link to the latest recommendations from NMFS is +% +% The SEL function |sound_exposure_level| has a few inputs. One is the pressure +% spectral density (PSD) in Pa^2/Hz, found from |sound_pressure_spectral_density|. +% The 'n_bin' input to that function should be the same length of time as you +% want to calculate SEL for. If you would like to calculate the 24 hr SEL, you +% need 24 hours worth of data - the best way to do this and save computer RAM +% would be to calculate the PSD from each 1-5 minute individual file and then +% concatenate the PSD from each file together. +% +% Set "group=None" to calculate the standard SEL. Set "group" to "LF" for 'low +% frequency' cetaceans; "HF" for 'high frequency' cetaceans; "VHF" for 'very high +% frequency' cetaceans; "PW" for phocid pinnepeds; and "OW" for otariid pinnepeds. + +% five minute SEL +spsd_300s = sound_pressure_spectral_density(V, V.fs, 300); +% calibrate PSD +fill_Sf = sense_curve(end, 2); % use last value in curve for higher frequencies +spsd_300s = apply_calibration(spsd_300s, sense_curve, fill_Sf); + +group = ["LF", "HF", "VHF", "PW", "OW"]; +sel = sound_exposure_level(spsd_300s, [], fmin, fmax); +fprintf('Group-Standard SEL = %.1f\n', sel.data) +for g = 1:length(group) + sel = sound_exposure_level(spsd_300s, group(g), fmin, fmax); + fprintf('Group-%s SEL = %.1f\n', group(g), sel.data) +end + +%% Marine Mammal Auditory Weighting Functions +% +% We can look at the specific marine mammal auditory weighting and noise +% exposure functions as well using nmfs_auditory_weighting, given a frequency +% vector and one of the mammal groups. It outputs the weighting function and +% the exposure function (the inverse of the former) in units of dB. To +% convert back to a unitless magnitude, use 10.^( / 10). The +% exposure function shows the SEL in dB at and above which temporary or +% permanent hearing damage can occur to an individual in the specified group. +% The minimum value in the exposure function is the known or estimated +% injury level for a given group. + +% Calculate weighting and exposure functions +[weight_func, exp_func] = nmfs_auditory_weighting(spsd_300s.freq, 'LF'); + +% Visualize weighting and exposure +figure('Position', [100, 100, 800, 400]); + +% Weighting Function subplot +subplot(1, 2, 1); +semilogx(spsd_300s.freq, weight_func, 'LineWidth', 2, 'DisplayName', 'Weighting Function'); +hold on; +yline(0, 'k--', 'LineWidth', 1.5, 'DisplayName', 'Highest Sensitivity'); +xlabel('Frequency [Hz]'); +ylabel('Transmission [dB]'); +ylim([-50, 20]); +xlim([10, 48000]); +legend('Location', 'northeast'); +grid on; +title('Auditory Weighting Function'); + +% Exposure Function subplot +subplot(1, 2, 2); +semilogx(spsd_300s.freq, exp_func, 'LineWidth', 2, 'DisplayName', 'Exposure Function'); +hold on; +yline(min(exp_func), 'k--', 'LineWidth', 1.5, 'DisplayName', 'Highest Sensitivity'); + +% Fill injury region using fill() +% Note: For large datasets, this downsamples for visualization performance +if length(spsd_300s.freq) > 10000 + % Downsample for plotting efficiency while preserving shape + downsample_target_size = 5000; % target size for plotting + downsample_factor = ceil(length(spsd_300s.freq) / downsample_target_size); + freq_plot = spsd_300s.freq(1:downsample_factor:end); + exp_plot = exp_func(1:downsample_factor:end); +else + freq_plot = spsd_300s.freq; + exp_plot = exp_func; +end + +% Create filled region FROM exp_func UP TO level fill_max_y_level +fill_max_y_level = 300; +fill([freq_plot; flipud(freq_plot)], ... + [exp_plot; repmat(fill_max_y_level, size(exp_plot))], ... + 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none', 'DisplayName', 'Injury Region'); + +xlabel('Frequency [Hz]'); +ylabel('Exposure Level [dB]'); +ylim([min(exp_func)-20, min(exp_func)+50]); +xlim([10, 48000]); +legend('Location', 'southeast'); +grid on; +title('Auditory Exposure Function'); + +set(gcf, 'Color', 'white'); +sgtitle('NMFS Marine Mammal Auditory Functions (LF Group)'); diff --git a/examples/acoustics_example.mlx b/examples/acoustics_example.mlx new file mode 100644 index 000000000..3d55caee9 Binary files /dev/null and b/examples/acoustics_example.mlx differ diff --git a/examples/data/acoustics/6247.230204150508.wav b/examples/data/acoustics/6247.230204150508.wav new file mode 100644 index 000000000..5341f49d0 Binary files /dev/null and b/examples/data/acoustics/6247.230204150508.wav differ diff --git a/examples/data/acoustics/6247_calibration.csv b/examples/data/acoustics/6247_calibration.csv new file mode 100644 index 000000000..0bcda8f62 --- /dev/null +++ b/examples/data/acoustics/6247_calibration.csv @@ -0,0 +1,41 @@ +Frequency,Analog Sensitivity, +1,-223.49, +1.183,-220.8, +1.399,-218.13, +1.655,-215.41, +1.958,-212.68, +2.316,-209.91, +2.74,-207.12, +3.241,-204.29, +3.834,-201.45, +4.535,-198.58, +5.364,-195.69, +6.345,-192.79, +7.506,-189.85, +8.879,-186.9, +10.5,-183.93, +12.42,-180.93, +14.7,-177.92, +17.38,-174.97, +20.56,-172.16, +24.33,-169.61, +28.78,-167.69, +34.04,-166.52, +40.26,-165.96, +47.63,-165.81, +56.34,-165.85, +66.65,-165.95, +78.84,-166.05, +93.26,-166.13, +110.3,-166.2, +130.5,-166.25, +154.4,-166.28, +182.6,-166.29, +216,-166.29, +255.5,-166.27, +302.2,-166.24, +357.5,-166.17, +422.9,-166.03, +500.3,-165.79, +591.8,-165.47, +700,-164.87, diff --git a/examples/data/acoustics/RBW_6661_20240601_053114.wav b/examples/data/acoustics/RBW_6661_20240601_053114.wav new file mode 100644 index 000000000..08fdb27e6 Binary files /dev/null and b/examples/data/acoustics/RBW_6661_20240601_053114.wav differ diff --git a/mhkit/acoustics/analysis/apply_calibration.m b/mhkit/acoustics/analysis/apply_calibration.m new file mode 100644 index 000000000..18c8fe611 --- /dev/null +++ b/mhkit/acoustics/analysis/apply_calibration.m @@ -0,0 +1,63 @@ +function spsd_cal = apply_calibration(spsd,sensitivity_curve, fill_value) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Apply custom calibration to spectral density values +% +% Parameters +% ------------ +% spsd: struct +% Mean square sound pressure spectral density in V^2/Hz +% spsd.data : Spectral density data [V^2/Hz] +% spsd.freq : Frequency vector [Hz] +% spsd.time : Time vector (if time series) +% sensitivity_curve: matrix +% Calibrated sensitivity curve in units of dB rel 1 V^2/uPa^2 +% First column should be frequency, second column should be calibration values +% fill_value: float +% Value with which to fill missing values from the calibration curve [dB rel 1 V^2/uPa^2] +% +% Returns +% --------- +% spsd_cal: struct +% Calibrated spectral density in Pa^2/Hz, indexed by time and frequency +% spsd_cal.data : Calibrated spectral density data [Pa^2/Hz] +% spsd_cal.freq : Frequency vector [Hz] +% spsd_cal.time : Time vector (if time series) +% spsd_cal.name : Data name +% spsd_cal.units : Data units +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsd struct + sensitivity_curve {mustBeNumeric, mustBeFinite} + fill_value {mustBeNumeric, mustBeFinite} +end + +arguments (Output) + spsd_cal struct +end + +% Validate spsd structure +validate_spsd_struct(spsd, 'apply_calibration', ... + 'required_fields', {{'freq'}}); + +spsd_cal = spsd; + +% interpolate calibration +calibration = interp1(sensitivity_curve(:,1), ... + sensitivity_curve(:,2), spsd.freq, "linear"); +% use first cal value to fill NaNs in lower freq +idx = find(spsd.freq < sensitivity_curve(1,1)); +calibration(idx) = fillmissing(calibration(idx),'constant',sensitivity_curve(1,2)); +% fill NaNs at high freq using specified fill_value +calibration = fillmissing(calibration,'constant',fill_value); + +% subtract from sound pressure spectral density +sense_ratio = 10.^(calibration / 10); % V^2/uPa^2 +spsd_cal.data = (spsd_cal.data ./ sense_ratio) / 1e12; % Pa^2/Hz +spsd_cal.name = "Calibrated Sound Pressure Spectral Density"; +spsd_cal.units = "Pa^2/Hz"; + +end diff --git a/mhkit/acoustics/analysis/band_aggregate.m b/mhkit/acoustics/analysis/band_aggregate.m new file mode 100644 index 000000000..2f15e2b93 --- /dev/null +++ b/mhkit/acoustics/analysis/band_aggregate.m @@ -0,0 +1,98 @@ +function out = band_aggregate(spsdl, octave, fmin, fmax, method) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Reorganize spectral density level frequency tensor into fractional octave bands and applies a function to them +% +% Parameters +% ------------ +% spsdl: struct +% Mean square sound pressure spectral density level in dB rel 1 uPa^2/Hz +% spsdl.data : Spectral density level data [dB rel 1 uPa^2/Hz] +% spsdl.freq : Frequency vector [Hz] +% spsdl.time : Time vector +% spsdl.fs : Sampling frequency [Hz] +% octave: (1,2) vector +% Octave and octave base to subdivide spectral density level by +% Set octave base to 2 for true octave band; set to base 10 for decidecade octave band +% Default = [3, 2] (true third octave) +% fmin: float +% Lower frequency band limit (lower limit of the hydrophone) [Hz] +% fmax: float +% Upper frequency band limit (Nyquist frequency) [Hz] +% method: string or cell +% Method to run on the binned data. Can be a string (e.g., "median") or a cell +% where the first element is the method and the second is its argument +% Options: [median, mean, min, max, sum, quantile, std, var, count, map] +% +% Returns +% --------- +% out: struct +% Frequency band-averaged sound pressure spectral density level +% out.data : Band-averaged spectral density level [dB re 1 uPa^2/Hz] +% out.freq : Center frequencies of bands [Hz] +% out.time : Time vector +% out.name : Data name +% out.units : Data units +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +arguments (Input) + spsdl struct + octave {mustBeVector} = [3,2] + fmin {mustBeNumeric} = 10 + fmax {mustBeNumeric} = 100000 + method = "median" +end + +arguments (Output) + out struct +end + +% Validate spsdl structure +validate_spsdl_struct(spsdl, 'band_aggregate', ... + 'required_fields', {{'data', 'freq', 'time', 'name', 'units', 'fs'}}); + +if ~isnumeric(octave) || numel(octave) ~= 2 + error('MHKiT:acoustics:band_aggregate:InvalidOctave', 'octave must be a vector of two positive integers'); +end +if any(octave <= 0) + error('MHKiT:acoustics:band_aggregate:InvalidOctave', 'octave must contain positive integers'); +end +if ~isnumeric(fmin) || fmin <= 0 + error('MHKiT:acoustics:band_aggregate:InvalidFrequency', 'fmin must be a positive number'); +end +if ~isnumeric(fmax) || fmax <= fmin + error('MHKiT:acoustics:band_aggregate:InvalidFrequency', 'fmax must be greater than fmin'); +end + +fmax = fmax_warning(spsdl.fs/2, fmax); + +% validate method +[method_name, method_arg] = validate_method(method); +if isempty(method_arg) + mfunc = str2func(method_name); +elseif method_arg == "quantile" + mfunc = @(x)quantile(x,method_arg); +else + mfunc = method_arg; +end + +% derive octave bins +[octave_bins, band] = create_frequency_bands(octave(1),octave(2),fmin,fmax); + +% groupby and apply method +idx_bin = discretize(spsdl.freq, octave_bins); +temp = zeros(length(band.center_freq), length(spsdl.time)); +for x=1:length(spsdl.time) + temp(:,x) = splitapply(mfunc,spsdl.data(:,x),idx_bin); +end + +out.data = temp; +out.freq = band.center_freq(min(idx_bin):max(idx_bin)); +out.time = spsdl.time; +out.name = spsdl.name; +out.units = spsdl.units; + +end diff --git a/mhkit/acoustics/analysis/band_sound_pressure_level.m b/mhkit/acoustics/analysis/band_sound_pressure_level.m new file mode 100644 index 000000000..acf94d2f8 --- /dev/null +++ b/mhkit/acoustics/analysis/band_sound_pressure_level.m @@ -0,0 +1,81 @@ +function mspl = band_sound_pressure_level(spsd, octave, base, fmin, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate band-averaged sound pressure levels from the mean square sound pressure spectral density (SPSD) +% +% Parameters +% ------------ +% spsd: struct +% Mean square sound pressure spectral density in [Pa^2/Hz] +% spsd.data : Spectral density data [Pa^2/Hz] +% spsd.freq : Frequency vector [Hz] +% spsd.time : Time vector +% spsd.fs : Sampling frequency [Hz] +% octave: float +% Octave subdivision (1 = full octave, 3 = third-octave, etc.) +% base: float +% Octave base subdivision (2 = true octave, 10 = decade octave, etc.) +% fmin: float +% Lower frequency band limit (lower limit of the hydrophone) [Hz] +% fmax: float +% Upper frequency band limit (Nyquist frequency) [Hz] +% +% Returns +% --------- +% mspl: struct +% Sound pressure level [dB re 1 uPa] indexed by time and frequency of specified bandwidth +% mspl.data : Sound pressure level data [dB re 1 uPa] +% mspl.freq : Center frequencies of bands [Hz] +% mspl.time : Time vector +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsd struct + octave {mustBeInteger} + base {mustBeInteger} = 2 + fmin {mustBeNumeric} = 10 + fmax {mustBeNumeric} = 100000 +end + +arguments (Output) + mspl struct +end + +% Validate spsd structure +validate_spsd_struct(spsd, 'band_sound_pressure_level', ... + 'required_fields', {{'freq', 'time', 'fs'}}); + +reference = 1e-12; % Pa^2, = 1 uPa^2 + +fmax = fmax_warning(spsd.fs/2, fmax); + +% Create frequency bands +[obins, band] = create_frequency_bands(octave, base, fmin, fmax); + +nBands = numel(band.center_freq); +nTime = numel(spsd.time); +mspl.data = zeros(nBands, nTime); + +for i = 1:nBands + band_range = [band.lower_limit(i), band.upper_limit(i)]; + idx = find(spsd.freq >= band_range(1) & spsd.freq <= band_range(2)); + if numel(idx) < 2 + % Interpolate if only one frequency in band + spsd_slc = interp1(spsd.freq, spsd.data, band_range, 'linear', 'extrap'); + x = band_range; + else + spsd_slc = spsd.data(idx, :); + x = spsd.freq(idx); + end + % Integrate spectral density by frequency for each time + for t = 1:nTime + mspl.data(i, t) = 10 * log10(trapz(x, spsd_slc(:, t)) / reference); + end +end + +mspl.freq = band.center_freq; +mspl.time = spsd.time; + +end diff --git a/mhkit/acoustics/analysis/create_frequency_bands.m b/mhkit/acoustics/analysis/create_frequency_bands.m new file mode 100644 index 000000000..851568b36 --- /dev/null +++ b/mhkit/acoustics/analysis/create_frequency_bands.m @@ -0,0 +1,60 @@ +function [octave_bins, band] = create_frequency_bands(octave, base, fmin, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate frequency bands based on the specified octave, minimum and maximum frequency limits +% +% Parameters +% ------------ +% octave: double +% Octave to subdivide spectral density level by +% base: double, optional +% Octave base. Set to 2 for the true octave band; set to base 10 for +% the decidecade octave band. Default: 2 +% fmin: double, optional +% Lower frequency band limit (lower limit of the hydrophone). Default is 10 Hz +% fmax: double, optional +% Upper frequency band limit (Nyquist frequency). Default is 100,000 Hz +% +% Returns +% --------- +% octave_bins: array +% Array of octave bin edges +% band: structure +% band.center_freq : Center frequencies [Hz] +% band.lower_limit : Lower frequency limits [Hz] +% band.upper_limit : Upper frequency limits [Hz] +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + octave double {mustBeFinite, mustBePositive} + base double {mustBeFinite, mustBePositive} = 2 + fmin double {mustBeFinite, mustBePositive} = 10 + fmax double {mustBeFinite, mustBePositive} = 100000 +end + +arguments (Output) + octave_bins double + band struct +end + +bandwidth = base^(1 / octave); +half_bandwidth = base^(1 / (octave * 2)); + +% Calculate center frequencies +log_fmin = log10(fmin); +log_fmax = log10(fmax * bandwidth); +step = log10(bandwidth); +center_freq = 10 .^ (log_fmin : step : log_fmax); + +lower_limit = center_freq / half_bandwidth; +upper_limit = center_freq * half_bandwidth; +octave_bins = [lower_limit, upper_limit(end)]; + +band = struct(); +band.center_freq = center_freq; +band.lower_limit = lower_limit; +band.upper_limit = upper_limit; + +end diff --git a/mhkit/acoustics/analysis/decidecade_sound_pressure_level.m b/mhkit/acoustics/analysis/decidecade_sound_pressure_level.m new file mode 100644 index 000000000..5f16a453f --- /dev/null +++ b/mhkit/acoustics/analysis/decidecade_sound_pressure_level.m @@ -0,0 +1,61 @@ +function mspl = decidecade_sound_pressure_level(spsd, fmin, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate sound pressure level in decidecade bands directly +% from mean square sound pressure spectral density (SPSD). +% +% Parameters +% ------------ +% spsd: structure +% spsd.data : Mean square sound pressure spectral density [Pa^2/Hz] +% spsd.freq : Frequency vector [Hz] +% spsd.time : Time vector [s] +% spsd.fs : Sampling frequency [Hz] +% fmin: double +% Lower frequency band limit [Hz]. Default: 10 +% fmax: double +% Upper frequency band limit [Hz]. Default: 100000 +% +% Returns +% --------- +% mspl: structure +% mspl.data : Sound pressure level [dB re 1 uPa] +% mspl.freq : Decidecade band center frequencies [Hz] +% mspl.time : Time vector [s] +% mspl.units : Units string +% mspl.name : Descriptive name string +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +arguments (Input) + spsd struct + fmin {mustBeNumeric, mustBePositive} = 10 + fmax {mustBeNumeric, mustBePositive} = 100000 +end + +arguments (Output) + mspl struct +end + +% Validate spsd structure +validate_spsd_struct(spsd, 'decidecade_sound_pressure_level', ... + 'required_fields', {{'freq', 'time', 'fs'}}); + +% Validate frequency range +if fmin >= fmax + error('MHKiT:acoustics:decidecade_sound_pressure_level:InvalidInput', ... + 'fmin must be less than fmax'); +end + +fmax = fmax_warning(spsd.fs/2, fmax); + +octave = 10; +base = 10; + +mspl = band_sound_pressure_level(spsd, octave, base, fmin, fmax); +mspl.units = 'dB re 1 uPa'; +mspl.name = 'Decidecade Sound Pressure Level'; + +end diff --git a/mhkit/acoustics/analysis/fft_hann.m b/mhkit/acoustics/analysis/fft_hann.m new file mode 100644 index 000000000..b2fe9650e --- /dev/null +++ b/mhkit/acoustics/analysis/fft_hann.m @@ -0,0 +1,64 @@ +function [f, spec] = fft_hann(fs, x, nfft) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Apply Fast-Fourier-Transform to a time series using a Hanning window +% +% Parameters +% ------------ +% fs: double +% Sample rate of data [Hz] +% x: double array +% Array of time series data +% nfft: double +% Number of elements in the FFT +% +% Returns +% --------- +% f: double array +% Frequency array vector [Hz] +% spec: complex double array +% Frequency spectra resulting from the FFT +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + fs {mustBeNumeric, mustBePositive} + x {mustBeNumeric} + nfft {mustBeNumeric, mustBePositive, mustBeInteger} +end + +arguments (Output) + f {mustBeNumeric} + spec {mustBeNumeric} +end + +% Pre-process the input signal +x0 = x - mean(x); % Demean the time series +N = length(x); % Get length of time series + +% Generate and apply Hanning window +w = 0.5 * (1 - cos(2*pi*(0:1:N-1) / (N-1))); % Generate Hanning window coefficients +xx = x0 .* w'; % Apply Hanning window + +% Calculate FFT parameters +Fn = fs / 2; % Nyquist frequency +NumFFT = nfft; % Number of FFT points +NumUniquePts = ceil(NumFFT / 2); % Number of unique FFT points + +% Compute FFT +TempFFT = fft(xx, NumFFT); % Take FFT, padding with zeros if needed +spec = TempFFT(1:NumUniquePts); % Keep only positive frequencies + +% Adjust spectrum for single-sided representation +spec = spec * 2; % Account for negative frequencies +spec(1) = spec(1) / 2; % DC component should not be doubled +if mod(NumFFT, 2) == 0 % For even NumFFT + spec(end) = spec(end) / 2; % Nyquist frequency should not be doubled +end + +% Normalize spectrum and generate frequency vector +spec = spec / length(xx); % Normalize by signal length +f = (1:NumUniquePts) * 2 * Fn / NumFFT; % Generate frequency vector + +end diff --git a/mhkit/acoustics/analysis/fmax_warning.m b/mhkit/acoustics/analysis/fmax_warning.m new file mode 100644 index 000000000..79d646c61 --- /dev/null +++ b/mhkit/acoustics/analysis/fmax_warning.m @@ -0,0 +1,50 @@ +function fmax_out = fmax_warning(fn, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Check maximum frequency limit isn't greater than the Nyquist frequency +% +% Parameters +% ------------ +% fn: double (scalar) +% The Nyquist frequency [Hz] +% fmax: double (scalar) +% The maximum frequency limit [Hz] +% +% Returns +% --------- +% fmax_out: double (scalar) +% The adjusted maximum frequency limit, ensuring it does not exceed +% the Nyquist frequency [Hz] +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + fn {mustBeNumeric} + fmax {mustBeNumeric} +end + +arguments (Output) + fmax_out {mustBeNumeric} +end + +% Validate inputs are positive +if fn <= 0 + error('MHKiT:acoustics:fmax_warning:InvalidNyquist', ... + 'Nyquist frequency must be positive. Got fn = %.3f', fn); +end + +if fmax <= 0 + error('MHKiT:acoustics:fmax_warning:InvalidFmax', ... + 'Maximum frequency must be positive. Got fmax = %.3f', fmax); +end + +if fmax > fn + warning('MHKiT:acoustics:fmax_warning:FrequencyExceedsNyquist', ... + 'fmax = %.1f is greater than Nyquist frequency. Setting fmax = %.1f', fmax, fn); + fmax = fn; +end + +fmax_out = fmax; + +end diff --git a/mhkit/acoustics/analysis/minimum_frequency.m b/mhkit/acoustics/analysis/minimum_frequency.m new file mode 100644 index 000000000..c70752500 --- /dev/null +++ b/mhkit/acoustics/analysis/minimum_frequency.m @@ -0,0 +1,57 @@ +function fmin = minimum_frequency(water_depth, c, c_seabed) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Estimate the shallow water cutoff frequency based on the speed of +% sound in the water column and the speed of sound in the seabed +% material (generally ranges from 1450 - 1800 m/s) +% +% Parameters +% ------------ +% water_depth: double +% Depth of the water column [m] +% c: double, optional +% Speed of sound in the water column [m/s] (default: 1500) +% c_seabed: double, optional +% Speed of sound in the seabed material [m/s] (default: 1700) +% +% Returns +% --------- +% fmin: double +% The minimum cutoff frequency [Hz] +% +% Reference +% --------- +% Jennings 2011 - Computational Ocean Acoustics, 2nd ed. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + water_depth {mustBeNumeric} + c {mustBeNumeric} = 1500 + c_seabed {mustBeNumeric} = 1700 +end + +arguments (Output) + fmin {mustBeNumeric} +end + +water_depth = double(water_depth); + +if any(water_depth <= 0) + error('MHKiT:acoustics:minimum_frequency', 'All elements of water_depth must be positive numbers.'); +end +if c <= 0 + error('MHKiT:acoustics:minimum_frequency', 'c must be a positive number.'); +end +if c_seabed <= 0 + error('MHKiT:acoustics:minimum_frequency', 'c_seabed must be a positive number.'); +end +if c_seabed <= c + error('MHKiT:acoustics:minimum_frequency', 'c_seabed must be greater than c.'); +end + +fmin = c ./ (4 .* water_depth .* sqrt(1 - (c ./ c_seabed).^2)); + +end + diff --git a/mhkit/acoustics/analysis/nmfs_auditory_weighting.m b/mhkit/acoustics/analysis/nmfs_auditory_weighting.m new file mode 100644 index 000000000..c12e465de --- /dev/null +++ b/mhkit/acoustics/analysis/nmfs_auditory_weighting.m @@ -0,0 +1,66 @@ +function [weighting_func,exposure_func] = nmfs_auditory_weighting(frequency, group) +% +% Calculates the auditory weighting and exposure functions for marine mammals +% based on the National Marine Fisheries Service (NMFS) guidelines. +% +% The weighting function is applied to sound exposure level to determine the +% auditory impact on marine mammals. The exposure function is the inverse of the +% weighting function and illustrates how the weighting function relates to marine +% mammal hearing thresholds. +% Both function are returned in their log10-transform, in units of dB. To transform +% back to linear units, use 10**(weighting_func/10). +% +% https://www.fisheries.noaa.gov/national/marine-mammal-protection/marine-mammal-acoustic-technical-guidance-other-acoustic-tools +% +% Parameters +% ---------- +% frequency: array +% Frequency vector in [Hz]. +% group: str +% Marine mammal group for which the auditory weighting function is applied. +% Options: 'LF' (low frequency cetaceans), 'HF' (high frequency cetaceans), +% 'VHF' (very high frequency cetaceans), 'PW' (phocid pinnepeds), +% 'OW' (otariid pinnepeds) +% +% Returns +% ------- +% weighting_func: struct +% Auditory weighting function [unitless] indexed by frequency +% exposure_func: struct +% Log-transformed auditory exposure function [dB] indexed by frequency + +arguments (Input) + frequency {mustBeVector} + group {mustBeText} +end + +arguments (Output) + weighting_func {mustBeVector} + exposure_func {mustBeVector} +end + +group = lower(group); +valid_groups = {'lf', 'hf', 'vhf', 'pw', 'ow'}; +if ~ismember(group, valid_groups) + error('Group must be one of: LF, HF, VHF, PW, OW'); +end + +params = struct( ... + 'lf', struct('a', 0.99, 'b', 5, 'f1', 0.168, 'f2', 26.6, 'c', 0.12, 'k', 177), ... + 'hf', struct('a', 1.55, 'b', 5, 'f1', 1.73, 'f2', 129, 'c', 0.32, 'k', 181), ... + 'vhf', struct('a', 2.23, 'b', 5, 'f1', 5.93, 'f2', 186, 'c', 0.91, 'k', 160), ... + 'pw', struct('a', 1.63, 'b', 5, 'f1', 0.81, 'f2', 68.3, 'c', 0.29, 'k', 175), ... + 'ow', struct('a', 1.58, 'b', 5, 'f1', 2.53, 'f2', 43.8, 'c', 1.37, 'k', 178) ... +); + +p = params.(group); +frequency_kHz = frequency / 1000; % Convert to kHz +ratio_a = frequency_kHz / p.f1; +ratio_b = frequency_kHz / p.f2; +band_filter = ratio_a .^ (2 * p.a) ./ ( ((1 + ratio_a.^2) .^ p.a) .* ((1 + ratio_b.^2) .^ p.b) ); + +weighting_func = p.c + 10 * log10(band_filter); % dB +exposure_func = p.k - 10 * log10(band_filter); % dB + + +end \ No newline at end of file diff --git a/mhkit/acoustics/analysis/sound_exposure_level.m b/mhkit/acoustics/analysis/sound_exposure_level.m new file mode 100644 index 000000000..e7e1fc23d --- /dev/null +++ b/mhkit/acoustics/analysis/sound_exposure_level.m @@ -0,0 +1,110 @@ +function sel = sound_exposure_level(spsd, group, fmin, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate sound exposure level (SEL) across a specified frequency band +% from the sound pressure spectral density (SPSD). If a marine mammal group is +% provided, the resulting SEL is weighted according to the U.S. National Marine +% Fisheries Service (NMFS) guidelines. +% +% Parameters +% ------------ +% spsd: struct +% spsd.data : Sound pressure spectral density [Pa^2/Hz] +% spsd.freq : Frequency vector [Hz] +% spsd.time : Time vector [s] +% spsd.fs : Sampling frequency [Hz] +% spsd.nfft : Number of FFT points +% spsd.nbin : Bin length for exposure computation +% group: char or string +% Marine mammal group for auditory weighting function. Options: +% 'LF' (low frequency cetaceans), +% 'HF' (high frequency cetaceans), +% 'VHF' (very high frequency cetaceans), +% 'PW' (phocid pinnepeds), +% 'OW' (otariid pinnepeds). +% Default: [] +% fmin: double +% Lower frequency band limit [Hz]. Default: 10 +% fmax: double +% Upper frequency band limit [Hz]. Default: 100000 +% +% Returns +% --------- +% sel: struct +% sel.data : Sound exposure level [dB re 1 uPa^2 s] +% sel.time : Time vector [s] +% sel.units : Units string +% sel.name : Descriptive name +% sel.group : Marine mammal group used +% sel.integration_time : Integration time [s] +% sel.freq_band_min : Lower frequency limit [Hz] +% sel.freq_band_max : Upper frequency limit [Hz] +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsd struct + group = [] + fmin double {mustBeFinite, mustBePositive} = 10 + fmax double {mustBeFinite, mustBePositive} = 100000 +end + +arguments (Output) + sel struct +end + +% Validate spsd structure +validate_spsd_struct(spsd, 'sound_exposure_level', ... + 'required_fields', {{'freq', 'time', 'fs', 'nfft', 'nbin'}}); + +% Validate frequency bounds +if fmax <= fmin + error('MHKiT:acoustics:sound_exposure_level:InvalidInput', 'fmax must be greater than fmin'); +end + +% Validate group parameter if provided +if ~isempty(group) + valid_groups = {'LF', 'HF', 'VHF', 'PW', 'OW'}; + if ~ismember(group, valid_groups) + error('MHKiT:acoustics:sound_exposure_level:InvalidInput', ... + 'group must be one of: %s', strjoin(valid_groups, ', ')); + end +end +fmax = fmax_warning(spsd.fs/2, fmax); + +% get weighting factor +if ~isempty(group) + [w, e] = nmfs_auditory_weighting(spsd.freq, group); + w = 10 .^ (w/10); + name = "Weighted Sound Exposure Level"; +else + w = ones(length(spsd.freq), 1); + name = "Sound Exposure Level"; +end + +% Reference value of sound pressure +reference = 1e-12 * 1; % Pa^2 s, = 1 uPa^2 s + +% Select frequency band +idx = spsd.freq >= fmin & spsd.freq <= fmax; +band = spsd.data(idx, :); +freqs = spsd.freq(idx); +w = w(idx); + +exposure = trapz(freqs, band.*w); + +% Sound exposure level (L_{E,p}) = (L_{p,rms} + 10log10(t)) +s = 10 * log10(exposure / reference) + 10 * log10(spsd.nfft / spsd.fs); % n_points / (n_points/s) + +sel.data = s; +sel.time = spsd.time; +sel.units = "dB re 1 uPa^2 s"; +sel.name = name; +sel.group = group; +sel.integration_time = spsd.nbin; +sel.freq_band_min = fmin; +sel.freq_band_max = fmax; + + +end diff --git a/mhkit/acoustics/analysis/sound_pressure_level.m b/mhkit/acoustics/analysis/sound_pressure_level.m new file mode 100644 index 000000000..a6183341a --- /dev/null +++ b/mhkit/acoustics/analysis/sound_pressure_level.m @@ -0,0 +1,75 @@ +function spl = sound_pressure_level(spsd, fmin, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate sound pressure level (SPL) in a specified frequency band from +% the mean square sound pressure spectral density (SPSD) +% +% Parameters +% ------------ +% spsd: struct +% spsd.data : Mean square sound pressure spectral density [Pa^2/Hz] +% spsd.freq : Frequency vector [Hz] +% spsd.time : Time vector +% spsd.fs : Sampling frequency [Hz] +% fmin: integer +% Lower frequency band limit [Hz] (default: 10) +% fmax: integer +% Upper frequency band limit [Hz] (default: 100000) +% +% Returns +% --------- +% spl: struct +% spl.data : Sound pressure level [dB re 1 uPa] +% spl.time : Time vector +% spl.name : Data name identifier +% spl.units : Data units string +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsd struct + fmin {mustBeInteger} = 10 + fmax {mustBeInteger} = 100000 +end + +arguments (Output) + spl struct +end + +% Validate spsd structure +validate_spsd_struct(spsd, 'sound_pressure_level', ... + 'required_fields', {{'freq', 'time', 'fs'}}); + +% Additional parameter validation +if fmin <= 0 + error('MHKiT:acoustics:sound_pressure_level:InvalidInput', ... + 'fmin must be positive'); +end +if fmax <= fmin + error('MHKiT:acoustics:sound_pressure_level:InvalidInput', ... + 'fmax must be greater than fmin'); +end +fmax = fmax_warning(spsd.fs/2, fmax); + +% Reference value of sound pressure +reference = 1e-12; % Pa^2, = 1 uPa^2 + +% Select frequency band +idx = spsd.freq >= fmin & spsd.freq <= fmax; +band = spsd.data(idx, :); +freqs = spsd.freq(idx); + +% Integrate mean square sound pressure over frequency band +pressure_squared = trapz(freqs, band); + +% Mean square sound pressure level +mspl = 10 * log10(pressure_squared / reference); + +% Output +spl.data = mspl; +spl.time = spsd.time; +spl.name = 'Sound Pressure Level'; +spl.units = 'dB re 1 uPa'; + +end diff --git a/mhkit/acoustics/analysis/sound_pressure_spectral_density.m b/mhkit/acoustics/analysis/sound_pressure_spectral_density.m new file mode 100644 index 000000000..7fe8cf1a5 --- /dev/null +++ b/mhkit/acoustics/analysis/sound_pressure_spectral_density.m @@ -0,0 +1,107 @@ +function spsd = sound_pressure_spectral_density(data, fs, bin_length) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate sound pressure spectral density (SPSD) from audio samples split into FFTs +% with a specified bin length in seconds, using Hanning windowing with 50% overlap. +% +% Parameters +% ------------ +% data: structure +% data.data : Sound pressure [Pa] or Voltage [V] +% data.time : Time vector [s] or datetime +% data.units : Data units string +% fs: double +% Data collection sampling rate [Hz] +% bin_length: double +% Length of time in seconds to create FFTs. Default: 1. +% +% Returns +% --------- +% spsd: structure +% spsd.data : Mean square sound pressure spectral density [Pa^2/Hz or V^2/Hz] +% spsd.time : Time vector for spectral density bins [s] or datetime +% spsd.freq : Frequency vector [Hz] +% spsd.name : Description string +% spsd.units : Units string [Pa^2/Hz or V^2/Hz] +% spsd.fs : Sampling frequency [Hz] +% spsd.nbin : Bin length string [s] +% spsd.overlap : Overlap percentage string +% spsd.nfft : Number of FFT points +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + data struct + fs {mustBeNumeric} + bin_length {mustBeNumeric} = 1 +end + +arguments (Output) + spsd struct +end + +% Validate data structure has required fields +if ~isfield(data, 'data') + error('MHKiT:acoustics:sound_pressure_spectral_density:MissingField', ... + 'data structure must contain data field'); +end + +if ~isfield(data, 'time') + error('MHKiT:acoustics:sound_pressure_spectral_density:MissingField', ... + 'data structure must contain time field'); +end + +% Validate field types +if ~isnumeric(data.data) + error('MHKiT:acoustics:sound_pressure_spectral_density:InvalidInput', ... + 'data.data must be numeric'); +end + +% Accept both numeric and datetime for time fields +if ~(isnumeric(data.time) || isdatetime(data.time)) + error('MHKiT:acoustics:sound_pressure_spectral_density:InvalidInput', ... + 'data.time must be numeric or datetime'); +end + +% Check dimensional consistency: data length must match time length +if length(data.data) ~= length(data.time) + error('MHKiT:acoustics:sound_pressure_spectral_density:DimensionMismatch', ... + 'data.data length (%d) must match data.time length (%d)', ... + length(data.data), length(data.time)); +end + +pressure = data.data; + +% cut into multiple 1-second samples, calculate PSD of each sample +win = bin_length*fs; % window length of each time series +step = 0.5*win; % overlap betewen each window +ns = floor((length(pressure)-win)/step)+1; % number of time series samples +nfft = win; % number of fft points +df = fs/nfft; % frequency resolution +nfreq= ceil((nfft)/2); % Next highest power of 2 greater than length(x). +freq = (1:nfreq)*2*fs/2/nfft; % frequency of spectrum +Pf2 = zeros(nfreq,ns); % mean-squared sound pressure spectral density + +for i=1:ns + sample = pressure((i-1)*step+1:(i-1)*step+win); + sample = sample-mean(sample); + t_power = sum(sample.^2)/length(sample); % mean squred sound pressure; power in time domain + [f,spec] = fft_hann(fs,sample, nfft); % spectrum + psd = spec.*conj(spec)/df/2; % PSD + f_power = sum(psd)*df; % power in frequency domain + psd_adj = psd*t_power/f_power; % adjust the amplitude of PSD according to Parseval's theorem + Pf2(:,i) = psd_adj; %mean-squared sound pressure spectral density +end + +spsd.data = Pf2; +spsd.time = linspace(data.time(1), data.time(end), ns); +spsd.freq = freq'; +spsd.name = "Mean Square Sound Pressure Spectral Density"; +spsd.units = strcat(data.units, "^2/Hz"); +spsd.fs = fs; +spsd.nbin = strcat(string(bin_length)," s"); +spsd.overlap = "50%"; +spsd.nfft = nfft; + +end diff --git a/mhkit/acoustics/analysis/sound_pressure_spectral_density_level.m b/mhkit/acoustics/analysis/sound_pressure_spectral_density_level.m new file mode 100644 index 000000000..9ba341bff --- /dev/null +++ b/mhkit/acoustics/analysis/sound_pressure_spectral_density_level.m @@ -0,0 +1,49 @@ +function spsdl = sound_pressure_spectral_density_level(spsd) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate sound pressure spectral density level from mean square sound pressure spectral density +% +% Parameters +% ------------ +% spsd: struct +% spsd.data : Mean square sound pressure spectral density [Pa^2/Hz] +% spsd.freq : Frequency vector [Hz] (optional) +% spsd.time : Time vector (optional) +% +% Returns +% --------- +% spsdl: struct +% spsdl.data : Sound pressure spectral density level [dB re 1 uPa^2/Hz] +% spsdl.freq : Frequency vector [Hz] (if present in input) +% spsdl.time : Time vector (if present in input) +% spsdl.name : Data name identifier +% spsdl.units : Data units string +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsd struct +end + +arguments (Output) + spsdl struct +end + +% Validate spsd structure +validate_spsd_struct(spsd, 'sound_pressure_spectral_density_level', ... + 'require_positive', true, 'check_dimensions', false); + +% reference value of sound pressure +reference = 1e-12; % Pa^2 to 1uPa^2 + +% Sound pressure spectral density level from mean square values +lpf = 10 * log10(spsd.data / reference); + +% update struct +spsdl = spsd; +spsdl.data = lpf; +spsdl.name = "Sound Pressure Spectral Density Level"; +spsdl.units = "dB re 1 uPa^2/Hz"; + +end diff --git a/mhkit/acoustics/analysis/third_octave_sound_pressure_level.m b/mhkit/acoustics/analysis/third_octave_sound_pressure_level.m new file mode 100644 index 000000000..09b502654 --- /dev/null +++ b/mhkit/acoustics/analysis/third_octave_sound_pressure_level.m @@ -0,0 +1,53 @@ +function mspl = third_octave_sound_pressure_level(spsd, fmin, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculate sound pressure level in third octave bands directly +% from the mean square sound pressure spectral density (SPSD). +% +% Parameters +% ------------ +% spsd : structure +% Mean square sound pressure spectral density structure containing: +% spsd.S : Mean square sound pressure spectral density [Pa^2/Hz] +% spsd.f : Frequency vector [Hz] +% spsd.fs : Sampling frequency [Hz] +% fmin : double +% Lower frequency band limit (typically hydrophone lower limit) [Hz] +% fmax : double +% Upper frequency band limit (typically Nyquist frequency) [Hz] +% +% Returns +% --------- +% mspl : structure +% Third octave band sound pressure level structure containing: +% mspl.spl : Sound pressure level [dB re 1 μPa] +% mspl.f : Center frequencies of third octave bands [Hz] +% mspl.units : Units string 'dB re 1 uPa' +% mspl.name : Descriptive name 'Third Octave Sound Pressure Level' +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsd struct + fmin {mustBeNumeric, mustBePositive, mustBeFinite} = 10 + fmax {mustBeNumeric, mustBePositive, mustBeFinite} = 100000 +end + +arguments (Output) + mspl struct +end + +validate_spsd_struct(spsd, 'third_octave_sound_pressure_level', ... + 'required_fields', {{'fs'}}, 'check_dimensions', false); + +fmax = fmax_warning(spsd.fs/2, fmax); + +octave = 3; +base = 2; + +mspl = band_sound_pressure_level(spsd, octave, base, fmin, fmax); +mspl.units = 'dB re 1 uPa'; +mspl.name = 'Third Octave Sound Pressure Level'; + +end diff --git a/mhkit/acoustics/analysis/time_aggregate.m b/mhkit/acoustics/analysis/time_aggregate.m new file mode 100644 index 000000000..4a5c6dea4 --- /dev/null +++ b/mhkit/acoustics/analysis/time_aggregate.m @@ -0,0 +1,74 @@ +function out = time_aggregate(spsdl, window, method) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Reorganize spectral density level frequency tensor into time windows and apply a function to them. +% +% Parameters +% ------------ +% spsdl: structure +% spsdl.data : Sound pressure spectral density level data [dB re 1 uPa^2/Hz] +% spsdl.time : Time vector +% spsdl.freq : Frequency vector [Hz] +% spsdl.units : Data units string +% spsdl.name : Data name string +% window: integer +% Time in seconds to subdivide spectral density level into. Default: 60 s. +% method: string or structure +% Method to run on the binned data. Can be a string (e.g., "median") or a structure +% where the field is the method and the value is its argument (e.g., struct('quantile', 0.25)). +% Options: [median, mean, min, max, sum, quantile, std, var, count] +% +% Returns +% --------- +% out: structure +% out.data : Time-averaged sound pressure spectral density level [dB re 1 uPa^2/Hz] +% out.time : Time vector for binned data +% out.freq : Frequency vector [Hz] +% out.units : Data units string +% out.name : Data name string +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsdl struct + window {mustBeInteger} = 60 % seconds + method = "median" +end + +arguments (Output) + out struct +end + +% Validate spsdl structure +validate_spsdl_struct(spsdl, 'time_aggregate'); + +if window <= 0 + error('MHKiT:acoustics:time_aggregate:InvalidInput', ... + 'window must be a positive integer'); +end + +% validate method +[method_name, method_arg] = validate_method(method); +if isempty(method_arg) + mfunc = str2func(method_name); +elseif method_arg == "quantile" + mfunc = @(x)quantile(x,method_arg); +else + mfunc = method_arg; +end + +% groupby and apply method +tbins = spsdl.time(1):seconds(window):spsdl.time(end); +idx_bin = discretize(spsdl.time, tbins); +temp = zeros(length(spsdl.freq),length(tbins)-1); +for x=1:length(spsdl.freq) + temp(x,:) = splitapply(mfunc,spsdl.data(x,:),idx_bin); +end +out.data = temp; +out.time = tbins(1:end-1) + seconds(window)/2; +out.freq = spsdl.freq; +out.units = spsdl.units; +out.name = spsdl.name; + +end diff --git a/mhkit/acoustics/analysis/validate_method.m b/mhkit/acoustics/analysis/validate_method.m new file mode 100644 index 000000000..4e8d50378 --- /dev/null +++ b/mhkit/acoustics/analysis/validate_method.m @@ -0,0 +1,69 @@ +function [method_name,method_arg] = validate_method(method) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Helper function for validating the input of possible statistical methods +% +% Parameters +% ------------ +% method: string or cell array +% Name of statistical method or cell array with method and argument +% +% Returns +% --------- +% method_name: string +% Name of validated statistical method +% method_arg: string or double +% Corresponding method argument, if applicable +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + method +end + +arguments (Output) + method_name {mustBeText} + method_arg +end + +allowed_methods = [ + "all"; + "any"; + "numel"; + "cumsum"; + "map"; + "max"; + "mean"; + "min"; + "median"; + "prod"; + "quantile"; + "sum"; + "std"; + "var";]; + +if ~isstring(method) && ~iscell(method) + error('MHKiT:acoustics:validate_method', "'method' must be a string or a cell array") +end + +if isstring(method) + if ~ismember(lower(method), allowed_methods) + error('MHKiT:acoustics:validate_method', "'method' is not supported. See list of supported methods.") + end + method_name = lower(method); + method_arg = []; +end + +if iscell(method) + if ~isstring(method{1}) && ~ischar(method{1}) + error('MHKiT:acoustics:validate_method', "'method' must be a string or a cell array") + end + method_name = lower(method{1}); + if ~ismember(method_name, allowed_methods) + error('MHKiT:acoustics:validate_method', "'method' is not supported. See list of supported methods") + end + method_arg = method{2}; +end + +end diff --git a/mhkit/acoustics/graphics/plot_spectra.m b/mhkit/acoustics/graphics/plot_spectra.m new file mode 100644 index 000000000..d7946ed73 --- /dev/null +++ b/mhkit/acoustics/graphics/plot_spectra.m @@ -0,0 +1,53 @@ +function plot_spectra(spsdl, fmin, fmax) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Plot spectra of the sound pressure spectral density level +% +% Parameters +% ------------ +% spsdl: struct +% Mean square sound pressure spectral density level in dB rel 1 uPa^2/Hz +% spsdl.data : Spectral density level data [dB rel 1 uPa^2/Hz] +% spsdl.freq : Frequency vector [Hz] +% spsdl.time : Time vector +% spsdl.name : Data name +% spsdl.units : Data units +% fmin: float +% Lower frequency band limit (lower limit of the hydrophone) [Hz] +% fmax: float +% Upper frequency band limit (Nyquist frequency) [Hz] +% +% Returns +% --------- +% This function creates a figure but does not return a value +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsdl struct + fmin {mustBeNumeric} = 10 + fmax {mustBeNumeric} = 100000 +end + +% Validate spsdl structure +validate_spsdl_struct(spsdl, 'plot_spectra', ... + 'required_fields', {{'data', 'freq', 'name', 'units'}}, ... + 'check_dimensions', false); + +% check fmax +fn = max(spsdl.freq); +fmax = fmax_warning(fn, fmax); + +% select freq range +idxf = find(spsdl.freq>=fmin & spsdl.freq<=fmax); + +figure; +plot(spsdl.freq(idxf), spsdl.data(idxf,:)) +xscale log; +xlim([fmin fmax]) +xlabel('Frequency [Hz]'); +ylabel({spsdl.name;strcat('[',spsdl.units,']')}); +grid on; + +end diff --git a/mhkit/acoustics/graphics/plot_spectrogram.m b/mhkit/acoustics/graphics/plot_spectrogram.m new file mode 100644 index 000000000..c055d9441 --- /dev/null +++ b/mhkit/acoustics/graphics/plot_spectrogram.m @@ -0,0 +1,106 @@ +function plot_spectrogram(spsdl, options) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Plot spectrogram of the sound pressure spectral density level +% +% Parameters +% ------------ +% spsdl: struct +% spsdl.data : Spectral density level data [dB rel 1 uPa^2/Hz] +% spsdl.freq : Frequency vector [Hz] +% spsdl.time : Time vector +% spsdl.units : Data units +% options: name-value pairs (optional) +% fmin: float - Lower frequency limit [Hz] (default: 10) +% fmax: float - Upper frequency limit [Hz] (default: 100000) +% cm: string - Colormap name (default: 'hot') +% cmin: float - Minimum colorbar value (default: []) +% cmax: float - Maximum colorbar value (default: []) +% use_smoothing: logical - Apply median smoothing (default: true) +% smoothing_stride: integer - Smoothing kernel size (default: 3) +% plot_method: string - 'pcolor' (slowest), 'surf' (medium), or 'imagesc' (fastest) (default: 'imagesc') +% +% Returns +% --------- +% This function creates a figure but does not return a value +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments + spsdl struct + options.fmin {mustBeNumeric} = 10 + options.fmax {mustBeNumeric} = 100000 + options.cm {mustBeText} = 'hot' + options.cmin {mustBeNumeric} = [] + options.cmax {mustBeNumeric} = [] + options.use_smoothing {mustBeNumericOrLogical} = true + options.smoothing_stride {mustBeInteger, mustBePositive} = 3 + options.plot_method {mustBeText} = 'imagesc' +end + +% Validate spsdl structure +validate_spsdl_struct(spsdl, 'plot_spectrogram', ... + 'required_fields', {{'data', 'freq', 'time', 'units'}}, ... + 'check_dimensions', false); + +% Validate plot method +if ~ismember(options.plot_method, {'pcolor', 'surf', 'imagesc'}) + error('MHKiT:acoustics:plot_spectrogram:InvalidInput', ... + 'plot_method must be pcolor, surf, or imagesc'); +end + +% Check fmax +fn = max(spsdl.freq); +fmax = fmax_warning(fn, options.fmax); + +% Select frequency range +idxf = find(spsdl.freq>=options.fmin & spsdl.freq<=fmax); +data_subset = spsdl.data(idxf,:); + +% Apply smoothing if requested +if options.use_smoothing + data_subset = movmedian(movmedian(data_subset, options.smoothing_stride, 1), options.smoothing_stride, 2); +end + +figure; + +% Plot using specified method +switch options.plot_method + case 'imagesc' + imagesc(spsdl.time, spsdl.freq(idxf), data_subset); + set(gca, 'YDir', 'normal'); + case 'surf' + surf(spsdl.time, spsdl.freq(idxf), data_subset, 'EdgeColor', 'none'); + view(0, 90); + case 'pcolor' + pcolor(spsdl.time, spsdl.freq(idxf), data_subset); + shading interp; +end + +ylim([options.fmin fmax]); +set(gca, 'YScale', 'log'); +colormap(options.cm); +cb = colorbar(); + +% Set colorbar limits based on what's defined +if ~isempty(options.cmin) && ~isempty(options.cmax) + % Both limits defined + clim([options.cmin options.cmax]); +elseif ~isempty(options.cmin) + % Only minimum defined + current_limits = clim; + clim([options.cmin current_limits(2)]); +elseif ~isempty(options.cmax) + % Only maximum defined + current_limits = clim; + clim([current_limits(1) options.cmax]); +end +% If neither defined, use MATLAB default behavior + +xlabel('Time'); +ylabel('Frequency [Hz]'); +cb.Label.String = spsdl.units; + +end + diff --git a/mhkit/acoustics/input_validation/validate_spsd_struct.m b/mhkit/acoustics/input_validation/validate_spsd_struct.m new file mode 100644 index 000000000..1fcc8c955 --- /dev/null +++ b/mhkit/acoustics/input_validation/validate_spsd_struct.m @@ -0,0 +1,88 @@ +function validate_spsd_struct(spsd_struct, function_name, options) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Validate SPSD (Sound Pressure Spectral Density) structure format and content +% +% Parameters +% ------------ +% spsd_struct: struct +% SPSD structure to validate +% function_name: string +% Name of calling function for error messages +% options: struct (optional) +% require_positive : logical - Check data > 0 for log operations (default: false) +% check_dimensions : logical - Check dimensional consistency (default: true) +% required_fields : cell array - Additional required fields beyond 'data' (default: {}) +% +% Returns +% --------- +% No return value - throws error if validation fails +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsd_struct struct + function_name {mustBeText} + options.require_positive {mustBeNumericOrLogical} = false + options.check_dimensions {mustBeNumericOrLogical} = true + options.required_fields cell = {} +end + +% Always require 'data' field +base_required_fields = {'data'}; +all_required_fields = [base_required_fields, options.required_fields]; + +% Validate required fields exist +for i = 1:length(all_required_fields) + field_name = all_required_fields{i}; + if ~isfield(spsd_struct, field_name) + error(sprintf('MHKiT:acoustics:%s:MissingField', function_name), ... + 'spsd structure must contain %s field', field_name); + end +end + +% Validate data field +if ~isnumeric(spsd_struct.data) || ~all(isfinite(spsd_struct.data), 'all') + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsd.data must be numeric and finite'); +end + +% Check positive values if required (for log operations) +if options.require_positive && any(spsd_struct.data <= 0, 'all') + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsd.data must contain positive values for log calculation'); +end + +% Validate freq field if present +if isfield(spsd_struct, 'freq') + if ~isnumeric(spsd_struct.freq) || ~isvector(spsd_struct.freq) || any(spsd_struct.freq <= 0) + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsd.freq must be a positive numeric vector'); + end +end + +% Validate time field if present +if isfield(spsd_struct, 'time') + if ~(isnumeric(spsd_struct.time) || isdatetime(spsd_struct.time)) + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsd.time must be numeric or datetime'); + end +end + +% Check dimensional consistency if requested +if options.check_dimensions + if isfield(spsd_struct, 'freq') && size(spsd_struct.data, 1) ~= length(spsd_struct.freq) + error(sprintf('MHKiT:acoustics:%s:DimensionMismatch', function_name), ... + 'spsd.data rows (%d) must match spsd.freq length (%d)', ... + size(spsd_struct.data, 1), length(spsd_struct.freq)); + end + + if isfield(spsd_struct, 'time') && size(spsd_struct.data, 2) ~= length(spsd_struct.time) + error(sprintf('MHKiT:acoustics:%s:DimensionMismatch', function_name), ... + 'spsd.data columns (%d) must match spsd.time length (%d)', ... + size(spsd_struct.data, 2), length(spsd_struct.time)); + end +end + +end diff --git a/mhkit/acoustics/input_validation/validate_spsdl_struct.m b/mhkit/acoustics/input_validation/validate_spsdl_struct.m new file mode 100644 index 000000000..3c83b8d12 --- /dev/null +++ b/mhkit/acoustics/input_validation/validate_spsdl_struct.m @@ -0,0 +1,84 @@ +function validate_spsdl_struct(spsdl_struct, function_name, options) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Validate SPSDL (Sound Pressure Spectral Density Level) structure format and content +% +% Parameters +% ------------ +% spsdl_struct: struct +% SPSDL structure to validate +% function_name: string +% Name of calling function for error messages +% options: struct (optional) +% require_positive : logical - Check data > 0 for log operations (default: false) +% check_dimensions : logical - Check dimensional consistency (default: true) +% required_fields : cell array - Override default required fields (default: {'data', 'freq', 'time', 'name', 'units'}) +% +% Returns +% --------- +% No return value - throws error if validation fails +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + spsdl_struct struct + function_name {mustBeText} + options.require_positive {mustBeNumericOrLogical} = false + options.check_dimensions {mustBeNumericOrLogical} = true + options.required_fields cell = {'data', 'freq', 'time', 'name', 'units'} +end + +% Validate required fields exist +for i = 1:length(options.required_fields) + field_name = options.required_fields{i}; + if ~isfield(spsdl_struct, field_name) + error(sprintf('MHKiT:acoustics:%s:MissingField', function_name), ... + 'spsdl structure must contain %s field', field_name); + end +end + +% Validate data field +if ~isnumeric(spsdl_struct.data) || ~all(isfinite(spsdl_struct.data), 'all') + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsdl.data must be numeric and finite'); +end + +% Check positive values if required (for log operations) +if options.require_positive && any(spsdl_struct.data <= 0, 'all') + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsdl.data must contain positive values for log calculation'); +end + +% Validate freq field if present +if isfield(spsdl_struct, 'freq') + if ~isnumeric(spsdl_struct.freq) || ~isvector(spsdl_struct.freq) || any(spsdl_struct.freq <= 0) + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsdl.freq must be a positive numeric vector'); + end +end + +% Validate time field if present +if isfield(spsdl_struct, 'time') + if ~(isnumeric(spsdl_struct.time) || isdatetime(spsdl_struct.time)) + error(sprintf('MHKiT:acoustics:%s:InvalidInput', function_name), ... + 'spsdl.time must be numeric or datetime'); + end +end + +% Check dimensional consistency if requested +if options.check_dimensions + if isfield(spsdl_struct, 'freq') && size(spsdl_struct.data, 1) ~= length(spsdl_struct.freq) + error(sprintf('MHKiT:acoustics:%s:DimensionMismatch', function_name), ... + 'spsdl.data rows (%d) must match spsdl.freq length (%d)', ... + size(spsdl_struct.data, 1), length(spsdl_struct.freq)); + end + + if isfield(spsdl_struct, 'time') && size(spsdl_struct.data, 2) ~= length(spsdl_struct.time) + error(sprintf('MHKiT:acoustics:%s:DimensionMismatch', function_name), ... + 'spsdl.data columns (%d) must match spsdl.time length (%d)', ... + size(spsdl_struct.data, 2), length(spsdl_struct.time)); + end +end + +end diff --git a/mhkit/acoustics/io/export_audio.m b/mhkit/acoustics/io/export_audio.m new file mode 100644 index 000000000..42a3ff317 --- /dev/null +++ b/mhkit/acoustics/io/export_audio.m @@ -0,0 +1,101 @@ +function export_audio(filename, pressure_data, gain) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Export human-scaled audio data from hydrophone recording into WAV file +% +% Parameters +% ------------ +% filename: string +% Output filename for the WAV file (without extension) +% pressure_data: struct +% Sound pressure data structure containing: +% pressure_data.data : Pressure data array [Pa] +% pressure_data.sensitivity : Sensitivity of the hydrophone [dB re 1 V/uPa] +% pressure_data.fs : Sampling frequency [Hz] +% gain: float +% Gain to multiply the original time series by. Default: 1 +% +% Returns +% --------- +% None (writes WAV file to disk) +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + filename {mustBeText} + pressure_data struct + gain {mustBeNumeric} = 1 +end + +% Validate pressure_data structure has required fields +if ~isfield(pressure_data, 'data') + error('MHKiT:acoustics:export_audio:MissingField', ... + 'pressure_data structure must contain data field'); +end + +if ~isfield(pressure_data, 'sensitivity') + error('MHKiT:acoustics:export_audio:MissingField', ... + 'pressure_data structure must contain sensitivity field'); +end + +if ~isfield(pressure_data, 'fs') + error('MHKiT:acoustics:export_audio:MissingField', ... + 'pressure_data structure must contain fs field'); +end + +% Validate field types +if ~isnumeric(pressure_data.data) + error('MHKiT:acoustics:export_audio:InvalidInput', ... + 'pressure_data.data must be numeric'); +end + +if ~isnumeric(pressure_data.sensitivity) || ~isscalar(pressure_data.sensitivity) + error('MHKiT:acoustics:export_audio:InvalidInput', ... + 'pressure_data.sensitivity must be a numeric scalar'); +end + +if ~isnumeric(pressure_data.fs) || ~isscalar(pressure_data.fs) || pressure_data.fs <= 0 + error('MHKiT:acoustics:export_audio:InvalidInput', ... + 'pressure_data.fs must be a positive numeric scalar'); +end + +if ~isnumeric(gain) || ~isscalar(gain) + error('MHKiT:acoustics:export_audio:InvalidInput', ... + 'gain must be a numeric scalar'); +end + +% Process audio data following Python implementation +% Convert from Pascals to microPascals +upa = pressure_data.data * 1e6; + +% Convert to voltage waveform using sensitivity +v = upa * 10^(pressure_data.sensitivity / 20); % in V + +% Normalize to maximum absolute value and apply gain +max_abs_v = max(abs(v)); +if max_abs_v > 0 + v = v / max_abs_v * gain; +else + warning('MHKiT:acoustics:export_audio:ZeroSignal', ... + 'Signal has zero amplitude - no normalization applied'); +end + +% Ensure filename has .wav extension +[filepath, name, ext] = fileparts(filename); +if isempty(ext) + output_filename = fullfile(filepath, strcat(name, '.wav')); +else + output_filename = filename; +end + +% Write as 16-bit WAV file +try + audiowrite(output_filename, v, pressure_data.fs, 'BitsPerSample', 16); + fprintf('Audio exported successfully to: %s\n', output_filename); +catch ME + error('MHKiT:acoustics:export_audio:WriteError', ... + 'Failed to write audio file: %s', ME.message); +end + +end diff --git a/mhkit/acoustics/io/read_hydrophone.m b/mhkit/acoustics/io/read_hydrophone.m new file mode 100644 index 000000000..ddfd928dd --- /dev/null +++ b/mhkit/acoustics/io/read_hydrophone.m @@ -0,0 +1,84 @@ +function out = read_hydrophone(filename, peak_voltage, sensitivity, gain, start_time) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Read .wav file from a hydrophone. Returns voltage timeseries if +% sensitivity not provided, returns pressure timeseries if it is provided +% +% Parameters +% ------------ +% filename: string +% Input filename +% peak_voltage: float +% Peak voltage supplied to the analog to digital converter (ADC) [V] +% (Or 1/2 of the peak to peak voltage) +% sensitivity: float +% Hydrophone calibration sensitivity in dB re 1 V/uPa +% Should be negative. Default: None +% gain: float +% Amplifier gain in dB re 1 V/uPa +% start_time: string +% Start time in the format yyyy-MM-ddTHH:mm:ss +% +% Returns +% --------- +% out: struct +% Contains Sound pressure [Pa] or Voltage [V] along with metadata +% out.data : Time series data [Pa] or [V] +% out.time : Time vector +% out.units : Data units ('Pa' or 'V') +% out.fs : Sampling frequency [Hz] +% out.filename : Original filename +% out.valid_min : Minimum valid data value +% out.valid_max : Maximum valid data value +% out.sensitivity : Sensitivity value (if provided) +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + filename {mustBeText} + peak_voltage {mustBeNumeric} + sensitivity = [] + gain {mustBeNumeric} = 0 + start_time {mustBeText} = "2024-01-01T00:00:00" +end + +arguments (Output) + out struct +end + +[y, fs] = audioread(filename); % read in file +v = y * peak_voltage; % scale based on peak voltage + +out = struct(); +len = length(v) / fs; +start = datetime(start_time); +ending = start + seconds(len); +out.time = linspace(start, ending, length(v)); + +if ~isempty(sensitivity) + if sensitivity > 0 + error('MHKiT:acoustics:read_hydrophone:InvalidSensitivity', 'Sensitivity must be negative numeric input'); + end + sense = sensitivity - gain; + Sf = 10^(sense/20); + Pu = v / Sf; % sound pressure in uPa + P = Pu/10^6; % sound pressure in Pa + out.data = P; + out.units = 'Pa'; + out.sensitivity = Sf; + out.valid_min = -peak_voltage / Sf / 1e6; + out.valid_max = peak_voltage / Sf / 1e6; +else + out.data = v; % scaled voltage output + out.units = 'V'; + out.valid_min = -peak_voltage; + out.valid_max = peak_voltage; +end + +% add metadata +out.fs = fs; +[~,ff,ext] = fileparts(filename); +out.filename = strcat(ff,ext); + +end diff --git a/mhkit/acoustics/io/read_iclisten.m b/mhkit/acoustics/io/read_iclisten.m new file mode 100644 index 000000000..37be6c384 --- /dev/null +++ b/mhkit/acoustics/io/read_iclisten.m @@ -0,0 +1,82 @@ +function out = read_iclisten(filename, sensitivity, use_metadata) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Read .wav file from an icListen hydrophone +% +% Returns voltage timeseries if sensitivity not provided, returns pressure timeseries if it is provided +% +% Parameters +% ------------ +% filename: string +% Input filename +% sensitivity: float +% Hydrophone calibration sensitivity in dB re 1 V/uPa +% Should be negative. Default: None +% use_metadata: logical +% If true and 'sensitivity' is empty, applies sensitivity value +% stored in the .wav file. If false and 'sensitivity' is empty, a +% sensitivity value is not applied +% +% Returns +% --------- +% out: struct +% Contains Sound pressure [Pa] or Voltage [V] along with metadata +% out.data : Time series data [Pa] or [V] +% out.time : Time vector +% out.units : Data units ('Pa' or 'V') +% out.fs : Sampling frequency [Hz] +% out.filename : Original filename +% out.serial_number : Device serial number +% out.bits_per_sample : Bits per sample +% out.peak_voltage : Peak voltage [V] +% out.humidity : Humidity reading +% out.temperature : Temperature reading +% out.accelerometer : Accelerometer data +% out.magnetometer : Magnetometer data +% out.count_at_peak_voltage : Count at peak voltage +% out.sequence_num : Sequence number +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + filename {mustBeText} + sensitivity = [] + use_metadata logical = true +end + +arguments (Output) + out struct +end + +% decode metadata +metaraw = audioinfo(filename); +% get start_time +slen = strlength(metaraw.Title); +stime = datetime(metaraw.Title(slen-14:end),"InputFormat",'yyyyMMdd_HHmmss'); +start_time = string(stime,'yyyy-MM-dd HH:mm:ss'); + +meta = split(metaraw.Comment, ', '); +% get peak voltage +temp = split(meta{1},' '); +peak_voltage = str2double(temp{1}); +if isempty(sensitivity) & use_metadata + % get sensitivity + temp = split(meta{2},' '); + sensitivity = str2double(temp{1}); +end + +out = read_hydrophone(filename,peak_voltage,sensitivity, 0, start_time); + +% add extra metadata +out.serial_number = metaraw.Artist; +out.bits_per_sample = metaraw.BitsPerSample; +out.peak_voltage = peak_voltage; +out.humidity = meta{3}; +out.temperature = meta{4}; +out.accelerometer = meta{5}; +out.magnetometer = meta{6}; +out.count_at_peak_voltage = meta{7}; +out.sequence_num = meta{8}; + +end diff --git a/mhkit/acoustics/io/read_soundtrap.m b/mhkit/acoustics/io/read_soundtrap.m new file mode 100644 index 000000000..b8fec20ad --- /dev/null +++ b/mhkit/acoustics/io/read_soundtrap.m @@ -0,0 +1,54 @@ +function out = read_soundtrap(filename, sensitivity, gain) + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Read .wav file from an Ocean Instruments SoundTrap hydrophone +% +% Returns voltage timeseries if sensitivity not provided, returns pressure timeseries if it is provided +% +% Parameters +% ------------ +% filename: string +% Input filename +% sensitivity: float +% Hydrophone calibration sensitivity in dB re 1 V/uPa +% Should be negative. Default: None +% gain: float +% Amplifier gain in dB re 1 V/μPa +% +% Returns +% --------- +% out: struct +% Contains Sound pressure [Pa] or Voltage [V] along with metadata +% out.data : Time series data [Pa] or [V] +% out.time : Time vector +% out.units : Data units ('Pa' or 'V') +% out.fs : Sampling frequency [Hz] +% out.filename : Original filename +% out.valid_min : Minimum valid data value +% out.valid_max : Maximum valid data value +% out.sensitivity : Sensitivity value (if provided) +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + filename {mustBeText} + sensitivity = [] + gain {mustBeNumeric} = 0 +end + +arguments (Output) + out struct +end + +% parse start time from filename +s1 = strsplit(filename, '.'); +s1 = s1(end-1); +stime = datetime(s1, InputFormat="yyMMddHHmmss"); +start_time = string(stime, 'yyyy-MM-dd HH:mm:ss'); + +peak_voltage = 1; % soundtrap uses peak voltage of 1 V + +out = read_hydrophone(filename, peak_voltage, sensitivity, gain, start_time); + +end diff --git a/mhkit/tests/Acoustics_TestAnalysis.m b/mhkit/tests/Acoustics_TestAnalysis.m new file mode 100644 index 000000000..e77c0c3ce --- /dev/null +++ b/mhkit/tests/Acoustics_TestAnalysis.m @@ -0,0 +1,143 @@ +classdef Acoustics_TestAnalysis < matlab.unittest.TestCase + properties + spsd + end + + methods (TestClassSetup) + function setupOnce(testCase) + testdir = fileparts(mfilename('fullpath')); + datadir = fullfile(testdir, '..', '..', 'examples', 'data', 'acoustics'); + file_name = fullfile(datadir, '6247.230204150508.wav'); + P = read_soundtrap(file_name, -177); + testCase.spsd = sound_pressure_spectral_density(P, P.fs, 1); + end + end + + methods (Test) + function test_sound_pressure_spectral_density(testCase) + time = 0:0.1:9.9; + data = sin(time); + pressure.data = data; + pressure.time = time; + pressure.fs = 100; + pressure.units = 'Pa'; + + fs = 100; + bin_length = 0.1; + win_samples = bin_length * fs; + + spsd = sound_pressure_spectral_density(pressure, fs, bin_length); + + testCase.assertTrue(isstruct(spsd)); + testCase.assertTrue(isfield(spsd, 'freq')); + testCase.assertTrue(isfield(spsd, 'time')); + testCase.assertEqual(spsd.units, "Pa^2/Hz"); + testCase.assertEqual(spsd.nbin, sprintf("%g s", bin_length)); + + overlap = 0.5; % 50% overlap + expected_segments = floor(length(pressure.data) / (win_samples*(1-overlap))); + testCase.assertEqual(size(spsd.data, 2), expected_segments - 1); + end + + function test_apply_calibration(testCase) + time = 0:0.1:9.9; + freq = linspace(10, 1000, numel(time)); + spsd_data = rand(numel(time), numel(freq)); + spsd.data = spsd_data; + spsd.time = time; + spsd.freq = freq; + spsd.units = 'V^2/Hz'; + + f_curve = rand(1, numel(freq))'; + fill_value = 0.0; + + calibrated_spsd = apply_calibration(spsd, [freq',f_curve], fill_value); + + testCase.assertTrue(isstruct(calibrated_spsd)); + testCase.assertEqual(calibrated_spsd.units, "Pa^2/Hz"); + testCase.assertEqual(size(calibrated_spsd.data), size(spsd.data)); + testCase.verifyLessThan(max(calibrated_spsd.data), max(spsd.data)); + end + + function test_freq_loss(testCase) + fmin = minimum_frequency(20, 1500, 1700); + testCase.assertEqual(fmin, 39.84375); + end + + function test_spsdl(testCase) + td_spsdl = sound_pressure_spectral_density_level(testCase.spsd); + + cc = datetime([ ... + '2023-02-04 15:05:08.000' + '2023-02-04 15:05:08.502' + '2023-02-04 15:05:09.003' + '2023-02-04 15:05:09.505' + '2023-02-04 15:05:10.007'], 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSS'); + cd_spsdl = [ + 35.772061577247996 53.885951515794524 54.313180208582452 47.946115821970807 24.179883661640915; + 61.602936043528160 62.793717041169742 64.208043551995516 52.348715690865937 60.092868343630663; + 60.430110697601087 62.127736199381765 63.592162829344915 53.767504875422254 65.278154534206323; + 61.021732667297471 63.382134116643563 56.682598279488126 60.444088944449106 68.497868657610923; + 62.106024136858018 53.869424935555941 55.637077987595177 63.060370296217918 66.973991003254199; + ]; + + testCase.assertLessThan(abs(td_spsdl.data(1:5,1:5) - cd_spsdl), 1e-6); + testCase.assertLessThan(abs(posixtime(td_spsdl.time(1:5))' - posixtime(cc)), 1e-3); + end + + function test_averaging(testCase) + td_spsdl = sound_pressure_spectral_density_level(testCase.spsd); + + octave = [3, 2]; + td_spsdl_mean = band_aggregate(td_spsdl, octave, 10, 100000); + + lbin = 30; + td_spsdl_50 = time_aggregate(td_spsdl_mean, lbin, "median"); + + cc = datetime([ ... + '2023-02-04 15:05:23.000' + '2023-02-04 15:05:53.000' + '2023-02-04 15:06:23.000' + '2023-02-04 15:06:53.000' + '2023-02-04 15:07:23.000'], 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSS'); + cd_spsdl_50 = [ + 63.452642368099660 62.682385486447949 62.890011366673811 63.464187980216835 61.577799520029686; + 64.402822675991317 63.982645344031987 65.334381885187696 62.433832731393835 64.002193070627015; + 65.879232812149041 65.367523459295711 65.470768947676163 64.462081277290977 63.062659821297906; + 65.589916391755892 65.179564656250449 66.532335459481544 66.836901111597541 64.724432690177494; + 71.712416396688070 71.542148577084191 71.397701080122886 71.198748953358432 70.985008130759127 + ]; + + testCase.assertLessThan(abs(td_spsdl_50.data(1:5,1:5) - cd_spsdl_50), 1e-5); + testCase.assertLessThan(abs(posixtime(td_spsdl_50.time(1:5))' - posixtime(cc)), 1); + end + + function test_fmax_warning(testCase) + fn = 1000; + fmax = 1500; + adjusted_fmax = fmax_warning(fn, fmax); + testCase.assertEqual(adjusted_fmax, fn); + + fmax = 500; + adjusted_fmax = fmax_warning(fn, fmax); + testCase.assertEqual(adjusted_fmax, fmax); + + testCase.verifyError(@() fmax_warning('not a number', fmax), 'MATLAB:validators:mustBeNumeric'); + testCase.verifyError(@() fmax_warning(fn, 'not a number'), 'MATLAB:validators:mustBeNumeric'); + end + + function test_validate_method(testCase) + [method_name, method_arg] = validate_method("median"); + testCase.assertEqual(method_name, "median"); + testCase.assertEmpty(method_arg); + + [method_name, method_arg] = validate_method({"quantile", 0.25}); + testCase.assertEqual(method_name, "quantile"); + testCase.assertEqual(method_arg, 0.25); + + testCase.verifyError(@() validate_method('unsupported_method'), 'MHKiT:acoustics:validate_method'); + testCase.verifyError(@() validate_method(struct('unsupported_method', [])), 'MHKiT:acoustics:validate_method'); + testCase.verifyError(@() validate_method({'unsupported_method', 0.5}), 'MHKiT:acoustics:validate_method'); + end + end +end \ No newline at end of file diff --git a/mhkit/tests/Acoustics_TestIO.m b/mhkit/tests/Acoustics_TestIO.m new file mode 100644 index 000000000..bd04a2a60 --- /dev/null +++ b/mhkit/tests/Acoustics_TestIO.m @@ -0,0 +1,169 @@ +% filepath: c:\Github\MHKiT-Python\mhkit\tests\acoustics\test_io.m +classdef Acoustics_TestIO < matlab.unittest.TestCase + properties + datadir + plotdir + end + + methods (TestClassSetup) + function setupOnce(testCase) + testdir = fileparts(mfilename('fullpath')); + testCase.plotdir = fullfile(testdir, 'plots'); + if ~isfolder(testCase.plotdir) + mkdir(testCase.plotdir); + end + testCase.datadir = fullfile(testdir, '..', '..', 'examples', 'data', 'acoustics'); + end + end + + methods (Test) + + function test_read_iclisten_metadata(testCase) + file_name = fullfile(testCase.datadir, 'RBW_6661_20240601_053114.wav'); + metadata = read_iclisten(file_name); + + expected_metadata = struct(... + 'bits_per_sample', 24, ... + 'filename', 'RBW_6661_20240601_053114.wav', ... + 'peak_voltage', 3.0, ... + 'humidity', '24.0 % RH', ... + 'temperature', '8.6 deg C', ... + 'accelerometer', 'Acc(-980,-18,141)', ... + 'magnetometer', 'Mag(3603,3223,-598)', ... + 'count_at_peak_voltage', '8388608 = Max Count', ... + 'sequence_num', '2589798400000 = Seq #' ... + ); + + fields = fieldnames(expected_metadata); + for i = 1:numel(fields) + key = fields{i}; + expected_value = expected_metadata.(key); + testCase.verifyTrue(isfield(metadata, key)); + if isnumeric(expected_value) + testCase.verifyEqual(metadata.(key), expected_value, 'AbsTol', 1e-6); + else + testCase.verifyEqual(metadata.(key), expected_value); + end + end + end + + function test_read_iclisten(testCase) + file_name = fullfile(testCase.datadir, 'RBW_6661_20240601_053114.wav'); + td_orig = read_iclisten(file_name); + td_wrap = read_hydrophone(file_name, 3, -177, 0, "2024-06-01T05:31:14"); + td_volt = read_iclisten(file_name, [], false); + td_ovrrd = read_iclisten(file_name, -180, false); + td_ovrrd2 = read_iclisten(file_name, -180, true); + + cc = datetime([... + "2024-06-01T05:31:14.000000000" + "2024-06-01T05:31:14.000001953" + "2024-06-01T05:31:14.000003906" + "2024-06-01T05:31:14.000005859" + "2024-06-01T05:31:14.000007812"], 'InputFormat','yyyy-MM-dd''T''HH:mm:ss.SSSSSSSSS'); + + cd_orig = [0.31546374, 0.30229832, 0.32229963, 0.3159701, 0.30356423]; + cd_volt = [0.0004456, 0.00042701, 0.00045526, 0.00044632, 0.0004288]; + cd_ovrrd = [0.44560438, 0.42700773, 0.45526033, 0.44631963, 0.42879587]; + + testCase.verifyEqual(td_orig.data(1:5), cd_orig', 'AbsTol', 1e-6); + testCase.verifyLessThanOrEqual(td_orig.time(1:5)-cc', seconds(1e-9)); + + testCase.verifyEqual(td_wrap.data(1:5), cd_orig', 'AbsTol', 1e-6); + testCase.verifyLessThanOrEqual(td_wrap.time(1:5)-cc', seconds(1e-9)); + + testCase.verifyEqual(td_volt.data(1:5), cd_volt', 'AbsTol', 1e-6); + testCase.verifyLessThanOrEqual(td_volt.time(1:5)-cc', seconds(1e-9)); + + testCase.verifyEqual(td_ovrrd.data(1:5), cd_ovrrd', 'AbsTol', 1e-6); + testCase.verifyLessThanOrEqual(td_ovrrd.time(1:5)-cc', seconds(1e-9)); + + testCase.verifyEqual(td_ovrrd.data(1:5), td_ovrrd2.data(1:5), 'AbsTol', 1e-6); + end + + function test_read_soundtrap(testCase) + file_name = fullfile(testCase.datadir, '6247.230204150508.wav'); + td_orig = read_soundtrap(file_name, -177); + td_wrap = read_hydrophone(file_name, 1, -177, 0, '2023-02-04T15:05:08'); + td_volt = read_soundtrap(file_name, []); + + cc = datetime([... + "2023-02-04T15:05:08.000000000" + "2023-02-04T15:05:08.000010416" + "2023-02-04T15:05:08.000020832" + "2023-02-04T15:05:08.000031249" + "2023-02-04T15:05:08.000041665"], 'InputFormat','yyyy-MM-dd''T''HH:mm:ss.SSSSSSSSS'); + + cd_orig = [0.929006, 0.929006, 0.929006, 0.929006, 1.01542517]; + cd_volt = [0.00131226, 0.00131226, 0.00131226, 0.00131226, 0.00143433]; + + testCase.verifyEqual(td_orig.data(1:5), cd_orig', 'AbsTol', 1e-6); + testCase.verifyLessThanOrEqual(td_orig.time(1:5)-cc', seconds(1e-8)); + + testCase.verifyEqual(td_wrap.data(1:5), cd_orig', 'AbsTol', 1e-6); + testCase.verifyLessThanOrEqual(td_wrap.time(1:5)-cc', seconds(1e-8)); + + testCase.verifyEqual(td_volt.data(1:5), cd_volt', 'AbsTol', 1e-6); + testCase.verifyLessThanOrEqual(td_volt.time(1:5)-cc', seconds(1e-8)); + end + + function test_calibration(testCase) + file_name = fullfile(testCase.datadir, '6247.230204150508.wav'); + td_volt = read_soundtrap(file_name, []); + td_spsd = sound_pressure_spectral_density(td_volt, td_volt.fs, 1); + + cal_name = fullfile(testCase.datadir, '6247_calibration.csv'); + calibration = readtable(cal_name); + calibration = table2array(calibration); + fill_Sf = calibration(end, 2); % use last value in curve for higher frequencies + + td_spsd_cal = apply_calibration(td_spsd, calibration, fill_Sf); + + cc = datetime([ ... + '2023-02-04 15:05:08.000' + '2023-02-04 15:05:08.502' + '2023-02-04 15:05:09.003' + '2023-02-04 15:05:09.505' + '2023-02-04 15:05:10.007'], 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSS'); + + cd_spsd = [... + 0.000168347 0.010904234 0.012031451 0.002777228 1.16678E-05; + 0.004963596 0.006529418 0.009042888 0.000589351 0.00350582; + 0.000809366 0.001196485 0.001676293 0.000174536 0.002471435; + 0.000301447 0.000519099 0.000110993 0.000263904 0.00168587; + 0.000160878 2.41456E-05 3.62746E-05 0.000200416 0.000493508 + ]; + + testCase.verifyEqual(td_spsd_cal.data(1:5,1:5), cd_spsd, 'AbsTol', 1e-6); + testCase.verifyEqual(posixtime(td_spsd_cal.time(1:5)), posixtime(cc'), 'AbsTol', 1e-3); + end + + function test_export_audio(testCase) + file_name = fullfile(testCase.datadir, '6247.230204150508.wav'); + td_pressure = read_soundtrap(file_name, -177); + + % Create output filename in plots directory + output_filename = fullfile(testCase.plotdir, 'test_export_audio'); + + % Test export_audio with default gain + export_audio(output_filename, td_pressure); + + % Verify the output file was created + expected_wav_file = strcat(output_filename, '.wav'); + testCase.verifyTrue(exist(expected_wav_file, 'file') == 2, ... + 'Export audio should create a WAV file'); + + % Verify the audio file can be read back + [audio_data, fs_read] = audioread(expected_wav_file); + testCase.verifyEqual(fs_read, td_pressure.fs, ... + 'Exported audio sample rate should match input'); + testCase.verifyEqual(length(audio_data), length(td_pressure.data), ... + 'Exported audio length should match input'); + + % Clean up test file + if exist(expected_wav_file, 'file') + delete(expected_wav_file); + end + end + end +end diff --git a/mhkit/tests/Acoustics_TestMetrics.m b/mhkit/tests/Acoustics_TestMetrics.m new file mode 100644 index 000000000..76937bdac --- /dev/null +++ b/mhkit/tests/Acoustics_TestMetrics.m @@ -0,0 +1,149 @@ +classdef Acoustics_TestMetrics < matlab.unittest.TestCase + properties + spsd + spsd_60s + end + + methods (TestClassSetup) + function setupOnce(testCase) + testdir = fileparts(mfilename('fullpath')); + datadir = fullfile(testdir, '..', '..', 'examples', 'data', 'acoustics'); + file_name = fullfile(datadir, '6247.230204150508.wav'); + P = read_soundtrap(file_name, -177); + testCase.spsd = sound_pressure_spectral_density(P, P.fs, 1); + testCase.spsd_60s = sound_pressure_spectral_density(P, P.fs, 60); + end + end + + methods (Test) + function test_spl(testCase) + td_spl = sound_pressure_level(testCase.spsd, 10, 100000); + td_spl10 = decidecade_sound_pressure_level(testCase.spsd, 10, 100000); + td_spl3 = third_octave_sound_pressure_level(testCase.spsd, 10, 100000); + + cc = datetime([ ... + "2023-02-04 15:05:08.000000" + "2023-02-04 15:05:08.501683" + "2023-02-04 15:05:09.003367" + "2023-02-04 15:05:09.505051" + "2023-02-04 15:05:10.006735"], 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + cd_spl_head = [98.1252, 98.3950, 98.6400, 98.3964, 97.6300]; + cd_spl_tail = [97.43068, 97.67616, 97.99467, 98.13005, 97.96256]; + + cd_spl10_freq_head = [10.0, 12.589254, 15.848932, 19.952623, 25.118864]; + cd_spl10_head = [ + 67.698427107731632 68.718504549486369 55.150239136134928 63.092171166170736 68.285316609698185; + 70.638117562215356 68.352358818755633 70.679291272194590 66.826965002510178 72.307235015440625; + 75.061041066120794 72.000060887346237 61.867764868621322 69.029564442104871 65.566121037342981; + 73.571073481923875 77.646357542405568 71.315874138184853 68.348127773527210 74.641385774888377; + 83.318133484432281 82.382487625370942 82.839583740048681 80.066178195328149 82.922058814906620 + ]; + cd_spl10_freq_tail = [19952.62315, 25118.864315, 31622.776602, 39810.717055, 50118.723363]; + cd_spl10_tail = [ + 81.263857874463184 81.623440496905090 81.705362821962368 82.055563841890347 80.905483744451274; + 81.411568612914579 81.499648040350351 81.424401084480735 81.612600277281615 81.397403019626026; + 83.528128225083123 83.623441861246164 83.454856374486340 83.690034791494767 83.367940400050202; + 81.817594809639502 81.759206900301180 81.471336683096752 81.647659460075118 81.573892228707166; + 74.153354148096412 74.155760294865715 73.855762818092344 74.059185929338497 74.344841337560410 + ]; + cd_spl3_freq_head = [10.0, 12.59921, 15.874011, 20.0, 25.198421]; + cd_spl3_head = [ + 67.698427107731632 68.718504549486369 55.150239136134928 63.092171166170736 68.285316609698185; + 70.638117562215356 68.352358818755633 70.679291272194590 66.826965002510178 72.307235015440625; + 75.061041066120794 72.000060887346237 61.867764868621322 69.029564442104871 65.566121037342981; + 73.571073481923875 77.646357542405568 71.315874138184853 68.348127773527210 74.641385774888377; + 83.318133484432281 82.382487625370942 82.839583740048681 80.066178195328149 82.922058814906620 + ]; + cd_spl3_freq_tail = [20480.0, 25803.183102, 32509.973544, 40960.0, 51606.366204]; + cd_spl3_tail = [ + 81.134544889533359 81.479943775144307 81.531385268188657 81.881263944278700 80.703903519567945; + 81.677810946590085 81.775935532186026 81.694836721439614 81.865021419436431 81.690848660849610; + 83.968848694011569 83.970789087350326 83.796709880206919 84.065714823337430 83.760770230467969; + 80.663402401843385 80.686137257221461 80.384396424319945 80.522388118409964 80.532964378490249; + 72.079946592718045 72.278741310892528 71.948646690849614 71.850495761875621 72.240261721430514 + ]; + + testCase.verifyEqual(td_spl.data(1:5), cd_spl_head, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl.data(end-4:end), cd_spl_tail, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl10.freq(1:5), cd_spl10_freq_head, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl10.data(1:5,1:5), cd_spl10_head, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl10.freq(end-4:end), cd_spl10_freq_tail, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl10.data(end-4:end,end-4:end), cd_spl10_tail, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl3.freq(1:5), cd_spl3_freq_head, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl3.data(1:5,1:5), cd_spl3_head, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl3.freq(end-4:end), cd_spl3_freq_tail, 'AbsTol', 1e-3); + testCase.verifyEqual(td_spl3.data(end-4:end,end-4:end), cd_spl3_tail, 'AbsTol', 1e-3); + testCase.verifyEqual(posixtime(td_spl.time(1:5)), posixtime(cc'), 'AbsTol', 1e-5); + end + + function test_nmfs_weighting(testCase) + freq = testCase.spsd.freq; + slc = 20:25; % MATLAB is 1-based + + [W_LF, E_LF] = nmfs_auditory_weighting(freq, 'LF'); + [W_HF, E_HF] = nmfs_auditory_weighting(freq, 'HF'); + [W_VHF, E_VHF] = nmfs_auditory_weighting(freq, 'VHF'); + [W_PW, E_PW] = nmfs_auditory_weighting(freq, 'PW'); + [W_OW, E_OW] = nmfs_auditory_weighting(freq, 'OW'); + + cd_W_LF = [-18.241247, -17.827854, -17.434275, -17.058767, -16.699821, -16.3561]; + cd_E_LF = [195.36125, 194.94786, 194.55428, 194.17877, 193.81982, 193.4761]; + cd_W_HF = [-59.7284, -59.071625, -58.44541, -57.847057, -57.274178, -56.724693]; + cd_E_HF = [241.0484, 240.39163, 239.76541, 239.16705, 238.59418, 238.0447]; + cd_W_VHF = [-109.34241, -108.397385, -107.49632, -106.635315, -105.81097, -105.02029]; + cd_E_VHF = [270.2524, 269.30737, 268.4063, 267.54532, 266.72098, 265.9303]; + cd_W_PW = [-52.117348, -51.427025, -50.768852, -50.13999, -49.537937, -48.96051]; + cd_E_PW = [227.40735, 226.71703, 226.05885, 225.43, 224.82794, 224.25052]; + cd_W_OW = [-65.056496, -64.386955, -63.748577, -63.138584, -62.55456, -61.99438]; + cd_E_OW = [244.4265, 243.75696, 243.11858, 242.50858, 241.92456, 241.36438]; + + testCase.verifyEqual(W_LF(slc), cd_W_LF', 'AbsTol', 1e-4); + testCase.verifyEqual(W_HF(slc), cd_W_HF', 'AbsTol', 1e-4); + testCase.verifyEqual(W_VHF(slc), cd_W_VHF', 'AbsTol', 1e-4); + testCase.verifyEqual(W_PW(slc), cd_W_PW', 'AbsTol', 1e-4); + testCase.verifyEqual(W_OW(slc), cd_W_OW', 'AbsTol', 1e-4); + + testCase.verifyEqual(E_LF(slc), cd_E_LF', 'AbsTol', 1e-4); + testCase.verifyEqual(E_HF(slc), cd_E_HF', 'AbsTol', 1e-4); + testCase.verifyEqual(E_VHF(slc), cd_E_VHF', 'AbsTol', 1e-4); + testCase.verifyEqual(E_PW(slc), cd_E_PW', 'AbsTol', 1e-4); + testCase.verifyEqual(E_OW(slc), cd_E_OW', 'AbsTol', 1e-4); + end + + function test_sel(testCase) + td_sel = sound_exposure_level(testCase.spsd_60s,[], 10, 100000); + td_sel_lf = sound_exposure_level(testCase.spsd_60s,'LF', 10, 100000); + td_sel_hf = sound_exposure_level(testCase.spsd_60s,'HF', 10, 100000); + td_sel_vhf = sound_exposure_level(testCase.spsd_60s,'VHF', 10, 100000 ); + td_sel_pw = sound_exposure_level(testCase.spsd_60s,'PW', 10, 100000); + td_sel_ow = sound_exposure_level(testCase.spsd_60s,'OW', 10, 100000); + + cc = datetime([ ... + "2023-02-04 15:05:08.000000" + "2023-02-04 15:05:45.500874" + "2023-02-04 15:06:23.001750" + "2023-02-04 15:07:00.502625" + "2023-02-04 15:07:38.003500"], 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS'); + cd_sel = [1.161832567214382e+02 1.177608180325696e+02 1.216988394734190e+02 1.322482028621648e+02 1.432811329170941e+02]; + cd_sel_lf = [1.123633263210942e+02 1.150560720648516e+02 1.201771340004273e+02 1.315286929637698e+02 1.427492968479198e+02]; + cd_sel_hf = [1.122226246068385e+02 1.143795287244090e+02 1.188811398942907e+02 1.293576472503867e+02 1.399411804663656e+02]; + cd_sel_vhf = [1.102333042133257e+02 1.111969207330859e+02 1.140075568256765e+02 1.228731328944110e+02 1.332000231806384e+02]; + cd_sel_pw = [1.122223727653538e+02 1.148390026535553e+02 1.198729727976822e+02 1.308165812210123e+02 1.416746474465080e+02]; + cd_sel_ow = [1.109456496065601e+02 1.133797705458589e+02 1.180640501702054e+02 1.285908810877421e+02 1.390843177149019e+02]; + + testCase.verifyEqual(td_sel.data(1:5), cd_sel, 'AbsTol', 1e-3); + testCase.verifyEqual(td_sel_lf.data(1:5), cd_sel_lf, 'AbsTol', 1e-3); + testCase.verifyEqual(td_sel_hf.data(1:5), cd_sel_hf, 'AbsTol', 1e-3); + testCase.verifyEqual(td_sel_vhf.data(1:5), cd_sel_vhf, 'AbsTol', 1e-3); + testCase.verifyEqual(td_sel_pw.data(1:5), cd_sel_pw, 'AbsTol', 1e-3); + testCase.verifyEqual(td_sel_ow.data(1:5), cd_sel_ow, 'AbsTol', 1e-3); + testCase.verifyEqual(posixtime(td_sel.time(1:5)), posixtime(cc'), 'AbsTol', 1e-6); + end + + function test_spl_vs_sel(testCase) + td_spl = sound_pressure_level(testCase.spsd, 10, 100000); + td_sel = sound_exposure_level(testCase.spsd, [], 10, 100000); + testCase.verifyEqual(td_spl.data, td_sel.data, 'AbsTol', 1e-6); + end + end +end \ No newline at end of file