[FEATURE] add the tonality of Aures 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) 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) 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_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

tests/golden/make_tonality_sqat.m

@@ -0,0 +1,100 @@
+function make_tonality_sqat(sqatRoot)
+% MAKE_TONALITY_SQAT  Golden data for tests/test_psychoacoustics.py (tonality).
+%
+%   make_tonality_sqat(sqatRoot) runs Tonality_Aures1985 from SQAT
+%   (https://github.com/ggrecow/SQAT) on short synthetic signals and writes
+%   tonality_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 papers. 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:
+%     'tone': sum over the frequencies f of sqrt(2)*2e-5*10^(L/20)*sin(2*pi*f*t);
+%     'tone+noise': the tone plus a multisine noise with components every
+%             1 Hz from f1 to f2, phases 2*pi*mod(g*k^2, 1) with
+%             g = 0.6180339887 (k = 1..K), which scatter the group delay
+%             over the 1 s period so that every frame holds the whole band,
+%             scaled to the rms level Ln.
+%
+%   Tonality: SQAT with the free-field loudness, time_skip 0.
+%   Stored per case (struct array "cases"): the signal parameters, the sum
+%   and sum of squares of the single precision input, and the outputs K(t),
+%   time and Kmean. Also stored: matlab_version and sqat_commit.
+%
+%   Usage:  matlab -singleCompThread -batch "make_tonality_sqat('/path/to/SQAT')"
+
+addpath(sqatRoot);
+startup_SQAT([strip(sqatRoot, "right", filesep) filesep]);  % expects a trailing separator
+
+tdefs = {
+%   name            kind          f            L   f1   f2    Ln  fs     dur
+    'tone1k_60'     'tone'        1000         60  0    0     0   48000  1.0
+    'tone500_60'    'tone'        500          60  0    0     0   48000  1.0
+    'tone4k_60'     'tone'        4000         60  0    0     0   48000  1.0
+    'tone1k_30'     'tone'        1000         30  0    0     0   48000  1.0
+    'tones1k_1k2'   'tone'        [1000 1200]  60  0    0     0   48000  1.0
+    'tone1k_noise60' 'tone+noise' 1000         60  200  5000  60  48000  1.0
+    'tone1k_noise50' 'tone+noise' 1000         60  200  5000  50  48000  1.0
+};
+cases = struct([]);
+for c = 1:size(tdefs, 1)
+    p = cell2struct(tdefs(c, :), {'name','kind','f','L','f1','f2','Ln','fs','dur'}, 2);
+    t = (0:round(p.dur * p.fs)-1)' / p.fs;
+    x = zeros(size(t));
+    for f = p.f
+        x = x + sqrt(2)*2e-5*10^(p.L/20) * sin(2*pi*f*t);
+    end
+    if strcmp(p.kind, 'tone+noise')
+        K = p.f2 - p.f1 + 1;
+        nz = zeros(size(t));
+        for k = 1:K
+            nz = nz + sin(2*pi*(p.f1 + k - 1)*t + 2*pi*mod(0.6180339887*k^2, 1));
+        end
+        x = x + nz * (2e-5*10^(p.Ln/20) / sqrt(mean(nz.^2)));
+    end
+    x = double(single(x));
+    OUT = Tonality_Aures1985(x, p.fs, 0, 0, 0);
+    p.sum_x  = sum(x);
+    p.sum_x2 = sum(x.^2);
+    p.K      = OUT.InstantaneousTonality;
+    p.time   = OUT.time;
+    p.Kmean  = OUT.Kmean;
+    cases = il_append(cases, p);
+    fprintf('%-14s Kmean = %.4f t.u. (%d frames)\n', p.name, p.Kmean, numel(p.K));
+end
+il_save('tonality_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