Skip to content

Commit 7e5ab9d

Browse files
simmsasimmsa
authored
Merge pull request #174 from hivanov-nrel/acoustics
Acoustics module
2 parents 2c76f49 + ac3132a commit 7e5ab9d

33 files changed

Lines changed: 3408 additions & 0 deletions

‎examples/acoustics_example.html‎

Lines changed: 622 additions & 0 deletions
Large diffs are not rendered by default.

‎examples/acoustics_example.m‎

Lines changed: 495 additions & 0 deletions
Large diffs are not rendered by default.

‎examples/acoustics_example.mlx‎

486 KB
Binary file not shown.

‎examples/data/acoustics/6247.230204150508.wav‎

54.9 MB
Binary file not shown.

‎examples/data/acoustics/6247_calibration.csv‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
Frequency,Analog Sensitivity,
2+
1,-223.49,
3+
1.183,-220.8,
4+
1.399,-218.13,
5+
1.655,-215.41,
6+
1.958,-212.68,
7+
2.316,-209.91,
8+
2.74,-207.12,
9+
3.241,-204.29,
10+
3.834,-201.45,
11+
4.535,-198.58,
12+
5.364,-195.69,
13+
6.345,-192.79,
14+
7.506,-189.85,
15+
8.879,-186.9,
16+
10.5,-183.93,
17+
12.42,-180.93,
18+
14.7,-177.92,
19+
17.38,-174.97,
20+
20.56,-172.16,
21+
24.33,-169.61,
22+
28.78,-167.69,
23+
34.04,-166.52,
24+
40.26,-165.96,
25+
47.63,-165.81,
26+
56.34,-165.85,
27+
66.65,-165.95,
28+
78.84,-166.05,
29+
93.26,-166.13,
30+
110.3,-166.2,
31+
130.5,-166.25,
32+
154.4,-166.28,
33+
182.6,-166.29,
34+
216,-166.29,
35+
255.5,-166.27,
36+
302.2,-166.24,
37+
357.5,-166.17,
38+
422.9,-166.03,
39+
500.3,-165.79,
40+
591.8,-165.47,
41+
700,-164.87,

‎examples/data/acoustics/RBW_6661_20240601_053114.wav‎

87.9 MB
Binary file not shown.

‎mhkit/acoustics/analysis/apply_calibration.m‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
function spsd_cal = apply_calibration(spsd,sensitivity_curve, fill_value)
2+
3+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4+
%
5+
% Apply custom calibration to spectral density values
6+
%
7+
% Parameters
8+
% ------------
9+
% spsd: struct
10+
% Mean square sound pressure spectral density in V^2/Hz
11+
% spsd.data : Spectral density data [V^2/Hz]
12+
% spsd.freq : Frequency vector [Hz]
13+
% spsd.time : Time vector (if time series)
14+
% sensitivity_curve: matrix
15+
% Calibrated sensitivity curve in units of dB rel 1 V^2/uPa^2
16+
% First column should be frequency, second column should be calibration values
17+
% fill_value: float
18+
% Value with which to fill missing values from the calibration curve [dB rel 1 V^2/uPa^2]
19+
%
20+
% Returns
21+
% ---------
22+
% spsd_cal: struct
23+
% Calibrated spectral density in Pa^2/Hz, indexed by time and frequency
24+
% spsd_cal.data : Calibrated spectral density data [Pa^2/Hz]
25+
% spsd_cal.freq : Frequency vector [Hz]
26+
% spsd_cal.time : Time vector (if time series)
27+
% spsd_cal.name : Data name
28+
% spsd_cal.units : Data units
29+
%
30+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31+
32+
arguments (Input)
33+
spsd struct
34+
sensitivity_curve {mustBeNumeric, mustBeFinite}
35+
fill_value {mustBeNumeric, mustBeFinite}
36+
end
37+
38+
arguments (Output)
39+
spsd_cal struct
40+
end
41+
42+
% Validate spsd structure
43+
validate_spsd_struct(spsd, 'apply_calibration', ...
44+
'required_fields', {{'freq'}});
45+
46+
spsd_cal = spsd;
47+
48+
% interpolate calibration
49+
calibration = interp1(sensitivity_curve(:,1), ...
50+
sensitivity_curve(:,2), spsd.freq, "linear");
51+
% use first cal value to fill NaNs in lower freq
52+
idx = find(spsd.freq < sensitivity_curve(1,1));
53+
calibration(idx) = fillmissing(calibration(idx),'constant',sensitivity_curve(1,2));
54+
% fill NaNs at high freq using specified fill_value
55+
calibration = fillmissing(calibration,'constant',fill_value);
56+
57+
% subtract from sound pressure spectral density
58+
sense_ratio = 10.^(calibration / 10); % V^2/uPa^2
59+
spsd_cal.data = (spsd_cal.data ./ sense_ratio) / 1e12; % Pa^2/Hz
60+
spsd_cal.name = "Calibrated Sound Pressure Spectral Density";
61+
spsd_cal.units = "Pa^2/Hz";
62+
63+
end

‎mhkit/acoustics/analysis/band_aggregate.m‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
function out = band_aggregate(spsdl, octave, fmin, fmax, method)
2+
3+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4+
%
5+
% Reorganize spectral density level frequency tensor into fractional octave bands and applies a function to them
6+
%
7+
% Parameters
8+
% ------------
9+
% spsdl: struct
10+
% Mean square sound pressure spectral density level in dB rel 1 uPa^2/Hz
11+
% spsdl.data : Spectral density level data [dB rel 1 uPa^2/Hz]
12+
% spsdl.freq : Frequency vector [Hz]
13+
% spsdl.time : Time vector
14+
% spsdl.fs : Sampling frequency [Hz]
15+
% octave: (1,2) vector
16+
% Octave and octave base to subdivide spectral density level by
17+
% Set octave base to 2 for true octave band; set to base 10 for decidecade octave band
18+
% Default = [3, 2] (true third octave)
19+
% fmin: float
20+
% Lower frequency band limit (lower limit of the hydrophone) [Hz]
21+
% fmax: float
22+
% Upper frequency band limit (Nyquist frequency) [Hz]
23+
% method: string or cell
24+
% Method to run on the binned data. Can be a string (e.g., "median") or a cell
25+
% where the first element is the method and the second is its argument
26+
% Options: [median, mean, min, max, sum, quantile, std, var, count, map]
27+
%
28+
% Returns
29+
% ---------
30+
% out: struct
31+
% Frequency band-averaged sound pressure spectral density level
32+
% out.data : Band-averaged spectral density level [dB re 1 uPa^2/Hz]
33+
% out.freq : Center frequencies of bands [Hz]
34+
% out.time : Time vector
35+
% out.name : Data name
36+
% out.units : Data units
37+
%
38+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39+
40+
41+
arguments (Input)
42+
spsdl struct
43+
octave {mustBeVector} = [3,2]
44+
fmin {mustBeNumeric} = 10
45+
fmax {mustBeNumeric} = 100000
46+
method = "median"
47+
end
48+
49+
arguments (Output)
50+
out struct
51+
end
52+
53+
% Validate spsdl structure
54+
validate_spsdl_struct(spsdl, 'band_aggregate', ...
55+
'required_fields', {{'data', 'freq', 'time', 'name', 'units', 'fs'}});
56+
57+
if ~isnumeric(octave) || numel(octave) ~= 2
58+
error('MHKiT:acoustics:band_aggregate:InvalidOctave', 'octave must be a vector of two positive integers');
59+
end
60+
if any(octave <= 0)
61+
error('MHKiT:acoustics:band_aggregate:InvalidOctave', 'octave must contain positive integers');
62+
end
63+
if ~isnumeric(fmin) || fmin <= 0
64+
error('MHKiT:acoustics:band_aggregate:InvalidFrequency', 'fmin must be a positive number');
65+
end
66+
if ~isnumeric(fmax) || fmax <= fmin
67+
error('MHKiT:acoustics:band_aggregate:InvalidFrequency', 'fmax must be greater than fmin');
68+
end
69+
70+
fmax = fmax_warning(spsdl.fs/2, fmax);
71+
72+
% validate method
73+
[method_name, method_arg] = validate_method(method);
74+
if isempty(method_arg)
75+
mfunc = str2func(method_name);
76+
elseif method_arg == "quantile"
77+
mfunc = @(x)quantile(x,method_arg);
78+
else
79+
mfunc = method_arg;
80+
end
81+
82+
% derive octave bins
83+
[octave_bins, band] = create_frequency_bands(octave(1),octave(2),fmin,fmax);
84+
85+
% groupby and apply method
86+
idx_bin = discretize(spsdl.freq, octave_bins);
87+
temp = zeros(length(band.center_freq), length(spsdl.time));
88+
for x=1:length(spsdl.time)
89+
temp(:,x) = splitapply(mfunc,spsdl.data(:,x),idx_bin);
90+
end
91+
92+
out.data = temp;
93+
out.freq = band.center_freq(min(idx_bin):max(idx_bin));
94+
out.time = spsdl.time;
95+
out.name = spsdl.name;
96+
out.units = spsdl.units;
97+
98+
end

‎mhkit/acoustics/analysis/band_sound_pressure_level.m‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
function mspl = band_sound_pressure_level(spsd, octave, base, fmin, fmax)
2+
3+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4+
%
5+
% Calculate band-averaged sound pressure levels from the mean square sound pressure spectral density (SPSD)
6+
%
7+
% Parameters
8+
% ------------
9+
% spsd: struct
10+
% Mean square sound pressure spectral density in [Pa^2/Hz]
11+
% spsd.data : Spectral density data [Pa^2/Hz]
12+
% spsd.freq : Frequency vector [Hz]
13+
% spsd.time : Time vector
14+
% spsd.fs : Sampling frequency [Hz]
15+
% octave: float
16+
% Octave subdivision (1 = full octave, 3 = third-octave, etc.)
17+
% base: float
18+
% Octave base subdivision (2 = true octave, 10 = decade octave, etc.)
19+
% fmin: float
20+
% Lower frequency band limit (lower limit of the hydrophone) [Hz]
21+
% fmax: float
22+
% Upper frequency band limit (Nyquist frequency) [Hz]
23+
%
24+
% Returns
25+
% ---------
26+
% mspl: struct
27+
% Sound pressure level [dB re 1 uPa] indexed by time and frequency of specified bandwidth
28+
% mspl.data : Sound pressure level data [dB re 1 uPa]
29+
% mspl.freq : Center frequencies of bands [Hz]
30+
% mspl.time : Time vector
31+
%
32+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33+
34+
arguments (Input)
35+
spsd struct
36+
octave {mustBeInteger}
37+
base {mustBeInteger} = 2
38+
fmin {mustBeNumeric} = 10
39+
fmax {mustBeNumeric} = 100000
40+
end
41+
42+
arguments (Output)
43+
mspl struct
44+
end
45+
46+
% Validate spsd structure
47+
validate_spsd_struct(spsd, 'band_sound_pressure_level', ...
48+
'required_fields', {{'freq', 'time', 'fs'}});
49+
50+
reference = 1e-12; % Pa^2, = 1 uPa^2
51+
52+
fmax = fmax_warning(spsd.fs/2, fmax);
53+
54+
% Create frequency bands
55+
[obins, band] = create_frequency_bands(octave, base, fmin, fmax);
56+
57+
nBands = numel(band.center_freq);
58+
nTime = numel(spsd.time);
59+
mspl.data = zeros(nBands, nTime);
60+
61+
for i = 1:nBands
62+
band_range = [band.lower_limit(i), band.upper_limit(i)];
63+
idx = find(spsd.freq >= band_range(1) & spsd.freq <= band_range(2));
64+
if numel(idx) < 2
65+
% Interpolate if only one frequency in band
66+
spsd_slc = interp1(spsd.freq, spsd.data, band_range, 'linear', 'extrap');
67+
x = band_range;
68+
else
69+
spsd_slc = spsd.data(idx, :);
70+
x = spsd.freq(idx);
71+
end
72+
% Integrate spectral density by frequency for each time
73+
for t = 1:nTime
74+
mspl.data(i, t) = 10 * log10(trapz(x, spsd_slc(:, t)) / reference);
75+
end
76+
end
77+
78+
mspl.freq = band.center_freq;
79+
mspl.time = spsd.time;
80+
81+
end

‎mhkit/acoustics/analysis/create_frequency_bands.m‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
function [octave_bins, band] = create_frequency_bands(octave, base, fmin, fmax)
2+
3+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4+
%
5+
% Calculate frequency bands based on the specified octave, minimum and maximum frequency limits
6+
%
7+
% Parameters
8+
% ------------
9+
% octave: double
10+
% Octave to subdivide spectral density level by
11+
% base: double, optional
12+
% Octave base. Set to 2 for the true octave band; set to base 10 for
13+
% the decidecade octave band. Default: 2
14+
% fmin: double, optional
15+
% Lower frequency band limit (lower limit of the hydrophone). Default is 10 Hz
16+
% fmax: double, optional
17+
% Upper frequency band limit (Nyquist frequency). Default is 100,000 Hz
18+
%
19+
% Returns
20+
% ---------
21+
% octave_bins: array
22+
% Array of octave bin edges
23+
% band: structure
24+
% band.center_freq : Center frequencies [Hz]
25+
% band.lower_limit : Lower frequency limits [Hz]
26+
% band.upper_limit : Upper frequency limits [Hz]
27+
%
28+
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
29+
30+
arguments (Input)
31+
octave double {mustBeFinite, mustBePositive}
32+
base double {mustBeFinite, mustBePositive} = 2
33+
fmin double {mustBeFinite, mustBePositive} = 10
34+
fmax double {mustBeFinite, mustBePositive} = 100000
35+
end
36+
37+
arguments (Output)
38+
octave_bins double
39+
band struct
40+
end
41+
42+
bandwidth = base^(1 / octave);
43+
half_bandwidth = base^(1 / (octave * 2));
44+
45+
% Calculate center frequencies
46+
log_fmin = log10(fmin);
47+
log_fmax = log10(fmax * bandwidth);
48+
step = log10(bandwidth);
49+
center_freq = 10 .^ (log_fmin : step : log_fmax);
50+
51+
lower_limit = center_freq / half_bandwidth;
52+
upper_limit = center_freq * half_bandwidth;
53+
octave_bins = [lower_limit, upper_limit(end)];
54+
55+
band = struct();
56+
band.center_freq = center_freq;
57+
band.lower_limit = lower_limit;
58+
band.upper_limit = upper_limit;
59+
60+
end

0 commit comments

Comments
 (0)