Skip to content

Commit 2fd5cd2

Browse files
hivanov-nrelhivanov-nrel
committed
acoustics module functions
1 parent c384d7e commit 2fd5cd2

24 files changed

Lines changed: 1597 additions & 9 deletions

‎mhkit/acoustics/analysis/apply_calibration.m‎

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

‎mhkit/acoustics/analysis/band_aggregate.m‎

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

‎mhkit/acoustics/analysis/band_sound_pressure_level.m‎

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

‎mhkit/acoustics/analysis/create_frequency_bands.m‎

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

‎mhkit/acoustics/analysis/decidecade_sound_pressure_level.m‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
function mspl = decidecade_sound_pressure_level(spsd, fmin, fmax)
2+
%
3+
% Calculates the sound pressure level in decidecade bands directly
4+
% from the mean square sound pressure spectral density (SPSD).
5+
%
6+
% Parameters
7+
% ----------
8+
% spsd: struct
9+
% Mean square sound pressure spectral density in [Pa^2/Hz].
10+
% fmin: int
11+
% Lower frequency band limit (lower limit of the hydrophone).
12+
% Default: 10 Hz
13+
% fmax: int
14+
% Upper frequency band limit (Nyquist frequency).
15+
% Default: 100000 Hz
16+
%
17+
% Returns
18+
% -------
19+
% mspl : struct
20+
% Sound pressure level [dB re 1 uPa] indexed by time and decidecade bands
21+
22+
23+
arguments (Input)
24+
spsd struct
25+
fmin {mustBeNumeric} = 10
26+
fmax {mustBeNumeric} = 100000
27+
end
28+
29+
arguments (Output)
30+
mspl struct
31+
end
32+
33+
fmax = fmax_warning(spsd.fs/2, fmax);
34+
35+
octave = 10;
36+
base = 10;
37+
38+
mspl = band_sound_pressure_level(spsd, octave, base, fmin, fmax);
39+
mspl.units = 'db re 1 uPa';
40+
mspl.name = 'Decidecade Sound Pressure Level';
41+
42+
end

‎mhkit/acoustics/analysis/fft_hann.m‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
function [f, spec] = fft_hann(fs, x, nfft)
2+
%
3+
% Apply Fast-Fourier-Transform to a time series using a Hanning window.
4+
%
5+
% Parameters
6+
% -----------
7+
% fs: int
8+
% Sample rate of data.
9+
% x: array
10+
% Array of timeseries data.
11+
% nfft: int
12+
% Number of elements in the FFT.
13+
%
14+
% Returns
15+
% -------
16+
% f: array
17+
% Frequency array vector.
18+
% spec: array
19+
% Frequency spectra resulting from the FFT.
20+
21+
arguments (Input)
22+
fs {mustBeNumeric}
23+
x {mustBeNumeric}
24+
nfft {mustBeNumeric}
25+
end
26+
27+
arguments (Output)
28+
f {mustBeNumeric}
29+
spec {mustBeNumeric}
30+
end
31+
32+
x0 = x-mean(x); % demean the time series
33+
N = length(x); % get length of time series
34+
w = 0.5 * (1-cos(2*pi*(0:1:N-1) / (N-1))); % generate hanning window
35+
xx = x0.*w'; % apply hamming window
36+
Fn = fs/2; % Nyquist frequency
37+
NumFFT = nfft;
38+
TempFFT = fft(xx,NumFFT); % Take FFT, padding with zeros.
39+
NumUniquePts = ceil((NumFFT)/2);
40+
spec = TempFFT(1:NumUniquePts); % FFT is symmetric, throw away second half
41+
spec = spec*2; % Multiply by 2 to take into account the fact that we
42+
% threw out second half of TempFFT above
43+
spec(1) = spec(1)/2; % Account for endpoint uniqueness
44+
spec(length(spec))= spec(length(spec))/2; % We know NumFFT is even
45+
spec = spec/length(xx); % Scale the FFT so that it is not a function of the length of x.
46+
f = (1:NumUniquePts)*2*Fn/NumFFT;
47+
48+
end

‎mhkit/acoustics/analysis/fmax_warning.m‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
function fmax_out = fmax_warning(fn, fmax)
2+
%
3+
% Checks that the maximum frequency limit isn't greater than the Nyquist frequency.
4+
%
5+
% Parameters
6+
% ----------
7+
% fn: double
8+
% The Nyquist frequency in Hz.
9+
% fmax: double
10+
% The maximum frequency limit in Hz.
11+
%
12+
% Returns
13+
% -------
14+
% fmax: double
15+
% The adjusted maximum frequency limit, ensuring it does not exceed
16+
% the Nyquist frequency.
17+
18+
19+
arguments (Input)
20+
fn {mustBeNumeric}
21+
fmax {mustBeNumeric}
22+
end
23+
24+
arguments (Output)
25+
fmax_out {mustBeNumeric}
26+
end
27+
28+
if fmax > fn
29+
warning('fmax = %.1f is greater than Nyquist frequency. Setting fmax = %.1f', fmax, fn)
30+
fmax = fn;
31+
end
32+
33+
fmax_out = fmax;
34+
35+
end

0 commit comments

Comments
 (0)