Acoustics module

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,

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

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

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

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