[FEATURE] add the fluctuation strength of Osses et al. to pytta.psychoacoustics

docs/pyttadoc/index.rst

@@ -57,11 +57,17 @@ Modules:
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.psychoacoustics:
+        loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel
+        and Weber 1997), fluctuation strength (Osses et al. 2016) and
+        tonality (Aures 1985) of a sound pressure signal;
+
 .. toctree::
    classes/index
    utils/index
    apps/index
    functions
    generate
    rooms
+   psychoacoustics
    properties

docs/pyttadoc/psychoacoustics.rst

@@ -0,0 +1,7 @@
+.. currentmodule:: pytta
+
+Psychoacoustics
+===============
+
+.. automodule:: pytta.psychoacoustics
+   :members:

pytta/__init__.py

@@ -72,6 +72,11 @@
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.psychoacoustics:
+        loudness (ISO 532-1), sharpness (DIN 45692), roughness (Daniel
+        and Weber 1997), fluctuation strength (Osses et al. 2016) and
+        tonality (Aures 1985) of a sound pressure signal;
+
 
 You can find usage examples in the toolbox's examples folder.
 
@@ -108,6 +113,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import psychoacoustics
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +136,7 @@
            'generate',
            'utils',
            'iso3741',
+           'psychoacoustics',
            'rooms',
 
            # Functions

tests/golden/make_fluctuation_sqat.m

@@ -0,0 +1,94 @@
+function make_fluctuation_sqat(sqatRoot)
+% MAKE_FLUCTUATION_SQAT  Golden data for tests/test_psychoacoustics.py
+% (fluctuation strength).
+%
+%   make_fluctuation_sqat(sqatRoot) runs FluctuationStrength_Osses2016 from
+%   SQAT (https://github.com/ggrecow/SQAT) on short synthetic signals and
+%   writes fluctuation_sqat.mat next to this file. SQAT is used here only as
+%   a test oracle: that SQAT file is GPL-3.0 and pytta.psychoacoustics is
+%   written from the paper. Python rebuilds the same signals from the
+%   parameters stored here (see the helpers of the test file).
+%
+%   Signals (sound pressure in Pa), rounded to single precision before the
+%   call, as SignalObj stores float32, t = (0:n-1)'/fs:
+%     'am':   A*(1 + m*sin(2*pi*fm*t)).*sin(2*pi*fc*t),
+%             A = sqrt(2)*2e-5*10^(L/20)/sqrt(1 + m^2/2) (L: rms level);
+%     'fm':   sqrt(2)*2e-5*10^(L/20)*sin(2*pi*fc*t - (df/fm)*cos(2*pi*fm*t))
+%             (frequency fc + df*sin(2*pi*fm*t));
+%
+%   Fluctuation strength: SQAT with method 1 (time-varying), time_skip 0,
+%   struct_opt.a0_type = 'fastl2007' (the a0 curve that pytta uses).
+%   Stored per case (struct array "cases"): the signal parameters, the sum
+%   and sum of squares of the single precision input, and the outputs FS(t),
+%   time, specific FS (47 x frames) and FSmean. Also stored: matlab_version and sqat_commit.
+%
+%   Usage:  matlab -singleCompThread -batch "make_fluctuation_sqat('/path/to/SQAT')"
+
+addpath(sqatRoot);
+startup_SQAT([strip(sqatRoot, "right", filesep) filesep]);  % expects a trailing separator
+
+fdefs = {
+%   name          kind   fc    fm  m   df   L   fs     dur
+    'am1k_4hz'    'am'   1000  4   1   0    70  48000  4.0
+    'am1k_16hz'   'am'   1000  16  1   0    70  48000  4.0
+    'am1k_4hz_441' 'am'  1000  4   1   0    60  44100  3.0
+    'fm1k5_4hz'   'fm'   1500  4   0   700  70  48000  4.0
+};
+opt.a0_type = 'fastl2007';
+cases = struct([]);
+for c = 1:size(fdefs, 1)
+    p = cell2struct(fdefs(c, :), {'name','kind','fc','fm','m','df','L','fs','dur'}, 2);
+    t = (0:round(p.dur * p.fs)-1)' / p.fs;
+    if strcmp(p.kind, 'am')
+        A = sqrt(2)*2e-5*10^(p.L/20) / sqrt(1 + p.m^2/2);
+        x = A * (1 + p.m*sin(2*pi*p.fm*t)) .* sin(2*pi*p.fc*t);
+    else
+        x = sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*p.fc*t - (p.df/p.fm)*cos(2*pi*p.fm*t));
+    end
+    x = double(single(x));
+    OUT = FluctuationStrength_Osses2016(x, p.fs, 1, 0, 0, opt);
+    p.sum_x  = sum(x);
+    p.sum_x2 = sum(x.^2);
+    p.FS     = OUT.InstantaneousFluctuationStrength;
+    p.time   = OUT.time;
+    p.FSspec = OUT.InstantaneousSpecificFluctuationStrength;
+    p.FSmean = OUT.FSmean;
+    cases = il_append(cases, p);
+    fprintf('%-14s FSmean = %.4f vacil (%d frames)\n', p.name, p.FSmean, numel(p.FS));
+end
+il_save('fluctuation_sqat.mat', cases, sqatRoot);
+
+end
+
+function cases = il_append(cases, p)
+if isempty(cases)
+    cases = p;
+else
+    cases(end+1) = p;
+end
+end
+
+function il_save(name, cases, sqatRoot) %#ok<INUSL>
+matlab_version = version; %#ok<NASGU>
+sqat_commit = il_sqat_commit(sqatRoot); %#ok<NASGU>
+out = fullfile(fileparts(mfilename('fullpath')), name);
+save(out, 'cases', 'matlab_version', 'sqat_commit', '-v7');
+fprintf('wrote %s (%d cases)\n', out, numel(cases));
+end
+
+function c = il_sqat_commit(sqatRoot)
+% Commit of the SQAT checkout, read from its .git folder ('' if unknown).
+c = '';
+g = fullfile(sqatRoot, '.git');
+if ~isfolder(g), return; end
+head = strtrim(fileread(fullfile(g, 'HEAD')));
+if ~startsWith(head, 'ref: '), c = head; return; end
+ref = head(6:end);
+if isfile(fullfile(g, ref))
+    c = strtrim(fileread(fullfile(g, ref)));
+elseif isfile(fullfile(g, 'packed-refs'))
+    lines = splitlines(fileread(fullfile(g, 'packed-refs')));
+    hit = lines(endsWith(lines, [' ' ref]));
+    if ~isempty(hit), c = extractBefore(hit{1}, ' '); end
+end
+end

tests/golden/make_loudness_sqat.m

@@ -0,0 +1,122 @@
+function make_loudness_sqat(sqatRoot)
+% MAKE_LOUDNESS_SQAT  Golden data for tests/test_psychoacoustics.py (loudness).
+%
+%   make_loudness_sqat(sqatRoot) runs Loudness_ISO532_1 from SQAT
+%   (https://github.com/ggrecow/SQAT) on short synthetic signals and writes
+%   loudness_sqat.mat next to this file. Python rebuilds the same signals
+%   from the parameters stored here (see the helpers of the test file).
+%
+%   Signals (sound pressure in Pa). Every signal is rounded to single
+%   precision before the call, as SignalObj stores float32, so MATLAB and
+%   Python see the same samples:
+%     tone      a*sin(2*pi*f*t), a = sqrt(2)*2e-5*10^(L/20)
+%     multitone the 28 one-third-octave centre frequencies 25 Hz..12.5 kHz,
+%               each at L dB, Schroeder phases pi*k*(k-1)/28, summed in
+%               order k = 1..28
+%     burst     tone that is non-zero for t0 <= t < t1 only
+%     ramp      tone whose level goes linearly from L to L2 dB over the
+%               signal: a(t) = sqrt(2)*2e-5*10.^((L+(L2-L)*t/T)/20)
+%   with t = (0:n-1)'/fs and T = n/fs.
+%
+%   Stored per case (cell array "cases" of structs): the signal
+%   parameters, the call arguments (field 0 free / 1 diffuse, method 1
+%   stationary / 2 time-varying, time_skip), sum and sum of squares of the
+%   single precision input (to check that Python rebuilt the same signal),
+%   and the outputs: stationary N, LN, Nspec [1x240]; time-varying N(t), LN(t),
+%   time, N5, Nmax, and the specific loudness at every 25th output frame
+%   (Nspec_t [frames x 240], frames in Nspec_idx, 1-based) to keep the file
+%   small.
+%
+%   Also stored: matlab_version and sqat_commit (the commit of the SQAT
+%   checkout, read from its .git folder).
+%
+%   Usage:  matlab -singleCompThread -batch "make_loudness_sqat('/path/to/SQAT')"
+
+addpath(sqatRoot);
+startup_SQAT([strip(sqatRoot, "right", filesep) filesep]);  % expects a trailing separator
+
+defs = {
+%   name             kind        f     L   L2  t0   t1   fs     dur  field method skip
+    'tone1k_40'      'tone'      1000  40  0   0    0    48000  1.0  0     1      0
+    'multitone45'    'multitone' 0     45  0   0    0    48000  1.0  0     1      0.2
+    'multitone45_D'  'multitone' 0     45  0   0    0    48000  1.0  1     1      0
+    'tone250_80_441' 'tone'      250   80  0   0    0    44100  1.0  0     1      0
+    'burst1k_70'     'burst'     1000  70  0   0.1  0.2  48000  0.5  0     2      0
+    'ramp1k_40_80_D' 'ramp'      1000  40  80  0    0    48000  1.0  1     2      0.2
+    'tone1k_60_441'  'tone'      1000  60  0   0    0    44100  0.5  0     2      0
+};
+
+cases = cell(1, size(defs, 1));  % cell array: stationary and time-varying cases store different fields
+for c = 1:size(defs, 1)
+    d = defs(c, :);
+    p = cell2struct(d, {'name','kind','f','L','L2','t0','t1','fs','dur', ...
+                        'field','method','time_skip'}, 2);
+    x = il_signal(p);
+    OUT = Loudness_ISO532_1(x, p.fs, p.field, p.method, p.time_skip, 0);
+    p.sum_x  = sum(x);
+    p.sum_x2 = sum(x.^2);
+    if p.method == 1
+        p.N     = OUT.Loudness;
+        p.LN    = OUT.LoudnessLevel;
+        p.Nspec = OUT.SpecificLoudness;
+    else
+        p.N    = OUT.InstantaneousLoudness;
+        p.LN   = OUT.InstantaneousLoudnessLevel;
+        p.time = OUT.time;
+        p.N5   = OUT.N5;
+        p.Nmax = OUT.Nmax;
+        idx = 1:25:numel(OUT.time);
+        p.Nspec_idx = idx;
+        p.Nspec_t   = OUT.InstantaneousSpecificLoudness(idx, :);
+    end
+    cases{c} = p;
+    fprintf('%-16s N = %s\n', p.name, mat2str(p.N(1:min(3,end))', 6));
+end
+
+matlab_version = version;
+sqat_commit = il_sqat_commit(sqatRoot);
+out = fullfile(fileparts(mfilename('fullpath')), 'loudness_sqat.mat');
+save(out, 'cases', 'matlab_version', 'sqat_commit', '-v7');
+fprintf('wrote %s (%d cases)\n', out, numel(cases));
+end
+
+function x = il_signal(p)
+n = round(p.dur * p.fs);
+t = (0:n-1)' / p.fs;
+switch p.kind
+    case 'tone'
+        x = sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*p.f*t);
+    case 'burst'
+        x = sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*p.f*t);
+        x(~(t >= p.t0 & t < p.t1)) = 0;
+    case 'ramp'
+        T = n / p.fs;
+        x = sqrt(2)*2e-5*10.^((p.L + (p.L2-p.L)*t/T)/20) .* sin(2*pi*p.f*t);
+    case 'multitone'
+        K = 28;
+        a = sqrt(2)*2e-5*10^(p.L/20);
+        x = zeros(n, 1);
+        for k = 1:K
+            fk = 1000 * 10^(((k-1)-16)/10);
+            x = x + a * sin(2*pi*fk*t + pi*k*(k-1)/K);
+        end
+end
+x = double(single(x));
+end
+
+function c = il_sqat_commit(sqatRoot)
+% Commit of the SQAT checkout, read from its .git folder ('' if unknown).
+c = '';
+g = fullfile(sqatRoot, '.git');
+if ~isfolder(g), return; end
+head = strtrim(fileread(fullfile(g, 'HEAD')));
+if ~startsWith(head, 'ref: '), c = head; return; end
+ref = head(6:end);
+if isfile(fullfile(g, ref))
+    c = strtrim(fileread(fullfile(g, ref)));
+elseif isfile(fullfile(g, 'packed-refs'))
+    lines = splitlines(fileread(fullfile(g, 'packed-refs')));
+    hit = lines(endsWith(lines, [' ' ref]));
+    if ~isempty(hit), c = extractBefore(hit{1}, ' '); end
+end
+end

tests/golden/make_roughness_sqat.m

@@ -0,0 +1,82 @@
+function make_roughness_sqat(sqatRoot)
+% MAKE_ROUGHNESS_SQAT  Golden data for tests/test_psychoacoustics.py (roughness).
+%
+%   make_roughness_sqat(sqatRoot) runs Roughness_Daniel1997 from SQAT
+%   (https://github.com/ggrecow/SQAT) on short synthetic signals and writes
+%   roughness_sqat.mat next to this file. Python rebuilds the same signals
+%   from the parameters stored here (see the helpers of the test file).
+%
+%   Signal (sound pressure in Pa), rounded to single precision before the
+%   call, as SignalObj stores float32:
+%     x = A*(1 + m*sin(2*pi*fm*t)).*sin(2*pi*fc*t),  t = (0:n-1)'/fs,
+%     A = sqrt(2)*2e-5*10^(L/20)/sqrt(1 + m^2/2),
+%   so that L is the level of the rms of the modulated signal (m = 0 gives
+%   a pure tone).
+%
+%   Stored per case (struct array "cases"): the signal parameters,
+%   time_skip, sum and sum of squares of the single precision input (to
+%   check that Python rebuilt the same signal), and the outputs R(t)
+%   (InstantaneousRoughness), time, the specific roughness
+%   (InstantaneousSpecificRoughness, 47 x frames) and Rmean.
+%
+%   Also stored: matlab_version and sqat_commit (the commit of the SQAT
+%   checkout, read from its .git folder).
+%
+%   Usage:  matlab -singleCompThread -batch "make_roughness_sqat('/path/to/SQAT')"
+
+addpath(sqatRoot);
+startup_SQAT([strip(sqatRoot, "right", filesep) filesep]);  % expects a trailing separator
+
+defs = {
+%   name              fc    fm   m    L   fs     dur  skip
+    'am1k_70_ref'     1000  70   1    60  48000  2.0  0.5
+    'am1k_70_441'     1000  70   1    60  44100  1.0  0
+    'am2k_50_m05_32k' 2000  50   0.5  70  32000  1.0  0
+    'tone1k_60'       1000  0    0    60  48000  1.0  0
+};
+
+cases = struct([]);
+for c = 1:size(defs, 1)
+    p = cell2struct(defs(c, :), {'name','fc','fm','m','L','fs','dur','time_skip'}, 2);
+    n = round(p.dur * p.fs);
+    t = (0:n-1)' / p.fs;
+    A = sqrt(2)*2e-5*10^(p.L/20) / sqrt(1 + p.m^2/2);
+    x = double(single(A * (1 + p.m*sin(2*pi*p.fm*t)) .* sin(2*pi*p.fc*t)));
+    OUT = Roughness_Daniel1997(x, p.fs, p.time_skip, 0);
+    p.sum_x  = sum(x);
+    p.sum_x2 = sum(x.^2);
+    p.R      = OUT.InstantaneousRoughness;
+    p.time   = OUT.time;
+    p.Rspec  = OUT.InstantaneousSpecificRoughness;
+    p.Rmean  = OUT.Rmean;
+    if isempty(cases)
+        cases = p;
+    else
+        cases(end+1) = p; %#ok<AGROW>
+    end
+    fprintf('%-16s Rmean = %.6f asper (%d frames)\n', p.name, p.Rmean, numel(p.R));
+end
+
+matlab_version = version;
+sqat_commit = il_sqat_commit(sqatRoot);
+out = fullfile(fileparts(mfilename('fullpath')), 'roughness_sqat.mat');
+save(out, 'cases', 'matlab_version', 'sqat_commit', '-v7');
+fprintf('wrote %s (%d cases)\n', out, numel(cases));
+end
+
+function c = il_sqat_commit(sqatRoot)
+% Commit of the SQAT checkout, read from its .git folder ('' if unknown).
+c = '';
+g = fullfile(sqatRoot, '.git');
+if ~isfolder(g), return; end
+head = strtrim(fileread(fullfile(g, 'HEAD')));
+if ~startsWith(head, 'ref: '), c = head; return; end
+ref = head(6:end);
+if isfile(fullfile(g, ref))
+    c = strtrim(fileread(fullfile(g, ref)));
+elseif isfile(fullfile(g, 'packed-refs'))
+    lines = splitlines(fileread(fullfile(g, 'packed-refs')));
+    hit = lines(endsWith(lines, [' ' ref]));
+    if ~isempty(hit), c = extractBefore(hit{1}, ' '); end
+end
+end