Add 2f model source code

tools/install_2fmodel.md

@@ -0,0 +1 @@
+../versa/sequence_metrics/2f_model/README.md

versa/sequence_metrics/2f_model/CB/PQCB.m

@@ -0,0 +1,99 @@
+function [Nc, fc, fl, fu, dz] = PQCB (Version)
+% Critical band parameters for the FFT model
+
+% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:32:57 $
+
+B = inline ('7 * asinh (f / 650)');
+BI = inline ('650 * sinh (z / 7)');
+
+fL = 80;
+fU = 18000;
+
+% Critical bands - set up the tables
+if (strcmp (Version, 'Basic'))
+    dz = 1/4;
+elseif (strcmp (Version, 'Advanced'))
+    dz = 1/2;
+else
+    error ('PQCB: Invalid version');
+end
+
+zL = B(fL);
+zU = B(fU);
+Nc = ceil((zU - zL) / dz);
+zl = zL + (0:Nc-1) * dz;
+zu = min (zL + (1:Nc) * dz, zU);
+zc = 0.5 * (zl + zu);
+
+fl = BI (zl);
+fc = BI (zc);
+fu = BI (zu);
+
+if (strcmp (Version, 'Basic'))
+    fl = [  80.000,   103.445,   127.023,   150.762,   174.694, ...
+           198.849,   223.257,   247.950,   272.959,   298.317, ...
+           324.055,   350.207,   376.805,   403.884,   431.478, ...
+           459.622,   488.353,   517.707,   547.721,   578.434, ...
+           609.885,   642.114,   675.161,   709.071,   743.884, ...
+           779.647,   816.404,   854.203,   893.091,   933.119, ...
+           974.336,  1016.797,  1060.555,  1105.666,  1152.187, ...
+          1200.178,  1249.700,  1300.816,  1353.592,  1408.094, ...
+          1464.392,  1522.559,  1582.668,  1644.795,  1709.021, ...
+          1775.427,  1844.098,  1915.121,  1988.587,  2064.590, ...
+          2143.227,  2224.597,  2308.806,  2395.959,  2486.169, ...
+          2579.551,  2676.223,  2776.309,  2879.937,  2987.238, ...
+          3098.350,  3213.415,  3332.579,  3455.993,  3583.817, ...
+          3716.212,  3853.817,  3995.399,  4142.547,  4294.979, ...
+          4452.890,  4616.482,  4785.962,  4961.548,  5143.463, ...
+          5331.939,  5527.217,  5729.545,  5939.183,  6156.396, ...
+          6381.463,  6614.671,  6856.316,  7106.708,  7366.166, ...
+          7635.020,  7913.614,  8202.302,  8501.454,  8811.450, ...
+          9132.688,  9465.574,  9810.536, 10168.013, 10538.460, ...
+         10922.351, 11320.175, 11732.438, 12159.670, 12602.412, ...
+         13061.229, 13536.710, 14029.458, 14540.103, 15069.295, ...
+         15617.710, 16186.049, 16775.035, 17385.420 ];
+    fc = [  91.708,   115.216,   138.870,   162.702,   186.742, ...
+           211.019,   235.566,   260.413,   285.593,   311.136, ...
+           337.077,   363.448,   390.282,   417.614,   445.479, ...
+           473.912,   502.950,   532.629,   562.988,   594.065, ...
+           625.899,   658.533,   692.006,   726.362,   761.644, ...
+           797.898,   835.170,   873.508,   912.959,   953.576, ...
+           995.408,  1038.511,  1082.938,  1128.746,  1175.995, ...
+          1224.744,  1275.055,  1326.992,  1380.623,  1436.014, ...
+          1493.237,  1552.366,  1613.474,  1676.641,  1741.946, ...
+          1809.474,  1879.310,  1951.543,  2026.266,  2103.573, ...
+          2183.564,  2266.340,  2352.008,  2440.675,  2532.456, ...
+          2627.468,  2725.832,  2827.672,  2933.120,  3042.309, ...
+          3155.379,  3272.475,  3393.745,  3519.344,  3649.432, ...
+          3784.176,  3923.748,  4068.324,  4218.090,  4373.237, ...
+          4533.963,  4700.473,  4872.978,  5051.700,  5236.866, ...
+          5428.712,  5627.484,  5833.434,  6046.825,  6267.931, ...
+          6497.031,  6734.420,  6980.399,  7235.284,  7499.397, ...
+          7773.077,  8056.673,  8350.547,  8655.072,  8970.639, ...
+          9297.648,  9636.520,  9987.683, 10351.586, 10728.695, ...
+         11119.490, 11524.470, 11944.149, 12379.066, 12829.775, ...
+         13294.850, 13780.887, 14282.503, 14802.338, 15341.057, ...
+         15899.345, 16477.914, 17077.504, 17690.045 ];
+    fu = [ 103.445,   127.023,   150.762,   174.694,   198.849, ...
+           223.257,   247.950,   272.959,   298.317,   324.055, ...
+           350.207,   376.805,   403.884,   431.478,   459.622, ...
+           488.353,   517.707,   547.721,   578.434,   609.885, ...
+           642.114,   675.161,   709.071,   743.884,   779.647, ...
+           816.404,   854.203,   893.091,   933.113,   974.336, ...
+          1016.797,  1060.555,  1105.666,  1152.187,  1200.178, ...
+          1249.700,  1300.816,  1353.592,  1408.094,  1464.392, ...
+          1522.559,  1582.668,  1644.795,  1709.021,  1775.427, ...
+          1844.098,  1915.121,  1988.587,  2064.590,  2143.227, ...
+          2224.597,  2308.806,  2395.959,  2486.169,  2579.551, ...
+          2676.223,  2776.309,  2879.937,  2987.238,  3098.350, ...
+          3213.415,  3332.579,  3455.993,  3583.817,  3716.212, ...
+          3853.348,  3995.399,  4142.547,  4294.979,  4452.890, ...
+          4643.482,  4785.962,  4961.548,  5143.463,  5331.939, ...
+          5527.217,  5729.545,  5939.183,  6156.396,  6381.463, ...
+          6614.671,  6856.316,  7106.708,  7366.166,  7635.020, ...
+          7913.614,  8202.302,  8501.454,  8811.450,  9132.688, ...
+          9465.574,  9810.536, 10168.013, 10538.460, 10922.351, ...
+         11320.175, 11732.438, 12159.670, 12602.412, 13061.229, ...
+         13536.710, 14029.458, 14540.103, 15069.295, 15617.710, ...
+         16186.049, 16775.035, 17385.420, 18000.000 ];
+end

versa/sequence_metrics/2f_model/CB/PQDFTFrame.m

@@ -0,0 +1,60 @@
+function X2 = PQDFTFrame (x)
+% Calculate the DFT of a frame of data (NF values), returning the
+% squared-magnitude DFT vector (NF/2 + 1 values)
+
+% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:32:57 $
+
+persistent hw
+
+NF = 2048;      % Frame size (samples)
+
+if (isempty (hw))
+    Amax = 32768;
+    fc = 1019.5;
+    Fs = 48000;
+    Lp = 92;
+    % Set up the window (including all gains)
+    GL = PQ_GL (NF, Amax, fc/Fs, Lp);
+    hw = GL * PQHannWin (NF);
+end
+
+% Window the data
+xw = hw .* x;
+
+% DFT (output is real followed by imaginary)
+X = PQRFFT (xw, NF, 1);
+
+% Squared magnitude
+X2 = PQRFFTMSq (X, NF);
+
+%----------------------------------------
+function GL = PQ_GL (NF, Amax, fcN, Lp)
+% Scaled Hann window, including loudness scaling
+
+% Calculate the gain for the Hann Window
+%  - level Lp (SPL) corresponds to a sine with normalized frequency
+%    fcN and a peak value of Amax
+
+W = NF - 1;
+gp = PQ_gp (fcN, NF, W);
+GL = 10^(Lp / 20) / (gp * Amax/4 * W);
+
+%----------
+function gp = PQ_gp (fcN, NF, W)
+% Calculate the peak factor. The signal is a sinusoid windowed with
+% a Hann window. The sinusoid frequency falls between DFT bins. The
+% peak of the frequency response (on a continuous frequency scale) falls
+% between DFT bins. The largest DFT bin value is the peak factor times
+% the peak of the continuous response.
+%  fcN - Normalized sinusoid frequency (0-1)
+%  NF  - Frame (DFT) length samples
+%  NW  - Window length samples
+
+% Distance to the nearest DFT bin
+df = 1 / NF;
+k = floor (fcN / df);
+dfN = min ((k+1) * df - fcN, fcN - k * df);
+
+dfW = dfN * W;
+gp = sin(pi * dfW) / (pi * dfW * (1 - dfW^2));
+

versa/sequence_metrics/2f_model/CB/PQeval.m

@@ -0,0 +1,111 @@
+function [MOVI, Fmem] = PQeval (xR, xT, Fmem)
+% PEAQ - Process one frame with the FFT model
+
+% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:32:58 $
+
+NF = 2048;
+Version = 'Basic';
+
+% Windowed DFT
+X2(1,:) = PQDFTFrame (xR);
+X2(2,:) = PQDFTFrame (xT);
+
+% Critical band grouping and frequency spreading
+[EbN, Es] = PQ_excitCB (X2);
+
+% Time domain smoothing => "Excitation patterns"
+[Ehs(1,:), Fmem.TDS.Ef(1,:)] = PQ_timeSpread (Es(1,:), Fmem.TDS.Ef(1,:));
+[Ehs(2,:), Fmem.TDS.Ef(2,:)] = PQ_timeSpread (Es(2,:), Fmem.TDS.Ef(2,:));
+
+% Level and pattern adaptation => "Spectrally adapted patterns"
+[EP, Fmem.Adap] = PQadapt (Ehs, Fmem.Adap, Version, 'FFT');
+
+% Modulation patterns
+[M, ERavg, Fmem.Env] = PQmodPatt (Es, Fmem.Env);
+% Loudness
+MOVI.Loud.NRef  = PQloud (Ehs(1,:), Version, 'FFT');
+MOVI.Loud.NTest = PQloud (Ehs(2,:), Version, 'FFT');
+
+% Modulation differences
+MOVI.MDiff = PQmovModDiffB (M, ERavg);
+
+% Noise Loudness
+MOVI.NLoud.NL = PQmovNLoudB (M, EP);
+
+% Bandwidth
+MOVI.BW = PQmovBW (X2);
+
+% Noise-to-mask ratios
+MOVI.NMR = PQmovNMRB (EbN, Ehs(1,:));
+
+% Probability of detection
+MOVI.PD = PQmovPD (Ehs);
+
+% Error harmonic structure
+MOVI.EHS.EHS = PQmovEHS (xR, xT, X2);
+
+%--------------------
+function [EbN, Es] = PQ_excitCB (X2)
+
+persistent W2 EIN
+
+NF = 2048;
+Version = 'Basic';
+if (isempty (W2))
+    Fs = 48000;
+    f = linspace (0, Fs/2, NF/2+1);
+    W2 = PQWOME (f);
+    [Nc, fc] = PQCB (Version);
+    EIN = PQIntNoise (fc);
+end
+
+% Allocate storage
+XwN2 = zeros (1, NF/2+1);
+
+% Outer and middle ear filtering
+Xw2(1,:) = W2 .* X2(1,1:NF/2+1);
+Xw2(2,:) = W2 .* X2(2,1:NF/2+1);
+
+% Form the difference magnitude signal
+for (k = 0:NF/2)
+    XwN2(k+1) = (Xw2(1,k+1) - 2 * sqrt (Xw2(1,k+1) * Xw2(2,k+1)) ...
+               + Xw2(2,k+1));
+end
+
+% Group into partial critical bands
+Eb(1,:) = PQgroupCB (Xw2(1,:), Version);
+Eb(2,:) = PQgroupCB (Xw2(2,:), Version);
+EbN     = PQgroupCB (XwN2, Version);
+
+% Add the internal noise term => "Pitch patterns"
+E(1,:) = Eb(1,:) + EIN;
+E(2,:) = Eb(2,:) + EIN;
+
+% Critical band spreading => "Unsmeared excitation patterns"
+Es(1,:) = PQspreadCB (E(1,:), Version);
+Es(2,:) = PQspreadCB (E(2,:), Version);
+
+%--------------------
+function [Ehs, Ef] = PQ_timeSpread (Es, Ef)
+
+persistent Nc a b
+
+if (isempty (Nc))
+    [Nc, fc] = PQCB ('Basic');
+    Fs = 48000;
+    NF = 2048;
+    Nadv = NF / 2;
+    Fss = Fs / Nadv;
+    t100 = 0.030;
+    tmin = 0.008;
+    [a, b] = PQtConst (t100, tmin, fc, Fss);
+end
+
+% Allocate storage
+Ehs = zeros (1, Nc);
+
+% Time domain smoothing
+for (m = 0:Nc-1)
+    Ef(m+1) = a(m+1) * Ef(m+1) + b(m+1) * Es(m+1);
+    Ehs(m+1) = max(Ef(m+1), Es(m+1));
+end

versa/sequence_metrics/2f_model/CB/PQgroupCB.m

@@ -0,0 +1,63 @@
+function Eb = PQgroupCB (X2, Ver)
+% Group a DFT energy vector into critical bands
+% X2 - Squared-magnitude vector (DFT bins)
+% Eb - Excitation vector (fractional critical bands)
+
+% P. Kabal $Revision: 1.2 $  $Date: 2004/02/05 04:25:46 $
+
+persistent Nc kl ku Ul Uu Version
+
+Emin = 1e-12;
+
+if (~ strcmp (Ver, Version))
+    Version = Ver;
+    % Set up the DFT bin to critical band mapping
+    NF = 2048;
+    Fs = 48000;
+    [Nc, kl, ku, Ul, Uu] = PQ_CBMapping (NF, Fs, Version);
+end
+
+% Allocate storage
+Eb = zeros (1, Nc);
+
+% Compute the excitation in each band
+for (i = 0:Nc-1)
+    Ea = Ul(i+1) * X2(kl(i+1)+1);       % First bin
+    for (k = (kl(i+1)+1):(ku(i+1)-1))
+        Ea = Ea + X2(k+1);              % Middle bins
+    end
+    Ea = Ea + Uu(i+1) * X2(ku(i+1)+1);  % Last bin
+    Eb(i+1) = max(Ea, Emin);
+end
+
+%---------------------------------------
+function [Nc, kl, ku, Ul, Uu] = PQ_CBMapping (NF, Fs, Version)
+
+[Nc, fc, fl, fu] = PQCB (Version);
+
+% Fill in the DFT bin to critical band mappings
+df = Fs / NF;
+for (i = 0:Nc-1)
+    fli = fl(i+1);
+    fui = fu(i+1);
+    for (k = 0:NF/2)
+        if ((k+0.5)*df > fli)
+            kl(i+1) = k;        % First bin in band i
+            Ul(i+1) = (min(fui, (k+0.5)*df) ...
+                       - max(fli, (k-0.5)*df)) / df;
+            break;
+        end
+    end
+    for (k = NF/2:-1:0)
+        if ((k-0.5)*df < fui)
+            ku(i+1) = k;        % Last bin in band i
+            if (kl(i+1) == ku(i+1))
+                Uu(i+1) = 0;       % Single bin in band
+            else
+                Uu(i+1) = (min(fui, (k+0.5)*df) ...
+                           - max(fli, (k-0.5)*df)) / df;
+            end
+            break;
+        end
+    end
+end