removing test
This commit is contained in:
Moliner Eloi
parent
ab2540ec3e
commit
0fe16d77be
@ -1,5 +0,0 @@
|
||||
function [ODG]=PEMOQ(ref,test)
|
||||
[refa, fs]=audioread(ref);
|
||||
[testa, fs]=audioread(test);
|
||||
[PSM, PSMt, ODG, PSM_inst] = audioqual(refa, testa, fs);
|
||||
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@ -1,128 +0,0 @@
|
||||
/*=================================================================
|
||||
*
|
||||
* toeplitzC.C Sample .MEX file corresponding to toeplitz.m
|
||||
* Solves simple 3 body orbit problem
|
||||
*
|
||||
* The calling syntax is:
|
||||
*
|
||||
* [yp] = yprime(t, y)
|
||||
* TOEPLITZ(C,R) is a non-symmetric Toeplitz matrix having C as its
|
||||
* first column and R as its first row.
|
||||
*
|
||||
* TOEPLITZ(R) is a symmetric Toeplitz matrix for real R.
|
||||
* For a complex vector R with a real first element, T = toeplitz(r)
|
||||
* returns the Hermitian Toeplitz matrix formed from R. When the
|
||||
* first element of R is not real, the resulting matrix is Hermitian
|
||||
* off the main diagonal, i.e., T_{i,j} = conj(T_{j,i}) for i ~= j.
|
||||
*
|
||||
* You may also want to look at the corresponding M-code, yprime.m.
|
||||
*
|
||||
* This is a MEX-file for MATLAB.
|
||||
* Copyright 1984-2006 The MathWorks, Inc.
|
||||
*
|
||||
*=================================================================*/
|
||||
/* $Revision: 1.10.6.4 $ */
|
||||
#include <math.h>
|
||||
#include "mex.h"
|
||||
|
||||
/* Input Arguments */
|
||||
|
||||
#define C_IN prhs[0]
|
||||
#define R_IN prhs[1]
|
||||
|
||||
|
||||
/* Output Arguments */
|
||||
|
||||
#define T_OUT plhs[0]
|
||||
|
||||
|
||||
static void toeplitzC(
|
||||
double t[],
|
||||
double c[],
|
||||
double r[],
|
||||
int m, int n
|
||||
)
|
||||
{
|
||||
int i,j,m0;
|
||||
|
||||
for (j=0;j<n;j++){
|
||||
m0 = j>m?m:j;
|
||||
for (i=0;i<=m0;i++)
|
||||
t[j*m+i] = r[j-i];
|
||||
for (i=j+1;i<m;i++)
|
||||
t[j*m+i] = c[i-j];
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void mexFunction( int nlhs, mxArray *plhs[],
|
||||
int nrhs, const mxArray*prhs[] )
|
||||
|
||||
{
|
||||
double *tr,*ti;
|
||||
double *cr,*ci,*rr,*ri;
|
||||
mwSize tm,tn;
|
||||
|
||||
/* Check for proper number of arguments */
|
||||
|
||||
if (nrhs > 2) {
|
||||
mexErrMsgTxt("More than 2 input arguments.");
|
||||
} else if (nrhs == 0) {
|
||||
mexErrMsgTxt("1 or 2 input arguments required.");
|
||||
} else if (nrhs == 1) {
|
||||
mexErrMsgTxt("1 input argument: not implemented (yet), use toeplitz.");
|
||||
} else if (nlhs > 1) {
|
||||
mexErrMsgTxt("Too many output arguments.");
|
||||
}
|
||||
|
||||
if (nrhs == 1) {
|
||||
mexErrMsgTxt("Not implemented (yet). Please use 2 input arguments or toeplitz.m with one input argument.");
|
||||
}
|
||||
|
||||
/* Check the dimensions of Y. Y can be 4 X 1 or 1 X 4. */
|
||||
|
||||
tm = mxGetNumberOfElements(C_IN);
|
||||
tn = mxGetNumberOfElements(R_IN);
|
||||
|
||||
/* Create a matrix for the return argument */
|
||||
if (mxIsComplex(C_IN) || mxIsComplex(R_IN))
|
||||
T_OUT = mxCreateDoubleMatrix(tm, tn, mxCOMPLEX);
|
||||
else
|
||||
T_OUT = mxCreateDoubleMatrix(tm, tn, mxREAL);
|
||||
|
||||
/* Assign pointers to the various parameters */
|
||||
tr = mxGetPr(T_OUT);
|
||||
|
||||
cr = mxGetPr(C_IN);
|
||||
rr = mxGetPr(R_IN);
|
||||
|
||||
/* Do the actual computations in a subroutine */
|
||||
toeplitzC(tr,cr,rr,tm,tn);
|
||||
|
||||
/* Imaginary part */
|
||||
if (mxIsComplex(C_IN) || mxIsComplex(R_IN)){
|
||||
/*if (!mxIsComplex(C_IN)){
|
||||
mexErrMsgTxt("Not implemented (yet). A");
|
||||
}
|
||||
else*/
|
||||
if (mxIsComplex(C_IN))
|
||||
ci = mxGetPi(C_IN);
|
||||
else
|
||||
ci = mxGetPr(mxCreateDoubleMatrix(tm, 1, mxREAL));
|
||||
/*if (!mxIsComplex(R_IN)){
|
||||
mexErrMsgTxt("Not implemented (yet). B");
|
||||
}
|
||||
else*/
|
||||
if (mxIsComplex(R_IN))
|
||||
ri = mxGetPi(R_IN);
|
||||
else
|
||||
ri = mxGetPr(mxCreateDoubleMatrix(tn, 1, mxREAL));
|
||||
ti = mxGetPi(T_OUT);
|
||||
toeplitzC(ti, ci, ri, tm, tn);
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -1,99 +0,0 @@
|
||||
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
|
||||
@ -1,60 +0,0 @@
|
||||
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));
|
||||
|
||||
@ -1,112 +0,0 @@
|
||||
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
|
||||
@ -1,63 +0,0 @@
|
||||
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
|
||||
@ -1,66 +0,0 @@
|
||||
function Es = PQspreadCB (E, Ver)
|
||||
% Spread an excitation vector (pitch pattern) - FFT model
|
||||
% Both E and Es are powers
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:32:58 $
|
||||
|
||||
persistent Bs Version
|
||||
|
||||
if (~ strcmp (Ver, Version))
|
||||
Version = Ver;
|
||||
Nc = length (E);
|
||||
Bs = PQ_SpreadCB (ones(1,Nc), ones(1,Nc), Version);
|
||||
end
|
||||
|
||||
Es = PQ_SpreadCB (E, Bs, Version);
|
||||
|
||||
%-------------------------
|
||||
function Es = PQ_SpreadCB (E, Bs, Ver);
|
||||
|
||||
persistent Nc dz fc aL aUC Version
|
||||
|
||||
% Power law for addition of spreading
|
||||
e = 0.4;
|
||||
|
||||
if (~ strcmp (Ver, Version))
|
||||
Version = Ver;
|
||||
[Nc, fc, fl, fu, dz] = PQCB (Version);
|
||||
end
|
||||
|
||||
% Allocate storage
|
||||
aUCEe = zeros (1, Nc);
|
||||
Ene = zeros (1, Nc);
|
||||
Es = zeros (1, Nc);
|
||||
|
||||
% Calculate energy dependent terms
|
||||
aL = 10^(-2.7 * dz);
|
||||
for (m = 0:Nc-1)
|
||||
aUC = 10^((-2.4 - 23 / fc(m+1)) * dz);
|
||||
aUCE = aUC * E(m+1)^(0.2 * dz);
|
||||
gIL = (1 - aL^(m+1)) / (1 - aL);
|
||||
gIU = (1 - aUCE^(Nc-m)) / (1 - aUCE);
|
||||
En = E(m+1) / (gIL + gIU - 1);
|
||||
aUCEe(m+1) = aUCE^e;
|
||||
Ene(m+1) = En^e;
|
||||
end
|
||||
|
||||
% Lower spreading
|
||||
Es(Nc-1+1) = Ene(Nc-1+1);
|
||||
aLe = aL^e;
|
||||
for (m = Nc-2:-1:0)
|
||||
Es(m+1) = aLe * Es(m+1+1) + Ene(m+1);
|
||||
end
|
||||
|
||||
% Upper spreading i > m
|
||||
for (m = 0:Nc-2)
|
||||
r = Ene(m+1);
|
||||
a = aUCEe(m+1);
|
||||
for (i = m+1:Nc-1)
|
||||
r = r * a;
|
||||
Es(i+1) = Es(i+1) + r;
|
||||
end
|
||||
end
|
||||
|
||||
for (i = 0:Nc-1)
|
||||
Es(i+1) = (Es(i+1))^(1/e) / Bs(i+1);
|
||||
end
|
||||
@ -1,263 +0,0 @@
|
||||
function MOV = PQavgMOVB (MOVC, Nchan, Nwup)
|
||||
% Time average MOV precursors
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $
|
||||
|
||||
Fs = 48000;
|
||||
NF = 2048;
|
||||
Nadv = NF / 2;
|
||||
Fss = Fs / Nadv;
|
||||
tdel = 0.5;
|
||||
tex = 0.050;
|
||||
|
||||
% BandwidthRefB, BandwidthTestB
|
||||
[MOV(0+1), MOV(1+1)] = PQ_avgBW (MOVC.BW);
|
||||
|
||||
% Total NMRB, RelDistFramesB
|
||||
[MOV(2+1), MOV(10+1)] = PQ_avgNMRB (MOVC.NMR);
|
||||
|
||||
% WinModDiff1B, AvgModDiff1B, AvgModDiff2B
|
||||
N500ms = ceil (tdel * Fss);
|
||||
Ndel = max (0, N500ms - Nwup);
|
||||
[MOV(3+1), MOV(6+1), MOV(7+1)] = PQ_avgModDiffB (Ndel, MOVC.MDiff);
|
||||
|
||||
% RmsNoiseLoudB
|
||||
N50ms = ceil (tex * Fss);
|
||||
Nloud = PQloudTest (MOVC.Loud);
|
||||
Ndel = max (Nloud + N50ms, Ndel);
|
||||
MOV(8+1) = PQ_avgNLoudB (Ndel, MOVC.NLoud);
|
||||
|
||||
% ADBB, MFPDB
|
||||
[MOV(4+1), MOV(9+1)] = PQ_avgPD (MOVC.PD);
|
||||
|
||||
% EHSB
|
||||
MOV(5+1) = PQ_avgEHS (MOVC.EHS);
|
||||
|
||||
%-----------------------------------------
|
||||
function EHSB = PQ_avgEHS (EHS)
|
||||
|
||||
[Nchan, Np] = size (EHS.EHS);
|
||||
|
||||
s = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
s = s + PQ_LinPosAvg (EHS.EHS(j+1,:));
|
||||
end
|
||||
EHSB = 1000 * s / Nchan;
|
||||
|
||||
|
||||
%-----------------------------------------
|
||||
function [ADBB, MFPDB] = PQ_avgPD (PD)
|
||||
|
||||
global PQopt
|
||||
|
||||
c0 = 0.9;
|
||||
if (isempty (PQopt))
|
||||
c1 = 1;
|
||||
else
|
||||
c1 = PQopt.PDfactor;
|
||||
end
|
||||
|
||||
N = length (PD.Pc);
|
||||
Phc = 0;
|
||||
Pcmax = 0;
|
||||
Qsum = 0;
|
||||
nd = 0;
|
||||
for (i = 0:N-1)
|
||||
Phc = c0 * Phc + (1 - c0) * PD.Pc(i+1);
|
||||
Pcmax = max (Pcmax * c1, Phc);
|
||||
|
||||
if (PD.Pc(i+1) > 0.5)
|
||||
nd = nd + 1;
|
||||
Qsum = Qsum + PD.Qc(i+1);
|
||||
end
|
||||
end
|
||||
|
||||
if (nd == 0)
|
||||
ADBB = 0;
|
||||
elseif (Qsum > 0)
|
||||
ADBB = log10 (Qsum / nd);
|
||||
else
|
||||
ADBB = -0.5;
|
||||
end
|
||||
|
||||
MFPDB = Pcmax;
|
||||
|
||||
%-----------------------------------------
|
||||
function [TotalNMRB, RelDistFramesB] = PQ_avgNMRB (NMR)
|
||||
|
||||
[Nchan, Np] = size (NMR.NMRavg);
|
||||
Thr = 10^(1.5 / 10);
|
||||
|
||||
s = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
s = s + 10 * log10 (PQ_LinAvg (NMR.NMRavg(j+1,:)));
|
||||
end
|
||||
TotalNMRB = s / Nchan;
|
||||
|
||||
s = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
s = s + PQ_FractThr (Thr, NMR.NMRmax(j+1,:));
|
||||
end
|
||||
RelDistFramesB = s / Nchan;
|
||||
|
||||
%-----------------------------------------
|
||||
function [BandwidthRefB, BandwidthTestB] = PQ_avgBW (BW)
|
||||
|
||||
[Nchan, Np] = size (BW.BWRef);
|
||||
|
||||
sR = 0;
|
||||
sT = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
sR = sR + PQ_LinPosAvg (BW.BWRef(j+1,:));
|
||||
sT = sT + PQ_LinPosAvg (BW.BWTest(j+1,:));
|
||||
end
|
||||
BandwidthRefB = sR / Nchan;
|
||||
BandwidthTestB = sT / Nchan;
|
||||
|
||||
%-----------------------------------------
|
||||
function [WinModDiff1B, AvgModDiff1B, AvgModDiff2B] = PQ_avgModDiffB (Ndel, MDiff)
|
||||
|
||||
NF = 2048;
|
||||
Nadv = NF / 2;
|
||||
Fs = 48000;
|
||||
|
||||
Fss = Fs / Nadv;
|
||||
tavg = 0.1;
|
||||
|
||||
[Nchan, Np] = size (MDiff.Mt1B);
|
||||
|
||||
% Sliding window average - delayed average
|
||||
L = floor (tavg * Fss); % 100 ms sliding window length
|
||||
s = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
s = s + PQ_WinAvg (L, MDiff.Mt1B(j+1,Ndel+1:Np-1+1));
|
||||
end
|
||||
WinModDiff1B = s / Nchan;
|
||||
|
||||
% Weighted linear average - delayed average
|
||||
s = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
s = s + PQ_WtAvg (MDiff.Mt1B(j+1,Ndel+1:Np-1+1), MDiff.Wt(j+1,Ndel+1:Np-1+1));
|
||||
end
|
||||
AvgModDiff1B = s / Nchan;
|
||||
|
||||
% Weighted linear average - delayed average
|
||||
s = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
s = s + PQ_WtAvg (MDiff.Mt2B(j+1,Ndel+1:Np-1+1), MDiff.Wt(j+1,Ndel+1:Np-1+1));
|
||||
end
|
||||
AvgModDiff2B = s / Nchan;
|
||||
|
||||
%-----------------------------------------
|
||||
function RmsNoiseLoudB = PQ_avgNLoudB (Ndel, NLoud)
|
||||
|
||||
[Nchan, Np] = size (NLoud.NL);
|
||||
|
||||
% RMS average - delayed average and loudness threshold
|
||||
s = 0;
|
||||
for (j = 0:Nchan-1)
|
||||
s = s + PQ_RMSAvg (NLoud.NL(j+1,Ndel+1:Np-1+1));
|
||||
end
|
||||
RmsNoiseLoudB = s / Nchan;
|
||||
|
||||
%-----------------------------------
|
||||
% Average values values, omitting values which are negative
|
||||
function s = PQ_LinPosAvg (x)
|
||||
|
||||
N = length(x);
|
||||
|
||||
Nv = 0;
|
||||
s = 0;
|
||||
for (i = 0:N-1)
|
||||
if (x(i+1) >= 0)
|
||||
s = s + x(i+1);
|
||||
Nv = Nv + 1;
|
||||
end
|
||||
end
|
||||
|
||||
if (Nv > 0)
|
||||
s = s / Nv;
|
||||
end
|
||||
|
||||
%----------
|
||||
% Fraction of values above a threshold
|
||||
function Fd = PQ_FractThr (Thr, x)
|
||||
|
||||
N = length (x);
|
||||
|
||||
Nv = 0;
|
||||
for (i = 0:N-1)
|
||||
if (x(i+1) > Thr)
|
||||
Nv = Nv + 1;
|
||||
end
|
||||
end
|
||||
|
||||
if (N > 0)
|
||||
Fd = Nv / N;
|
||||
else
|
||||
Fd = 0;
|
||||
end
|
||||
|
||||
%-----------
|
||||
% Sliding window (L samples) average
|
||||
function s = PQ_WinAvg (L, x)
|
||||
|
||||
N = length (x);
|
||||
|
||||
s = 0;
|
||||
for (i = L-1:N-1)
|
||||
t = 0;
|
||||
for (m = 0:L-1)
|
||||
t = t + sqrt (x(i-m+1));
|
||||
end
|
||||
s = s + (t / L)^4;
|
||||
end
|
||||
|
||||
if (N >= L)
|
||||
s = sqrt (s / (N - L + 1));
|
||||
end
|
||||
|
||||
%----------
|
||||
% Weighted average
|
||||
function s = PQ_WtAvg (x, W)
|
||||
|
||||
N = length (x);
|
||||
|
||||
s = 0;
|
||||
sW = 0;
|
||||
for (i = 0:N-1)
|
||||
s = s + W(i+1) * x(i+1);
|
||||
sW = sW + W(i+1);
|
||||
end
|
||||
|
||||
if (N > 0)
|
||||
s = s / sW;
|
||||
end
|
||||
|
||||
%----------
|
||||
% Linear average
|
||||
function LinAvg = PQ_LinAvg (x)
|
||||
|
||||
N = length (x);
|
||||
s = 0;
|
||||
for (i = 0:N-1)
|
||||
s = s + x(i+1);
|
||||
end
|
||||
|
||||
LinAvg = s / N;
|
||||
|
||||
%----------
|
||||
% Square root of average of squared values
|
||||
function RMSAvg = PQ_RMSAvg (x)
|
||||
|
||||
N = length (x);
|
||||
s = 0;
|
||||
for (i = 0:N-1)
|
||||
s = s + x(i+1)^2;
|
||||
end
|
||||
|
||||
if (N > 0)
|
||||
RMSAvg = sqrt(s / N);
|
||||
else
|
||||
RMSAvg = 0;
|
||||
end
|
||||
@ -1,65 +0,0 @@
|
||||
function PQframeMOV (i, MOVI)
|
||||
% Copy instantaneous MOV values to a new structure
|
||||
% The output struct MOVC is a global.
|
||||
% For most MOV's, they are just copied to the output structure.
|
||||
% The exception is for the probability of detection, where the
|
||||
% MOV's measure the maximum frequency-by-frequecy between channels.
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $
|
||||
|
||||
global MOVC
|
||||
|
||||
[Nchan,Nc] = size (MOVC.MDiff.Mt1B);
|
||||
|
||||
for (j = 1:Nchan)
|
||||
|
||||
% Modulation differences
|
||||
MOVC.MDiff.Mt1B(j,i+1) = MOVI(j).MDiff.Mt1B;
|
||||
MOVC.MDiff.Mt2B(j,i+1) = MOVI(j).MDiff.Mt2B;
|
||||
MOVC.MDiff.Wt(j,i+1) = MOVI(j).MDiff.Wt;
|
||||
|
||||
% Noise loudness
|
||||
MOVC.NLoud.NL(j,i+1) = MOVI(j).NLoud.NL;
|
||||
|
||||
% Total loudness
|
||||
MOVC.Loud.NRef(j,i+1) = MOVI(j).Loud.NRef;
|
||||
MOVC.Loud.NTest(j,i+1) = MOVI(j).Loud.NTest;
|
||||
|
||||
% Bandwidth
|
||||
MOVC.BW.BWRef(j,i+1) = MOVI(j).BW.BWRef;
|
||||
MOVC.BW.BWTest(j,i+1) = MOVI(j).BW.BWTest;
|
||||
|
||||
% Noise-to-mask ratio
|
||||
MOVC.NMR.NMRavg(j,i+1) = MOVI(j).NMR.NMRavg;
|
||||
MOVC.NMR.NMRmax(j,i+1) = MOVI(j).NMR.NMRmax;
|
||||
|
||||
% Error harmonic structure
|
||||
MOVC.EHS.EHS(j,i+1) = MOVI(j).EHS.EHS;
|
||||
end
|
||||
|
||||
% Probability of detection (collapse frequency bands)
|
||||
[MOVC.PD.Pc(i+1), MOVC.PD.Qc(i+1)] = PQ_ChanPD (MOVI);
|
||||
|
||||
%----------------------------------------
|
||||
function [Pc, Qc] = PQ_ChanPD (MOVI)
|
||||
|
||||
Nc = length (MOVI(1).PD.p);
|
||||
Nchan = length (MOVI);
|
||||
|
||||
Pr = 1;
|
||||
Qc = 0;
|
||||
if (Nchan > 1)
|
||||
for (m = 0:Nc-1)
|
||||
pbin = max (MOVI(1).PD.p(m+1), MOVI(2).PD.p(m+1));
|
||||
qbin = max (MOVI(1).PD.q(m+1), MOVI(2).PD.q(m+1));
|
||||
Pr = Pr * (1 - pbin);
|
||||
Qc = Qc + qbin;
|
||||
end
|
||||
else
|
||||
for (m = 0:Nc-1)
|
||||
Pr = Pr * (1 - MOVI.PD.p(m+1));
|
||||
Qc = Qc + MOVI.PD.q(m+1);
|
||||
end
|
||||
end
|
||||
|
||||
Pc = 1 - Pr;
|
||||
@ -1,29 +0,0 @@
|
||||
function Ndel = PQloudTest (Loud)
|
||||
% Loudness threshold
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $
|
||||
|
||||
[Nchan, Np] = size (Loud.NRef);
|
||||
|
||||
Thr = 0.1;
|
||||
|
||||
% Loudness threshold
|
||||
Ndel = Np;
|
||||
for (j = 0:Nchan-1)
|
||||
Ndel = min (Ndel, PQ_LThresh (Thr, Loud.NRef(j+1,:), Loud.NTest(j+1,:)));
|
||||
end
|
||||
|
||||
%-----------
|
||||
function it = PQ_LThresh (Thr, NRef, NTest)
|
||||
% Loudness check: Look for the first time, the loudness exceeds a threshold
|
||||
% for both the test and reference signals.
|
||||
|
||||
Np = length (NRef);
|
||||
|
||||
it = Np;
|
||||
for (i = 0:Np-1)
|
||||
if (NRef(i+1) > Thr & NTest(i+1) > Thr)
|
||||
it = i;
|
||||
break;
|
||||
end
|
||||
end
|
||||
@ -1,45 +0,0 @@
|
||||
function BW = PQmovBW (X2)
|
||||
% Bandwidth tests
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $
|
||||
|
||||
persistent kx kl FR FT N
|
||||
|
||||
if (isempty (kx))
|
||||
NF = 2048;
|
||||
Fs = 48000;
|
||||
fx = 21586;
|
||||
kx = round (fx / Fs * NF); % 921
|
||||
fl = 8109;
|
||||
kl = round (fl / Fs * NF); % 346
|
||||
FRdB = 10;
|
||||
FR = 10^(FRdB / 10);
|
||||
FTdB = 5;
|
||||
FT = 10^(FTdB / 10);
|
||||
N = NF / 2; % Limit from pseudo-code
|
||||
end
|
||||
|
||||
Xth = X2(2,kx+1);
|
||||
for (k = kx+1:N-1)
|
||||
Xth = max (Xth, X2(2,k+1));
|
||||
end
|
||||
|
||||
% BWRef and BWTest remain negative if the BW of the test signal
|
||||
% does not exceed FR * Xth for kx-1 <= k <= kl+1
|
||||
BW.BWRef = -1;
|
||||
XthR = FR * Xth;
|
||||
for (k = kx-1:-1:kl+1)
|
||||
if (X2(1,k+1) >= XthR)
|
||||
BW.BWRef = k + 1;
|
||||
break;
|
||||
end
|
||||
end
|
||||
|
||||
BW.BWTest = -1;
|
||||
XthT = FT * Xth;
|
||||
for (k = BW.BWRef-1:-1:0)
|
||||
if (X2(2,k+1) >= XthT)
|
||||
BW.BWTest = k + 1;
|
||||
break;
|
||||
end
|
||||
end
|
||||
@ -1,135 +0,0 @@
|
||||
function EHS = PQmovEHS (xR, xT, X2)
|
||||
% Calculate the EHS MOV values
|
||||
|
||||
% P. Kabal $Revision: 1.2 $ $Date: 2004/02/05 04:26:19 $
|
||||
|
||||
persistent NF Nadv NL M Hw
|
||||
|
||||
if (isempty (NL))
|
||||
NF = 2048;
|
||||
Nadv = NF / 2;
|
||||
Fs = 48000;
|
||||
Fmax = 9000;
|
||||
NL = 2^(PQ_log2(NF * Fmax / Fs));
|
||||
M = NL;
|
||||
Hw = (1 / M) * sqrt(8 / 3) * PQHannWin (M);
|
||||
end
|
||||
|
||||
EnThr = 8000;
|
||||
kmax = NL + M - 1;
|
||||
|
||||
EnRef = xR(Nadv+1:NF-1+1) * xR(Nadv+1:NF-1+1)';
|
||||
EnTest = xT(Nadv+1:NF-1+1) * xT(Nadv+1:NF-1+1)';
|
||||
|
||||
% Set the return value to be negative for small energy frames
|
||||
if (EnRef < EnThr & EnTest < EnThr)
|
||||
EHS = -1;
|
||||
return;
|
||||
end
|
||||
|
||||
% Allocate storage
|
||||
D = zeros (1, kmax);
|
||||
|
||||
% Differences of log values
|
||||
for (k = 0:kmax-1)
|
||||
D(k+1) = log (X2(2,k+1) / X2(1,k+1));
|
||||
end
|
||||
|
||||
% Correlation computation
|
||||
C = PQ_Corr (D, NL, M);
|
||||
|
||||
% Normalize the correlations
|
||||
Cn = PQ_NCorr (C, D, NL, M);
|
||||
Cnm = (1 / NL) * sum (Cn(1:NL));
|
||||
|
||||
% Window the correlation
|
||||
Cw = Hw .* (Cn - Cnm);
|
||||
|
||||
% DFT
|
||||
cp = PQRFFT (Cw, NL, 1);
|
||||
|
||||
% Squared magnitude
|
||||
c2 = PQRFFTMSq (cp, NL);
|
||||
|
||||
% Search for a peak after a valley
|
||||
EHS = PQ_FindPeak (c2, NL/2+1);
|
||||
|
||||
%----------------------------------------
|
||||
function log2 = PQ_log2 (a)
|
||||
|
||||
log2 = 0;
|
||||
m = 1;
|
||||
while (m < a)
|
||||
log2 = log2 + 1;
|
||||
m = 2 * m;
|
||||
end
|
||||
log2 = log2 - 1;
|
||||
|
||||
%----------
|
||||
function C = PQ_Corr (D, NL, M)
|
||||
% Correlation calculation
|
||||
|
||||
% Direct computation of the correlation
|
||||
% for (i = 0:NL-1)
|
||||
% s = 0;
|
||||
% for (j = 0:M-1)
|
||||
% s = s + D(j+1) * D(i+j+1);
|
||||
% end
|
||||
% C(i+1) = s;
|
||||
% end
|
||||
|
||||
% Calculate the correlation indirectly
|
||||
NFFT = 2 * NL;
|
||||
D0 = [D(1:M) zeros(1,NFFT-M)];
|
||||
D1 = [D(1:M+NL-1) zeros(1,NFFT-(M+NL-1))];
|
||||
|
||||
% DFTs of the zero-padded sequences
|
||||
d0 = PQRFFT (D0, NFFT, 1);
|
||||
d1 = PQRFFT (D1, NFFT, 1);
|
||||
|
||||
% Multiply (complex) sequences
|
||||
dx(0+1) = d0(0+1) * d1(0+1);
|
||||
for (n = 1:NFFT/2-1)
|
||||
m = NFFT/2 + n;
|
||||
dx(n+1) = d0(n+1) * d1(n+1) + d0(m+1) * d1(m+1);
|
||||
dx(m+1) = d0(n+1) * d1(m+1) - d0(m+1) * d1(n+1);
|
||||
end
|
||||
dx(NFFT/2+1) = d0(NFFT/2+1) * d1(NFFT/2+1);
|
||||
|
||||
% Inverse DFT
|
||||
Cx = PQRFFT (dx, NFFT, -1);
|
||||
C = Cx(1:NL);
|
||||
|
||||
%----------
|
||||
function Cn = PQ_NCorr (C, D, NL, M)
|
||||
% Normalize the correlation
|
||||
|
||||
Cn = zeros (1, NL);
|
||||
|
||||
s0 = C(0+1);
|
||||
sj = s0;
|
||||
Cn(0+1) = 1;
|
||||
for (i = 1:NL-1)
|
||||
sj = sj + (D(i+M-1+1)^2 - D(i-1+1)^2);
|
||||
d = s0 * sj;
|
||||
if (d <= 0)
|
||||
Cn(i+1) = 1;
|
||||
else
|
||||
Cn(i+1) = C(i+1) / sqrt (d);
|
||||
end
|
||||
end
|
||||
|
||||
%----------
|
||||
function EHS = PQ_FindPeak (c2, N)
|
||||
% Search for a peak after a valley
|
||||
|
||||
cprev = c2(0+1);
|
||||
cmax = 0;
|
||||
for (n = 1:N-1)
|
||||
if (c2(n+1) > cprev) % Rising from a valley
|
||||
if (c2(n+1) > cmax)
|
||||
cmax = c2(n+1);
|
||||
end
|
||||
end
|
||||
end
|
||||
EHS = cmax;
|
||||
@ -1,43 +0,0 @@
|
||||
function MDiff = PQmovModDiffB (M, ERavg)
|
||||
% Modulation difference related MOV precursors (Basic version)
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:46 $
|
||||
|
||||
persistent Nc Ete
|
||||
|
||||
if (isempty (Nc))
|
||||
e = 0.3;
|
||||
[Nc, fc] = PQCB ('Basic');
|
||||
Et = PQIntNoise (fc);
|
||||
for (m = 0:Nc-1)
|
||||
Ete(m+1) = Et(m+1)^e;
|
||||
end
|
||||
end
|
||||
|
||||
% Parameters
|
||||
negWt2B = 0.1;
|
||||
offset1B = 1.0;
|
||||
offset2B = 0.01;
|
||||
levWt = 100;
|
||||
|
||||
s1B = 0;
|
||||
s2B = 0;
|
||||
Wt = 0;
|
||||
for (m = 0:Nc-1)
|
||||
if (M(1,m+1) > M(2,m+1))
|
||||
num1B = M(1,m+1) - M(2,m+1);
|
||||
num2B = negWt2B * num1B;
|
||||
else
|
||||
num1B = M(2,m+1) - M(1,m+1);
|
||||
num2B = num1B;
|
||||
end
|
||||
MD1B = num1B / (offset1B + M(1,m+1));
|
||||
MD2B = num2B / (offset2B + M(1,m+1));
|
||||
s1B = s1B + MD1B;
|
||||
s2B = s2B + MD2B;
|
||||
Wt = Wt + ERavg(m+1) / (ERavg(m+1) + levWt * Ete(m+1));
|
||||
end
|
||||
|
||||
MDiff.Mt1B = (100 / Nc) * s1B;
|
||||
MDiff.Mt2B = (100 / Nc) * s2B;
|
||||
MDiff.Wt = Wt;
|
||||
@ -1,33 +0,0 @@
|
||||
function NL = PQmovNLoudB (M, EP)
|
||||
% Noise Loudness
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $
|
||||
|
||||
persistent Nc Et
|
||||
|
||||
if (isempty (Nc))
|
||||
[Nc, fc] = PQCB ('Basic');
|
||||
Et = PQIntNoise (fc);
|
||||
end
|
||||
|
||||
% Parameters
|
||||
alpha = 1.5;
|
||||
TF0 = 0.15;
|
||||
S0 = 0.5;
|
||||
NLmin = 0;
|
||||
e = 0.23;
|
||||
|
||||
s = 0;
|
||||
for (m = 0:Nc-1)
|
||||
sref = TF0 * M(1,m+1) + S0;
|
||||
stest = TF0 * M(2,m+1) + S0;
|
||||
beta = exp (-alpha * (EP(2,m+1) - EP(1,m+1)) / EP(1,m+1));
|
||||
a = max (stest * EP(2,m+1) - sref * EP(1,m+1), 0);
|
||||
b = Et(m+1) + sref * EP(1,m+1) * beta;
|
||||
s = s + (Et(m+1) / stest)^e * ((1 + a / b)^e - 1);
|
||||
end
|
||||
|
||||
NL = (24 / Nc) * s;
|
||||
if (NL < NLmin)
|
||||
NL = 0;
|
||||
end
|
||||
@ -1,36 +0,0 @@
|
||||
function NMR = PQmovNMRB (EbN, Ehs)
|
||||
% Noise-to-mask ratio - Basic version
|
||||
% NMR(1) average NMR
|
||||
% NMR(2) max NMR
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $
|
||||
|
||||
persistent Nc gm
|
||||
|
||||
if (isempty (Nc))
|
||||
[Nc, fc, fl, fu, dz] = PQCB ('Basic');
|
||||
gm = PQ_MaskOffset (dz, Nc);
|
||||
end
|
||||
|
||||
NMR.NMRmax = 0;
|
||||
s = 0;
|
||||
for (m = 0:Nc-1)
|
||||
NMRm = EbN(m+1) / (gm(m+1) * Ehs(m+1));
|
||||
s = s + NMRm;
|
||||
if (NMRm > NMR.NMRmax)
|
||||
NMR.NMRmax = NMRm;
|
||||
end
|
||||
end
|
||||
NMR.NMRavg = s / Nc;
|
||||
|
||||
%----------------------------------------
|
||||
function gm = PQ_MaskOffset (dz, Nc)
|
||||
|
||||
for (m = 0:Nc-1)
|
||||
if (m <= 12 / dz)
|
||||
mdB = 3;
|
||||
else
|
||||
mdB = 0.25 * m * dz;
|
||||
end
|
||||
gm(m+1) = 10^(-mdB / 10);
|
||||
end
|
||||
@ -1,42 +0,0 @@
|
||||
function PD = PQmovPD (Ehs)
|
||||
% Probability of detection
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $
|
||||
|
||||
Nc = length (Ehs);
|
||||
|
||||
% Allocate storage
|
||||
PD.p = zeros (1, Nc);
|
||||
PD.q = zeros (1, Nc);
|
||||
|
||||
persistent c g d1 d2 bP bM
|
||||
|
||||
if (isempty (c))
|
||||
c = [-0.198719 0.0550197 -0.00102438 5.05622e-6 9.01033e-11];
|
||||
d1 = 5.95072;
|
||||
d2 = 6.39468;
|
||||
g = 1.71332;
|
||||
bP = 4;
|
||||
bM = 6;
|
||||
end
|
||||
|
||||
for (m = 0:Nc-1)
|
||||
EdBR = 10 * log10 (Ehs(1,m+1));
|
||||
EdBT = 10 * log10 (Ehs(2,m+1));
|
||||
edB = EdBR - EdBT;
|
||||
if (edB > 0)
|
||||
L = 0.3 * EdBR + 0.7 * EdBT;
|
||||
b = bP;
|
||||
else
|
||||
L = EdBT;
|
||||
b = bM;
|
||||
end
|
||||
if (L > 0)
|
||||
s = d1 * (d2 / L)^g ...
|
||||
+ c(1) + L * (c(2) + L * (c(3) + L * (c(4) + L * c(5))));
|
||||
else
|
||||
s = 1e30;
|
||||
end
|
||||
PD.p(m+1) = 1 - 0.5^((edB / s)^b); % Detection probability
|
||||
PD.q(m+1) = abs (fix(edB)) / s; % Steps above threshold
|
||||
end
|
||||
@ -1,35 +0,0 @@
|
||||
function PQprtMOV (MOV, ODG)
|
||||
% Print MOV values (PEAQ Basic version)
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $
|
||||
|
||||
N = length (MOV);
|
||||
PQ_NMOV_B = 11;
|
||||
PQ_NMOV_A = 5;
|
||||
|
||||
fprintf ('Model Output Variables:\n');
|
||||
if (N == PQ_NMOV_B)
|
||||
fprintf (' BandwidthRefB: %g\n', MOV(1));
|
||||
fprintf (' BandwidthTestB: %g\n', MOV(2));
|
||||
fprintf (' Total NMRB: %g\n', MOV(3));
|
||||
fprintf (' WinModDiff1B: %g\n', MOV(4));
|
||||
fprintf (' ADBB: %g\n', MOV(5));
|
||||
fprintf (' EHSB: %g\n', MOV(6));
|
||||
fprintf (' AvgModDiff1B: %g\n', MOV(7));
|
||||
fprintf (' AvgModDiff2B: %g\n', MOV(8));
|
||||
fprintf (' RmsNoiseLoudB: %g\n', MOV(9));
|
||||
fprintf (' MFPDB: %g\n', MOV(10));
|
||||
fprintf (' RelDistFramesB: %g\n', MOV(11));
|
||||
elseif (N == NMOV_A)
|
||||
fprintf (' RmsModDiffA: %g\n', MOV(1));
|
||||
fprintf (' RmsNoiseLoudAsymA: %g\n', MOV(2));
|
||||
fprintf (' Segmental NMRB: %g\n', MOV(3));
|
||||
fprintf (' EHSB: %g\n', MOV(4));
|
||||
fprintf (' AvgLinDistA: %g\n', MOV(5));
|
||||
else
|
||||
error ('Invalid number of MOVs');
|
||||
end
|
||||
|
||||
fprintf ('Objective Difference Grade: %.3f\n', ODG);
|
||||
|
||||
return;
|
||||
@ -1,40 +0,0 @@
|
||||
function PQprtMOVCi (Nchan, i, MOVC)
|
||||
% Print MOV precursors
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:47 $
|
||||
|
||||
fprintf ('Frame: %d\n', i);
|
||||
|
||||
if (Nchan == 1)
|
||||
fprintf (' Ntot : %g %g\n', ...
|
||||
MOVC.Loud.NRef(1,i+1), MOVC.Loud.NTest(1,i+1));
|
||||
fprintf (' ModDiff: %g %g %g\n', ...
|
||||
MOVC.MDiff.Mt1B(1,i+1), MOVC.MDiff.Mt2B(1,i+1), MOVC.MDiff.Wt(1,i+1));
|
||||
fprintf (' NL : %g\n', MOVC.NLoud.NL(1,i+1));
|
||||
fprintf (' BW : %g %g\n', ...
|
||||
MOVC.BW.BWRef(1,i+1), MOVC.BW.BWTest(1,i+1));
|
||||
fprintf (' NMR : %g %g\n', ...
|
||||
MOVC.NMR.NMRavg(1,i+1), MOVC.NMR.NMRmax(1,i+1));
|
||||
fprintf (' PD : %g %g\n', MOVC.PD.Pc(i+1), MOVC.PD.Qc(i+1));
|
||||
fprintf (' EHS : %g\n', 1000 * MOVC.EHS.EHS(1,i+1));
|
||||
else
|
||||
fprintf (' Ntot : %g %g // %g %g\n', ...
|
||||
MOVC.Loud.NRef(1,i+1), MOVC.Loud.NTest(1,i+1), ...
|
||||
MOVC.Loud.NRef(2,i+1), MOVC.Loud.NTest(2,i+1));
|
||||
fprintf (' ModDiff: %g %g %g // %g %g %g\n', ...
|
||||
MOVC.MDiff.Mt1B(1,i+1), MOVC.MDiff.Mt2B(1,i+1), MOVC.MDiff.Wt(1,i+1), ...
|
||||
MOVC.MDiff.Mt1B(2,i+1), MOVC.MDiff.Mt2B(2,i+1), MOVC.MDiff.Wt(2,i+1));
|
||||
fprintf (' NL : %g // %g\n', ...
|
||||
MOVC.NLoud.NL(1,i+1), ...
|
||||
MOVC.NLoud.NL(2,i+1));
|
||||
fprintf (' BW : %g %g // %g %g\n', ...
|
||||
MOVC.BW.BWRef(1,i+1), MOVC.BW.BWTest(1,i+1), ...
|
||||
MOVC.BW.BWRef(2,i+1), MOVC.BW.BWTest(2,i+1));
|
||||
fprintf (' NMR : %g %g // %g %g\n', ...
|
||||
MOVC.NMR.NMRavg(1,i+1), MOVC.NMR.NMRmax(1,i+1), ...
|
||||
MOVC.NMR.NMRavg(2,i+1), MOVC.NMR.NMRmax(2,i+1));
|
||||
fprintf (' PD : %g %g\n', MOVC.PD.Pc(i+1), MOVC.PD.Qc(i+1));
|
||||
fprintf (' EHS : %g // %g\n', ...
|
||||
1000 * MOVC.EHS.EHS(1,i+1), ...
|
||||
1000 * MOVC.EHS.EHS(2,i+1));
|
||||
end
|
||||
@ -1,10 +0,0 @@
|
||||
function hw = PQHannWin (NF)
|
||||
% Hann window
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $
|
||||
|
||||
hw = zeros (1, NF);
|
||||
|
||||
for (n = 0:NF-1)
|
||||
hw(n+1) = 0.5 * (1 - cos(2 * pi * n / (NF-1)));
|
||||
end
|
||||
@ -1,10 +0,0 @@
|
||||
function EIN = PQIntNoise (f)
|
||||
% Generate the internal noise energy vector
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $
|
||||
|
||||
N = length (f);
|
||||
for (m = 0:N-1)
|
||||
INdB = 1.456 * (f(m+1) / 1000)^(-0.8);
|
||||
EIN(m+1) = 10^(INdB / 10);
|
||||
end
|
||||
@ -1,33 +0,0 @@
|
||||
function X = PQRFFT (x, N, ifn)
|
||||
% Calculate the DFT of a real N-point sequence or the inverse
|
||||
% DFT corresponding to a real N-point sequence.
|
||||
% ifn > 0, forward transform
|
||||
% input x(n) - N real values
|
||||
% output X(k) - The first N/2+1 points are the real
|
||||
% parts of the transform, the next N/2-1 points
|
||||
% are the imaginary parts of the transform. However
|
||||
% the imaginary part for the first point and the
|
||||
% middle point which are known to be zero are not
|
||||
% stored.
|
||||
% ifn < 0, inverse transform
|
||||
% input X(k) - The first N/2+1 points are the real
|
||||
% parts of the transform, the next N/2-1 points
|
||||
% are the imaginary parts of the transform. However
|
||||
% the imaginary part for the first point and the
|
||||
% middle point which are known to be zero are not
|
||||
% stored.
|
||||
% output x(n) - N real values
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $
|
||||
|
||||
if (ifn > 0)
|
||||
X = fft (x, N);
|
||||
XR = real(X(0+1:N/2+1));
|
||||
XI = imag(X(1+1:N/2-1+1));
|
||||
X = [XR XI];
|
||||
else
|
||||
xR = [x(0+1:N/2+1)];
|
||||
xI = [0 x(N/2+1+1:N-1+1) 0];
|
||||
x = complex ([xR xR(N/2-1+1:-1:1+1)], [xI -xI(N/2-1+1:-1:1+1)]);
|
||||
X = real (ifft (x, N));
|
||||
end
|
||||
@ -1,13 +0,0 @@
|
||||
function X2 = PQRFFTMSq (X, N)
|
||||
% Calculate the magnitude squared frequency response from the
|
||||
% DFT values corresponding to a real signal (assumes N is even)
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $
|
||||
|
||||
X2 = zeros (1, N/2+1);
|
||||
|
||||
X2(0+1) = X(0+1)^2;
|
||||
for (k = 1:N/2-1)
|
||||
X2(k+1) = X(k+1)^2 + X(N/2+k+1)^2;
|
||||
end
|
||||
X2(N/2+1) = X(N/2+1)^2;
|
||||
@ -1,13 +0,0 @@
|
||||
function W2 = PQWOME (f)
|
||||
% Generate the weighting for the outer & middle ear filtering
|
||||
% Note: The output is a magnitude-squared vector
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $
|
||||
|
||||
N = length (f);
|
||||
for (k = 0:N-1)
|
||||
fkHz = f(k+1) / 1000;
|
||||
AdB = -2.184 * fkHz^(-0.8) + 6.5 * exp(-0.6 * (fkHz - 3.3)^2) ...
|
||||
- 0.001 * fkHz^(3.6);
|
||||
W2(k+1) = 10^(AdB / 10);
|
||||
end
|
||||
@ -1,106 +0,0 @@
|
||||
function Lim = PQdataBoundary (WAV, Nchan, StartS, Ns)
|
||||
% Search for the data boundaries in a file
|
||||
% StartS - starting sample frame
|
||||
% Ns - Number of sample frames
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $
|
||||
|
||||
PQ_L = 5;
|
||||
Amax = 32768;
|
||||
NBUFF = 2048;
|
||||
PQ_ATHR = 200 * (Amax / 32768);
|
||||
|
||||
% Search from the beginning of the file
|
||||
Lim(1) = -1;
|
||||
is = StartS;
|
||||
EndS = StartS + Ns - 1;
|
||||
while (is <= EndS)
|
||||
Nf = min (EndS - is + 1, NBUFF);
|
||||
x = PQgetData (WAV, is, Nf);
|
||||
for (k = 0:Nchan-1)
|
||||
Lim(1) = max (Lim(1), PQ_DataStart (x(k+1,:), Nf, PQ_L, PQ_ATHR));
|
||||
end
|
||||
if (Lim(1) >= 0)
|
||||
Lim(1) = Lim(1) + is;
|
||||
break
|
||||
end
|
||||
is = is + NBUFF - (PQ_L-1);
|
||||
end
|
||||
|
||||
% Search from the end of the file
|
||||
% This loop is written as if it is going in a forward direction
|
||||
% - When the "forward" position is i, the "backward" position is
|
||||
% EndS - (i - StartS + 1) + 1
|
||||
Lim(2) = -1;
|
||||
is = StartS;
|
||||
while (is <= EndS)
|
||||
Nf = min (EndS - is + 1, NBUFF);
|
||||
ie = is + Nf - 1; % Forward limits [is, ie]
|
||||
js = EndS - (ie - StartS + 1) + 1; % Backward limits [js, js+Nf-1]
|
||||
x = PQgetData (WAV, js, Nf);
|
||||
for (k = 0:Nchan-1)
|
||||
Lim(2) = max (Lim(2), PQ_DataEnd (x(k+1,:), Nf, PQ_L, PQ_ATHR));
|
||||
end
|
||||
if (Lim(2) >= 0)
|
||||
Lim(2) = Lim(2) + js;
|
||||
break
|
||||
end
|
||||
is = is + NBUFF - (PQ_L-1);
|
||||
end
|
||||
|
||||
% Sanity checks
|
||||
if (~ ((Lim(1) >= 0 & Lim(2) >= 0) | (Lim(1) < 0 & Lim(2) < 0)))
|
||||
error ('>>> PQdataBoundary: limits have difference signs');
|
||||
end
|
||||
if (~(Lim(1) <= Lim(2)))
|
||||
error ('>>> PQdataBoundary: Lim(1) > Lim(2)');
|
||||
end
|
||||
|
||||
if (Lim(1) < 0)
|
||||
Lim(1) = 0;
|
||||
Lim(2) = 0;
|
||||
end
|
||||
|
||||
%----------
|
||||
function ib = PQ_DataStart (x, N, L, Thr)
|
||||
|
||||
ib = -1;
|
||||
s = 0;
|
||||
M = min (N, L);
|
||||
for (i = 0:M-1)
|
||||
s = s + abs (x(i+1));
|
||||
end
|
||||
if (s > Thr)
|
||||
ib = 0;
|
||||
return
|
||||
end
|
||||
|
||||
for (i = 1:N-L) % i is the first sample
|
||||
s = s + (abs (x(i+L-1+1)) - abs (x(i-1+1))); % L samples apart
|
||||
if (s > Thr)
|
||||
ib = i;
|
||||
return
|
||||
end
|
||||
end
|
||||
|
||||
%----------
|
||||
function ie = PQ_DataEnd (x, N, L, Thr)
|
||||
|
||||
ie = -1;
|
||||
s = 0;
|
||||
M = min (N, L);
|
||||
for (i = N-M:N-1)
|
||||
s = s + abs (x(i+1));
|
||||
end
|
||||
if (s > Thr)
|
||||
ie = N-1;
|
||||
return
|
||||
end
|
||||
|
||||
for (i = N-2:-1:L-1) % i is the last sample
|
||||
s = s + (abs (x(i-L+1+1)) - abs (x(i+1+1))); % L samples apart
|
||||
if (s > Thr)
|
||||
ie = i;
|
||||
return
|
||||
end
|
||||
end
|
||||
@ -1,59 +0,0 @@
|
||||
function x = PQgetData (WAV, i, N)
|
||||
% Get data from internal buffer or file
|
||||
% i - file position
|
||||
% N - number of samples
|
||||
% x - output data (scaled to the range -32768 to +32767)
|
||||
|
||||
% Only two files can be "active" at a time.
|
||||
% N = 0 resets the buffer
|
||||
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $
|
||||
|
||||
persistent Buff
|
||||
|
||||
iB = WAV.iB + 1;
|
||||
if (N == 0)
|
||||
Buff(iB).N = 20 * 1024; % Fixed size
|
||||
Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N);
|
||||
Buff(iB).i = i;
|
||||
end
|
||||
|
||||
if (N > Buff(iB).N)
|
||||
error ('>>> PQgetData: Request exceeds buffer size');
|
||||
end
|
||||
|
||||
% Check if requested data is not already in the buffer
|
||||
is = i - Buff(iB).i;
|
||||
if (is < 0 | is + N - 1 > Buff(iB).N - 1)
|
||||
Buff(iB).x = PQ_ReadWAV (WAV, i, Buff(iB).N);
|
||||
Buff(iB).i = i;
|
||||
end
|
||||
|
||||
% Copy the data
|
||||
Nchan = WAV.Nchan;
|
||||
is = i - Buff(iB).i;
|
||||
x = Buff(iB).x(1:Nchan,is+1:is+N-1+1);
|
||||
|
||||
%------
|
||||
function x = PQ_ReadWAV (WAV, i, N)
|
||||
% This function considers the data to extended with zeros before and
|
||||
% after the data in the file. If the starting offset i is negative,
|
||||
% zeros are filled in before the data starts at offset 0. If the request
|
||||
% extends beyond the end of data in the file, zeros are appended.
|
||||
|
||||
Amax = 32768;
|
||||
Nchan = WAV.Nchan;
|
||||
|
||||
x = zeros (Nchan, N);
|
||||
|
||||
Nz = 0;
|
||||
if (i < 0)
|
||||
Nz = min (-i, N);
|
||||
i = i + Nz;
|
||||
end
|
||||
|
||||
Ns = min (N - Nz, WAV.Nframe - i);
|
||||
if (i >= 0 & Ns > 0)
|
||||
x(1:Nchan,Nz+1:Nz+Ns-1+1) = Amax * (audioread (WAV.Fname, [i+1 i+Ns-1+1]))';
|
||||
end
|
||||
@ -1,13 +0,0 @@
|
||||
function Fmem = PQinitFMem (Nc, PCinit)
|
||||
% Initialize the filter memories
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:10 $
|
||||
|
||||
Fmem.TDS.Ef(1:2,1:Nc) = 0;
|
||||
Fmem.Adap.P(1:2,1:Nc) = 0;
|
||||
Fmem.Adap.Rn(1:Nc) = 0;
|
||||
Fmem.Adap.Rd(1:Nc) = 0;
|
||||
Fmem.Adap.PC(1:2,1:Nc) = PCinit;
|
||||
Fmem.Env.Ese(1:2,1:Nc) = 0;
|
||||
Fmem.Env.DE(1:2,1:Nc) = 0;
|
||||
Fmem.Env.Eavg(1:2,1:Nc) = 0;
|
||||
@ -1,13 +0,0 @@
|
||||
function [a, b] = PQtConst (t100, tmin, f , Fs)
|
||||
% Calculate the difference equation parameters. The time
|
||||
% constant of the difference equation depends on the center
|
||||
% frequencies.
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $
|
||||
|
||||
N = length (f);
|
||||
for (m = 0:N-1)
|
||||
t = tmin + (100 / f(m+1)) * (t100 - tmin);
|
||||
a(m+1) = exp (-1 / (Fs * t));
|
||||
b(m+1) = (1 - a(m+1));
|
||||
end
|
||||
@ -1,32 +0,0 @@
|
||||
function WAV = PQwavFilePar (File)
|
||||
% Print a WAVE file header, pick up the file parameters
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:34:11 $
|
||||
|
||||
persistent iB
|
||||
|
||||
if (isempty (iB))
|
||||
iB = 0;
|
||||
else
|
||||
iB = mod (iB + 1, 2); % Only two files can be "active" at a time
|
||||
end
|
||||
|
||||
%[size WAV.Fs Nbit] = wavread (File, 'size');
|
||||
[sound, WAV.Fs] = audioread(File);
|
||||
Nbit = 16;
|
||||
WAV.Fname = File;
|
||||
WAV.Nframe = size(sound, 1);
|
||||
WAV.Nchan = size(sound, 2);
|
||||
WAV.iB = iB; % Buffer number
|
||||
|
||||
% Initialize the buffer
|
||||
PQgetData (WAV, 0, 0);
|
||||
|
||||
fprintf (' WAVE file: %s\n', File);
|
||||
if (WAV.Nchan == 1)
|
||||
fprintf (' Number of samples : %d (%.4g s)\n', WAV.Nframe, WAV.Nframe / WAV.Fs);
|
||||
else
|
||||
fprintf (' Number of frames : %d (%.4g s)\n', WAV.Nframe, WAV.Nframe / WAV.Fs);
|
||||
end
|
||||
fprintf (' Sampling frequency: %g\n', WAV.Fs);
|
||||
fprintf (' Number of channels: %d (%d-bit integer)\n', WAV.Nchan, Nbit);
|
||||
@ -1,9 +0,0 @@
|
||||
clear; clc;
|
||||
%addpath('PQevalAudio', 'PQevalAudio/CB','PQevalAudio/Misc','PQevalAudio/MOV', 'PQevalAudio/Patt')
|
||||
|
||||
ref = '2_clean_48.wav' % 16 times oversampling
|
||||
test = '2_noise_48.wav' % 4 times oversampling
|
||||
%test = '../2_output_48.wav' % 4 times oversampling
|
||||
|
||||
|
||||
[odg, movb] = PQevalAudio(ref, test)
|
||||
@ -1,177 +0,0 @@
|
||||
function [ODG, MOVB]= PQevalAudio (Fref, Ftest, StartS, EndS)
|
||||
% Perceptual evaluation of audio quality.
|
||||
|
||||
% - StartS shifts the frames, so that the first frame starts at that sample.
|
||||
% This is a two element array, one element for each input file. If StartS is
|
||||
% a scalar, it applies to both files.
|
||||
% - EndS marks the end of data. The processing stops with the last frame that
|
||||
% contains that sample. This is a two element array, one element for each
|
||||
% input file. If EndS is as scalar, it applies to both files.
|
||||
|
||||
% P. Kabal $Revision: 1.2 $ $Date: 2004/02/05 04:25:24 $
|
||||
|
||||
% Globals (to save on copying in/out of functions)
|
||||
global MOVC PQopt
|
||||
|
||||
% Analysis parameters
|
||||
NF = 2048;
|
||||
Nadv = NF / 2;
|
||||
Version = 'Basic';
|
||||
|
||||
% Options
|
||||
PQopt.ClipMOV = 0;
|
||||
PQopt.PCinit = 0;
|
||||
PQopt.PDfactor = 1;
|
||||
PQopt.Ni = 1;
|
||||
PQopt.DelayOverlap = 1;
|
||||
PQopt.DataBounds = 1;
|
||||
PQopt.EndMin = NF / 2;
|
||||
|
||||
%addpath ('CB', 'MOV', 'Misc', 'Patt');
|
||||
|
||||
if (nargin < 3)
|
||||
StartS = [0, 0];
|
||||
end
|
||||
if (nargin < 4)
|
||||
EndS = [];
|
||||
end
|
||||
|
||||
% Get the number of samples and channels for each file
|
||||
WAV(1) = PQwavFilePar (Fref);
|
||||
WAV(2) = PQwavFilePar (Ftest);
|
||||
|
||||
% Reconcile file differences
|
||||
PQ_CheckWAV (WAV);
|
||||
if (WAV(1).Nframe ~= WAV(2).Nframe)
|
||||
disp ('>>> Number of samples differ: using the minimum');
|
||||
end
|
||||
|
||||
% Data boundaries
|
||||
Nchan = WAV(1).Nchan;
|
||||
[StartS, Fstart, Fend] = PQ_Bounds (WAV, Nchan, StartS, EndS, PQopt);
|
||||
|
||||
% Number of PEAQ frames
|
||||
Np = Fend - Fstart + 1;
|
||||
if (PQopt.Ni < 0)
|
||||
PQopt.Ni = ceil (Np / abs(PQopt.Ni));
|
||||
end
|
||||
|
||||
% Initialize the MOV structure
|
||||
MOVC = PQ_InitMOVC (Nchan, Np);
|
||||
|
||||
% Initialize the filter memory
|
||||
Nc = PQCB (Version);
|
||||
for (j = 0:Nchan-1)
|
||||
Fmem(j+1) = PQinitFMem (Nc, PQopt.PCinit);
|
||||
end
|
||||
|
||||
is = 0;
|
||||
for (i = -Fstart:Np-1)
|
||||
|
||||
% Read a frame of data
|
||||
xR = PQgetData (WAV(1), StartS(1) + is, NF); % Reference file
|
||||
xT = PQgetData (WAV(2), StartS(2) + is, NF); % Test file
|
||||
is = is + Nadv;
|
||||
|
||||
% Process a frame
|
||||
for (j = 0:Nchan-1)
|
||||
[MOVI(j+1), Fmem(j+1)] = PQeval (xR(j+1,:), xT(j+1,:), Fmem(j+1));
|
||||
end
|
||||
|
||||
if (i >= 0)
|
||||
% Move the MOV precursors into a new structure
|
||||
PQframeMOV (i, MOVI); % Output is in global MOVC
|
||||
|
||||
% Print the MOV precursors
|
||||
%if (PQopt.Ni ~= 0 & mod (i, PQopt.Ni) == 0)
|
||||
% PQprtMOVCi (Nchan, i, MOVC);
|
||||
%nd
|
||||
end
|
||||
end
|
||||
|
||||
% Time average of the MOV values
|
||||
if (PQopt.DelayOverlap)
|
||||
Nwup = Fstart;
|
||||
else
|
||||
Nwup = 0;
|
||||
end
|
||||
MOVB = PQavgMOVB (MOVC, Nchan, Nwup);
|
||||
|
||||
% Neural net
|
||||
ODG = PQnNet (MOVB);
|
||||
|
||||
% Summary printout
|
||||
%PQprtMOV (MOVB, ODG);
|
||||
|
||||
%----------
|
||||
function PQ_CheckWAV (WAV)
|
||||
% Check the file parameters
|
||||
|
||||
Fs = 48000;
|
||||
|
||||
if (WAV(1).Nchan ~= WAV(2).Nchan)
|
||||
error ('>>> Number of channels differ');
|
||||
end
|
||||
if (WAV(1).Nchan > 2)
|
||||
error ('>>> Too many input channels');
|
||||
end
|
||||
if (WAV(1).Nframe ~= WAV(2).Nframe)
|
||||
disp ('>>> Number of samples differ');
|
||||
end
|
||||
if (WAV(1).Fs ~= WAV(2).Fs)
|
||||
error ('>>> Sampling frequencies differ');
|
||||
end
|
||||
if (WAV(1).Fs ~= Fs)
|
||||
error ('>>> Invalid Sampling frequency: only 48 kHz supported');
|
||||
end
|
||||
|
||||
%----------
|
||||
function [StartS, Fstart, Fend] = PQ_Bounds (WAV, Nchan, StartS, EndS, PQopt)
|
||||
|
||||
PQ_NF = 2048;
|
||||
PQ_NADV = (PQ_NF / 2);
|
||||
|
||||
if (isempty (StartS))
|
||||
StartS(1) = 0;
|
||||
StartS(2) = 0;
|
||||
elseif (length (StartS) == 1)
|
||||
StartS(2) = StartS(1);
|
||||
end
|
||||
Ns = WAV(1).Nframe;
|
||||
|
||||
% Data boundaries (determined from the reference file)
|
||||
if (PQopt.DataBounds)
|
||||
Lim = PQdataBoundary (WAV(1), Nchan, StartS(1), Ns);
|
||||
%fprintf ('PEAQ Data Boundaries: %ld (%.3f s) - %ld (%.3f s)\n', ...
|
||||
% Lim(1), Lim(1)/WAV(1).Fs, Lim(2), Lim(2)/WAV(1).Fs);
|
||||
else
|
||||
Lim = [Starts(1), StartS(1) + Ns - 1];
|
||||
end
|
||||
|
||||
% Start frame number
|
||||
Fstart = floor ((Lim(1) - StartS(1)) / PQ_NADV);
|
||||
|
||||
% End frame number
|
||||
Fend = floor ((Lim(2) - StartS(1) + 1 - PQopt.EndMin) / PQ_NADV);
|
||||
|
||||
%----------
|
||||
function MOVC = PQ_InitMOVC (Nchan, Np)
|
||||
MOVC.MDiff.Mt1B = zeros (Nchan, Np);
|
||||
MOVC.MDiff.Mt2B = zeros (Nchan, Np);
|
||||
MOVC.MDiff.Wt = zeros (Nchan, Np);
|
||||
|
||||
MOVC.NLoud.NL = zeros (Nchan, Np);
|
||||
|
||||
MOVC.Loud.NRef = zeros (Nchan, Np);
|
||||
MOVC.Loud.NTest = zeros (Nchan, Np);
|
||||
|
||||
MOVC.BW.BWRef = zeros (Nchan, Np);
|
||||
MOVC.BW.BWTest = zeros (Nchan, Np);
|
||||
|
||||
MOVC.NMR.NMRavg = zeros (Nchan, Np);
|
||||
MOVC.NMR.NMRmax = zeros (Nchan, Np);
|
||||
|
||||
MOVC.PD.Pc = zeros (1, Np);
|
||||
MOVC.PD.Qc = zeros (1, Np);
|
||||
|
||||
MOVC.EHS.EHS = zeros (Nchan, Np);
|
||||
@ -1,100 +0,0 @@
|
||||
function ODG = PQnNetB (MOV)
|
||||
% Neural net to get the final ODG
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:27:44 $
|
||||
|
||||
persistent amin amax wx wxb wy wyb bmin bmax I J CLIPMOV
|
||||
global PQopt
|
||||
|
||||
if (isempty (amin))
|
||||
I = length (MOV);
|
||||
if (I == 11)
|
||||
[amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar ('Basic');
|
||||
else
|
||||
[amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar ('Advanced');
|
||||
end
|
||||
[I, J] = size (wx);
|
||||
end
|
||||
|
||||
sigmoid = inline ('1 / (1 + exp(-x))');
|
||||
|
||||
% Scale the MOV's
|
||||
Nclip = 0;
|
||||
MOVx = zeros (1, I);
|
||||
for (i = 0:I-1)
|
||||
MOVx(i+1) = (MOV(i+1) - amin(i+1)) / (amax(i+1) - amin(i+1));
|
||||
if (~ isempty (PQopt) & PQopt.ClipMOV ~= 0)
|
||||
if (MOVx(i+1) < 0)
|
||||
MOVx(i+1) = 0;
|
||||
Nclip = Nclip + 1;
|
||||
elseif (MOVx(i+1) > 1)
|
||||
MOVx(i+1) = 1;
|
||||
Nclip = Nclip + 1;
|
||||
end
|
||||
end
|
||||
end
|
||||
if (Nclip > 0)
|
||||
fprintf ('>>> %d MOVs clipped\n', Nclip);
|
||||
end
|
||||
|
||||
% Neural network
|
||||
DI = wyb;
|
||||
for (j = 0:J-1)
|
||||
arg = wxb(j+1);
|
||||
for (i = 0:I-1)
|
||||
arg = arg + wx(i+1,j+1) * MOVx(i+1);
|
||||
end
|
||||
DI = DI + wy(j+1) * sigmoid (arg);
|
||||
end
|
||||
|
||||
ODG = bmin + (bmax - bmin) * sigmoid (DI);
|
||||
|
||||
function [amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar (Version)
|
||||
|
||||
if (strcmp (Version, 'Basic'))
|
||||
amin = ...
|
||||
[393.916656, 361.965332, -24.045116, 1.110661, -0.206623, ...
|
||||
0.074318, 1.113683, 0.950345, 0.029985, 0.000101, ...
|
||||
0];
|
||||
amax = ...
|
||||
[921, 881.131226, 16.212030, 107.137772, 2.886017, ...
|
||||
13.933351, 63.257874, 1145.018555, 14.819740, 1, ...
|
||||
1];
|
||||
wx = ...
|
||||
[ [ -0.502657, 0.436333, 1.219602 ];
|
||||
[ 4.307481, 3.246017, 1.123743 ];
|
||||
[ 4.984241, -2.211189, -0.192096 ];
|
||||
[ 0.051056, -1.762424, 4.331315 ];
|
||||
[ 2.321580, 1.789971, -0.754560 ];
|
||||
[ -5.303901, -3.452257, -10.814982 ];
|
||||
[ 2.730991, -6.111805, 1.519223 ];
|
||||
[ 0.624950, -1.331523, -5.955151 ];
|
||||
[ 3.102889, 0.871260, -5.922878 ];
|
||||
[ -1.051468, -0.939882, -0.142913 ];
|
||||
[ -1.804679, -0.503610, -0.620456 ] ];
|
||||
wxb = ...
|
||||
[ -2.518254, 0.654841, -2.207228 ];
|
||||
wy = ...
|
||||
[ -3.817048, 4.107138, 4.629582 ];
|
||||
wyb = -0.307594;
|
||||
bmin = -3.98;
|
||||
bmax = 0.22;
|
||||
else
|
||||
amin = ...
|
||||
[ 13.298751, 0.041073, -25.018791, 0.061560, 0.024523 ];
|
||||
amax = ...
|
||||
[ 2166.5, 13.24326, 13.46708, 10.226771, 14.224874 ];
|
||||
wx = ...
|
||||
[ [ 21.211773, -39.913052, -1.382553, -14.545348, -0.320899 ];
|
||||
[ -8.981803, 19.956049, 0.935389, -1.686586, -3.238586 ];
|
||||
[ 1.633830, -2.877505, -7.442935, 5.606502, -1.783120 ];
|
||||
[ 6.103821, 19.587435, -0.240284, 1.088213, -0.511314 ];
|
||||
[ 11.556344, 3.892028, 9.720441, -3.287205, -11.031250 ] ];
|
||||
wxb = ...
|
||||
[ 1.330890, 2.686103, 2.096598, -1.327851, 3.087055 ];
|
||||
wy = ...
|
||||
[ -4.696996, -3.289959, 7.004782, 6.651897, 4.009144 ];
|
||||
wyb = -1.360308;
|
||||
bmin = -3.98;
|
||||
bmax = 0.22;
|
||||
end
|
||||
@ -1,107 +0,0 @@
|
||||
function [EP, Fmem] = PQadapt (Ehs, Fmem, Ver, Mod)
|
||||
% Level and pattern adaptation
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:35:08 $
|
||||
|
||||
persistent a b Nc M1 M2 Version Model
|
||||
|
||||
if (~strcmp (Ver, Version) | ~strcmp (Mod, Model))
|
||||
Version = Ver;
|
||||
Model = Mod;
|
||||
if (strcmp (Model, 'FFT'))
|
||||
[Nc, fc] = PQCB (Version);
|
||||
NF = 2048;
|
||||
Nadv = NF / 2;
|
||||
else
|
||||
[Nc, fc] = PQFB;
|
||||
Nadv = 192;
|
||||
end
|
||||
Version = Ver;
|
||||
Model = Mod;
|
||||
Fs = 48000;
|
||||
Fss = Fs / Nadv;
|
||||
t100 = 0.050;
|
||||
tmin = 0.008;
|
||||
[a b] = PQtConst (t100, tmin, fc, Fss);
|
||||
[M1, M2] = PQ_M1M2 (Version, Model);
|
||||
end
|
||||
|
||||
% Allocate memory
|
||||
EP = zeros (2, Nc);
|
||||
R = zeros (2, Nc);
|
||||
|
||||
% Smooth the excitation patterns
|
||||
% Calculate the correlation terms
|
||||
sn = 0;
|
||||
sd = 0;
|
||||
for (m = 0:Nc-1)
|
||||
Fmem.P(1,m+1) = a(m+1) * Fmem.P(1,m+1) + b(m+1) * Ehs(1,m+1);
|
||||
Fmem.P(2,m+1) = a(m+1) * Fmem.P(2,m+1) + b(m+1) * Ehs(2,m+1);
|
||||
sn = sn + sqrt (Fmem.P(2,m+1) * Fmem.P(1,m+1));
|
||||
sd = sd + Fmem.P(2,m+1);
|
||||
end
|
||||
|
||||
% Level correlation
|
||||
CL = (sn / sd)^2;
|
||||
|
||||
for (m = 0:Nc-1)
|
||||
|
||||
% Scale one of the signals to match levels
|
||||
if (CL > 1)
|
||||
EP(1,m+1) = Ehs(1,m+1) / CL;
|
||||
EP(2,m+1) = Ehs(2,m+1);
|
||||
else
|
||||
EP(1,m+1) = Ehs(1,m+1);
|
||||
EP(2,m+1) = Ehs(2,m+1) * CL;
|
||||
end
|
||||
|
||||
% Calculate a pattern match correction factor
|
||||
Fmem.Rn(m+1) = a(m+1) * Fmem.Rn(m+1) + EP(2,m+1) * EP(1,m+1);
|
||||
Fmem.Rd(m+1) = a(m+1) * Fmem.Rd(m+1) + EP(1,m+1) * EP(1,m+1);
|
||||
if (Fmem.Rd(m+1) <= 0 | Fmem.Rn(m+1) <= 0)
|
||||
error ('>>> PQadap: Rd or Rn is zero');
|
||||
end
|
||||
if (Fmem.Rn(m+1) >= Fmem.Rd(m+1))
|
||||
R(1,m+1) = 1;
|
||||
R(2,m+1) = Fmem.Rd(m+1) / Fmem.Rn(m+1);
|
||||
else
|
||||
R(1,m+1) = Fmem.Rn(m+1) / Fmem.Rd(m+1);
|
||||
R(2,m+1) = 1;
|
||||
end
|
||||
end
|
||||
|
||||
% Average the correction factors over M channels and smooth with time
|
||||
% Create spectrally adapted patterns
|
||||
for (m = 0:Nc-1)
|
||||
iL = max (m - M1, 0);
|
||||
iU = min (m + M2, Nc-1);
|
||||
s1 = 0;
|
||||
s2 = 0;
|
||||
for (i = iL:iU)
|
||||
s1 = s1 + R(1,i+1);
|
||||
s2 = s2 + R(2,i+1);
|
||||
end
|
||||
Fmem.PC(1,m+1) = a(m+1) * Fmem.PC(1,m+1) + b(m+1) * s1 / (iU-iL+1);
|
||||
Fmem.PC(2,m+1) = a(m+1) * Fmem.PC(2,m+1) + b(m+1) * s2 / (iU-iL+1);
|
||||
|
||||
% Final correction factor => spectrally adapted patterns
|
||||
EP(1,m+1) = EP(1,m+1) * Fmem.PC(1,m+1);
|
||||
EP(2,m+1) = EP(2,m+1) * Fmem.PC(2,m+1);
|
||||
end
|
||||
|
||||
%--------------------------------------
|
||||
function [M1, M2] = PQ_M1M2 (Version, Model)
|
||||
% Return band averaging parameters
|
||||
|
||||
if (strcmp (Version, 'Basic'))
|
||||
M1 = 3;
|
||||
M2 = 4;
|
||||
elseif (strcmp (Version, 'Advanced'))
|
||||
if (strcmp (Model, 'FFT'))
|
||||
M1 = 1;
|
||||
M2 = 2;
|
||||
else
|
||||
M1 = 1;
|
||||
M2 = 1;
|
||||
end
|
||||
end
|
||||
@ -1,53 +0,0 @@
|
||||
function Ntot = PQloud (Ehs, Ver, Mod)
|
||||
% Calculate the loudness
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:35:09 $
|
||||
|
||||
e = 0.23;
|
||||
|
||||
persistent Nc s Et Ets Version Model
|
||||
|
||||
if (~strcmp (Ver, Version) | ~strcmp (Mod, Model))
|
||||
Version = Ver;
|
||||
Model = Mod;
|
||||
if (strcmp (Model, 'FFT'))
|
||||
[Nc, fc] = PQCB (Version);
|
||||
c = 1.07664;
|
||||
else
|
||||
[Nc, fc] = PQFB;
|
||||
c = 1.26539;
|
||||
end
|
||||
E0 = 1e4;
|
||||
Et = PQ_enThresh (fc);
|
||||
s = PQ_exIndex (fc);
|
||||
for (m = 0:Nc-1)
|
||||
Ets(m+1) = c * (Et(m+1) / (s(m+1) * E0))^e;
|
||||
end
|
||||
end
|
||||
|
||||
sN = 0;
|
||||
for (m = 0:Nc-1)
|
||||
Nm = Ets(m+1) * ((1 - s(m+1) + s(m+1) * Ehs(m+1) / Et(m+1))^e - 1);
|
||||
sN = sN + max(Nm, 0);
|
||||
end
|
||||
Ntot = (24 / Nc) * sN;
|
||||
|
||||
%====================
|
||||
function s = PQ_exIndex (f)
|
||||
% Excitation index
|
||||
|
||||
N = length (f);
|
||||
for (m = 0:N-1)
|
||||
sdB = -2 - 2.05 * atan(f(m+1) / 4000) - 0.75 * atan((f(m+1) / 1600)^2);
|
||||
s(m+1) = 10^(sdB / 10);
|
||||
end
|
||||
|
||||
%--------------------
|
||||
function Et = PQ_enThresh (f)
|
||||
% Excitation threshold
|
||||
|
||||
N = length (f);
|
||||
for (m = 0:N-1)
|
||||
EtdB = 3.64 * (f(m+1) / 1000)^(-0.8);
|
||||
Et(m+1) = 10^(EtdB / 10);
|
||||
end
|
||||
@ -1,34 +0,0 @@
|
||||
function [M, ERavg, Fmem] = PQmodPatt (Es, Fmem)
|
||||
% Modulation pattern processing
|
||||
|
||||
% P. Kabal $Revision: 1.1 $ $Date: 2003/12/07 13:35:09 $
|
||||
|
||||
persistent Nc a b Fss
|
||||
|
||||
if (isempty (Nc))
|
||||
Fs = 48000;
|
||||
NF = 2048;
|
||||
Fss = Fs / (NF/2);
|
||||
[Nc, fc] = PQCB ('Basic');
|
||||
t100 = 0.050;
|
||||
t0 = 0.008;
|
||||
[a, b] = PQtConst (t100, t0, fc, Fss);
|
||||
end
|
||||
|
||||
% Allocate memory
|
||||
M = zeros (2, Nc);
|
||||
|
||||
e = 0.3;
|
||||
for (i = 1:2)
|
||||
for (m = 0:Nc-1)
|
||||
Ee = Es(i,m+1)^e;
|
||||
Fmem.DE(i,m+1) = a(m+1) * Fmem.DE(i,m+1) ...
|
||||
+ b(m+1) * Fss * abs (Ee - Fmem.Ese(i,m+1));
|
||||
Fmem.Eavg(i,m+1) = a(m+1) * Fmem.Eavg(i,m+1) + b(m+1) * Ee;
|
||||
Fmem.Ese(i,m+1) = Ee;
|
||||
|
||||
M(i,m+1) = Fmem.DE(i,m+1) / (1 + Fmem.Eavg(i,m+1)/0.3);
|
||||
end
|
||||
end
|
||||
|
||||
ERavg = Fmem.Eavg(1,:);
|
||||
@ -1,7 +1,5 @@
|
||||
defaults:
|
||||
- dset: dataset
|
||||
- hydra/job_logging: colorlog
|
||||
- hydra/hydra_logging: colorlog
|
||||
|
||||
path_experiment: "non-specified_name" #there should be a better way to do this
|
||||
tensorboard_logs: "/scratch/work/molinee2/tensorboard_logs/unet_historical" #path with tensorboard
|
||||
@ -87,26 +85,5 @@ hydra:
|
||||
kv_sep: '='
|
||||
item_sep: ','
|
||||
# Remove all paths, as the / in them would mess up things
|
||||
# Remove params that would not impact the training itself
|
||||
# Remove all slurm and submit params.
|
||||
# This is ugly I know...
|
||||
exclude_keys: ['path_experiment',
|
||||
'hydra.job_logging.handles.file.filename']
|
||||
job_logging:
|
||||
handlers:
|
||||
file:
|
||||
class: logging.FileHandler
|
||||
mode: w
|
||||
formatter: colorlog
|
||||
filename: trainer.log
|
||||
console:
|
||||
class: logging.StreamHandler
|
||||
formatter: colorlog
|
||||
stream: ext://sys.stderr
|
||||
|
||||
hydra_logging:
|
||||
handlers:
|
||||
console:
|
||||
class: logging.StreamHandler
|
||||
formatter: colorlog
|
||||
stream: ext://sys.stderr
|
||||
|
||||
151
test.py
151
test.py
@ -1,151 +0,0 @@
|
||||
import os
|
||||
import hydra
|
||||
import logging
|
||||
'''
|
||||
Script used for the objective experiments
|
||||
WARNING: it calls MATLAB to calculate PEAQ and PEMO-Q. The whole process may be very slow
|
||||
'''
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
def run(args):
|
||||
import unet
|
||||
import dataset_loader
|
||||
import tensorflow as tf
|
||||
import pandas as pd
|
||||
|
||||
path_experiment=str(args.path_experiment)
|
||||
|
||||
print(path_experiment)
|
||||
if not os.path.exists(path_experiment):
|
||||
os.makedirs(path_experiment)
|
||||
|
||||
unet_model = unet.build_model_denoise(stereo=stereo,unet_args=args.unet)
|
||||
|
||||
ckpt=os.path.join(path_experiment, 'checkpoint')
|
||||
unet_model.load_weights(ckpt)
|
||||
|
||||
|
||||
path_pianos_test=args.dset.path_piano_test
|
||||
path_strings_test=args.dset.path_strings_test
|
||||
path_orchestra_test=args.dset.path_orchestra_test
|
||||
path_opera_test=args.dset.path_opera_test
|
||||
path_noise=args.dset.path_noise
|
||||
fs=args.fs
|
||||
seg_len_s=20
|
||||
numsamples=1000//seg_len_s
|
||||
|
||||
def do_stft(noisy, clean=None):
|
||||
|
||||
if args.stft.window=="hamming":
|
||||
window_fn = tf.signal.hamming_window
|
||||
elif args.stft.window=="hann":
|
||||
window_fn=tf.signal.hann_window
|
||||
elif args.stft.window=="kaiser_bessel":
|
||||
window_fn=tf.signal.kaiser_bessel_derived_window
|
||||
|
||||
win_size=args.stft.win_size
|
||||
hop_size=args.stft.hop_size
|
||||
|
||||
|
||||
stft_signal_noisy=tf.signal.stft(noisy,frame_length=win_size, window_fn=window_fn, frame_step=hop_size)
|
||||
stft_noisy_stacked=tf.stack( values=[tf.math.real(stft_signal_noisy), tf.math.imag(stft_signal_noisy)], axis=-1)
|
||||
|
||||
if clean!=None:
|
||||
|
||||
stft_signal_clean=tf.signal.stft(clean,frame_length=win_size, window_fn=window_fn, frame_step=hop_size)
|
||||
stft_clean_stacked=tf.stack( values=[tf.math.real(stft_signal_clean), tf.math.imag(stft_signal_clean)], axis=-1)
|
||||
|
||||
|
||||
return stft_noisy_stacked, stft_clean_stacked
|
||||
else:
|
||||
|
||||
return stft_noisy_stacked
|
||||
|
||||
from tester import Tester
|
||||
|
||||
testPath=os.path.join(path_experiment,"final_test")
|
||||
if not os.path.exists(testPath):
|
||||
os.makedirs(testPath)
|
||||
|
||||
tester=Tester(unet_model, testPath, args)
|
||||
|
||||
PEAQ_dir="/scratch/work/molinee2/unet_dir/unet_historical_music/PQevalAudio"
|
||||
PEMOQ_dir="/scratch/work/molinee2/unet_dir/unet_historical_music/PEMOQ"
|
||||
|
||||
dataset_test_pianos=dataset_loader.load_data_formal( path_pianos_test, path_noise, noise_amount="mid_snr",num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_pianos=dataset_test_pianos.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_pianos,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("pianos_midsnr")
|
||||
|
||||
dataset_test_strings=dataset_loader.load_data_formal( path_strings_test, path_noise,noise_amount="mid_snr",num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_strings=dataset_test_strings.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_strings,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("strings_midsnr")
|
||||
|
||||
dataset_test_orchestra=dataset_loader.load_data_formal( path_orchestra_test, path_noise, noise_amount="mid_snr", num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_orchestra=dataset_test_orchestra.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_orchestra,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("orchestra_midsnr")
|
||||
|
||||
dataset_test_opera=dataset_loader.load_data_formal( path_opera_test, path_noise, noise_amount="mid_snr",num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_opera=dataset_test_opera.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_opera,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("opera_midsnr")
|
||||
|
||||
dataset_test_strings=dataset_loader.load_data_formal( path_strings_test, path_noise,noise_amount="low_snr",num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_strings=dataset_test_strings.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_strings,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("strings_lowsnr")
|
||||
|
||||
dataset_test_orchestra=dataset_loader.load_data_formal( path_orchestra_test, path_noise,noise_amount="low_snr", num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_orchestra=dataset_test_orchestra.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_orchestra,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("orchestra_lowsnr")
|
||||
|
||||
dataset_test_opera=dataset_loader.load_data_formal( path_opera_test, path_noise, noise_amount="low_snr",num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_opera=dataset_test_opera.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_opera,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("opera_lowsnr")
|
||||
|
||||
|
||||
dataset_test_pianos=dataset_loader.load_data_formal( path_pianos_test, path_noise, noise_amount="low_snr",num_samples=numsamples, fs=fs, seg_len_s=seg_len_s, stereo=stereo)
|
||||
dataset_test_pianos=dataset_test_pianos.map(do_stft, num_parallel_calls=args.num_workers, deterministic=None)
|
||||
tester.init_inference(dataset_test_pianos,numsamples,fs,args.stft, PEAQ_dir, PEMOQ_dir=PEMOQ_dir)
|
||||
metrics=tester.inference("pianos_lowsnr")
|
||||
|
||||
|
||||
names=["strings_midsnr","strings_lowsnr","opera_midsnr","opera_lowsnr","pianos_midsnr","pianos_lowsnr","orchestra_midsnr","orchestra_lowsnr"]
|
||||
for n in names:
|
||||
a=pd.read_csv(os.path.join(testPath,n,"metrics.csv"))
|
||||
meanPEAQ=a["PEAQ(ODG)_diff"].sum()/50
|
||||
meanPEMOQ=a["PEMOQ(ODG)_diff"].sum()/50
|
||||
meanSDR=a["SDR_diff"].sum()/50
|
||||
print(n,": PEAQ ",str(meanPEAQ), "PEMOQ ", str(meanPEMOQ), "SDR ", str(meanSDR))
|
||||
|
||||
def _main(args):
|
||||
global __file__
|
||||
|
||||
__file__ = hydra.utils.to_absolute_path(__file__)
|
||||
|
||||
run(args)
|
||||
|
||||
|
||||
@hydra.main(config_path="conf/conf.yaml")
|
||||
def main(args):
|
||||
try:
|
||||
_main(args)
|
||||
except Exception:
|
||||
logger.exception("Some error happened")
|
||||
# Hydra intercepts exit code, fixed in beta but I could not get the beta to work
|
||||
os._exit(1)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
391
tester.py
391
tester.py
@ -1,391 +0,0 @@
|
||||
|
||||
import os
|
||||
import numpy as np
|
||||
import cv2
|
||||
import librosa
|
||||
import imageio
|
||||
import tensorflow as tf
|
||||
import soundfile as sf
|
||||
import subprocess
|
||||
from tqdm import tqdm
|
||||
from vggish.vgg_distance import process_wav
|
||||
import pandas as pd
|
||||
from scipy.io import loadmat
|
||||
|
||||
class Tester():
|
||||
def __init__(self, model, path_experiment, args):
|
||||
if model !=None:
|
||||
self.model=model
|
||||
print(self.model.summary())
|
||||
self.args=args
|
||||
self.path_experiment=path_experiment
|
||||
|
||||
def init_inference(self, dataset_test=None,num_test_segments=0 , fs=44100, stft_args=None, PEAQ_dir=None, alg_dir=None, PEMOQ_dir=None):
|
||||
|
||||
self.num_test_segments=num_test_segments
|
||||
self.dataset_test=dataset_test
|
||||
|
||||
if self.dataset_test!=None:
|
||||
self.dataset_test=self.dataset_test.take(self.num_test_segments)
|
||||
|
||||
self.fs=fs
|
||||
self.stft_args=stft_args
|
||||
self.win_size=stft_args.win_size
|
||||
self.hop_size=stft_args.hop_size
|
||||
self.window=stft_args.window
|
||||
self.PEAQ_dir=PEAQ_dir
|
||||
self.PEMOQ_dir=PEMOQ_dir
|
||||
self.alg_dir=alg_dir
|
||||
|
||||
|
||||
|
||||
|
||||
def generate_inverse_window(self, stft_args):
|
||||
if stft_args.window=="hamming":
|
||||
return tf.signal.inverse_stft_window_fn(stft_args.hop_size, forward_window_fn=tf.signal.hamming_window)
|
||||
elif stft_args.window=="hann":
|
||||
return tf.signal.inverse_stft_window_fn(stft_args.hop_size, forward_window_fn=tf.signal.hann_window)
|
||||
elif stft_args.window=="kaiser_bessel":
|
||||
return tf.signal.inverse_stft_window_fn(stft_args.hop_size, forward_window_fn=tf.signal.kaiser_bessel_derived_window)
|
||||
def do_istft(self,data):
|
||||
|
||||
window_fn = self.generate_inverse_window(self.stft_args)
|
||||
win_size=self.win_size
|
||||
hop_size=self.hop_size
|
||||
pred_cpx=data[...,0] + 1j * data[...,1]
|
||||
pred_time=tf.signal.inverse_stft(pred_cpx, win_size, hop_size, window_fn=window_fn)
|
||||
return pred_time
|
||||
|
||||
def generate_images(self,cpx,name):
|
||||
spectro=np.clip((np.flipud(np.transpose(10*np.log10(np.sqrt(np.power(cpx[...,0],2)+np.power(cpx[...,1],2)))))+30)/50,0,1)
|
||||
spectrorgb=np.zeros(shape=(spectro.shape[0],spectro.shape[1],3))
|
||||
spectrorgb[...,0]=np.clip((np.flipud(np.transpose(10*np.log10(np.abs(cpx[...,0])+0.001)))+30)/50,0,1)
|
||||
spectrorgb[...,1]=np.clip((np.flipud(np.transpose(10*np.log10(np.abs(cpx[...,1])+0.001)))+30)/50,0,1)
|
||||
cmap=cv2.COLORMAP_JET
|
||||
spectro = np.array((1-spectro)* 255, dtype = np.uint8)
|
||||
spectro = cv2.applyColorMap(spectro, cmap)
|
||||
imageio.imwrite(os.path.join(self.test_results_filepath, name+".png"),spectro)
|
||||
spectrorgb = np.array(spectrorgb* 255, dtype = np.uint8)
|
||||
imageio.imwrite(os.path.join(self.test_results_filepath, name+"_ir.png"),spectrorgb)
|
||||
|
||||
def generate_image_diff(self,clean , pred,name):
|
||||
difference=np.sqrt((clean[...,0]-pred[...,0])**2+(clean[...,1]-pred[...,1])**2)
|
||||
dif=np.clip(np.flipud(np.transpose(difference)),0,1)
|
||||
cmap=cv2.COLORMAP_JET
|
||||
dif = np.array((1-dif)* 255, dtype = np.uint8)
|
||||
dif = cv2.applyColorMap(dif, cmap)
|
||||
imageio.imwrite(os.path.join(self.test_results_filepath, name+"_diff.png"),dif)
|
||||
|
||||
def inference_inner_classical(self, folder_name, method):
|
||||
nums=[]
|
||||
|
||||
PEAQ_odg_noisy=[]
|
||||
PEAQ_odg_output=[]
|
||||
PEAQ_odg_diff=[]
|
||||
|
||||
PEMOQ_odg_noisy=[]
|
||||
PEMOQ_odg_output=[]
|
||||
PEMOQ_odg_diff=[]
|
||||
|
||||
SDR_noisy=[]
|
||||
SDR_output=[]
|
||||
SDR_diff=[]
|
||||
|
||||
VGGish_noisy=[]
|
||||
VGGish_output=[]
|
||||
VGGish_diff=[]
|
||||
|
||||
self.test_results_filepath = os.path.join(self.path_experiment,folder_name)
|
||||
if not os.path.exists(self.test_results_filepath):
|
||||
os.makedirs(self.test_results_filepath)
|
||||
num=0
|
||||
for element in tqdm(self.dataset_test.take(self.num_test_segments)):
|
||||
test_element=tf.data.Dataset.from_tensors(element)
|
||||
noisy_time=element[0].numpy()
|
||||
#noisy_time=self.do_istft(noisy)
|
||||
name_noisy=str(num)+'_noisy'
|
||||
clean_time=element[1].numpy()
|
||||
#clean_time=self.do_istft(clean)
|
||||
name_clean=str(num)+'_clean'
|
||||
print("inferencing")
|
||||
|
||||
|
||||
nums.append(num)
|
||||
|
||||
print("generating wavs")
|
||||
#noisy_time=noisy_time.numpy().astype(np.float32)
|
||||
noisy_time=noisy_time.astype(np.float32)
|
||||
wav_noisy_name_pre=os.path.join(self.test_results_filepath, name_noisy+"pre.wav")
|
||||
sf.write(wav_noisy_name_pre, noisy_time, 44100)
|
||||
|
||||
#pred = self.model.predict(test_element.batch(1))
|
||||
name_pred=str(num)+'_output'
|
||||
wav_output_name_proc=os.path.join(self.test_results_filepath, name_pred+"proc.wav")
|
||||
self.process_in_matlab(wav_noisy_name_pre, wav_output_name_proc, method)
|
||||
|
||||
noisy_time=noisy_time[44100::] #remove pre noise
|
||||
|
||||
#clean_time=clean_time.numpy().astype(np.float32)
|
||||
clean_time=clean_time.astype(np.float32)
|
||||
clean_time=clean_time[44100::] #remove pre noise
|
||||
|
||||
#change that !!!!
|
||||
#pred_time=self.do_istft(pred[0])
|
||||
#pred_time=pred_time.numpy().astype(np.float32)
|
||||
#pred_time=librosa.resample(np.transpose(pred_time),self.fs, 48000)
|
||||
#sf.write(wav_output_name, pred_time, 48000)
|
||||
#LOAD THE AUDIO!!!
|
||||
pred_time, sr=sf.read(wav_output_name_proc)
|
||||
assert sr==44100
|
||||
pred_time=pred_time[44100::] #remove prenoise
|
||||
|
||||
#I am computing here the SDR at 48k, whle I was doing it before at 44.1k. I hope this won't cause any problem in the results. Consider resampling???
|
||||
SDR_t_noisy=10*np.log10(np.mean(np.square(clean_time))/np.mean(np.square(noisy_time-clean_time)))
|
||||
SDR_noisy.append(SDR_t_noisy)
|
||||
SDR_t_output=10*np.log10(np.mean(np.square(clean_time))/np.mean(np.square(pred_time-clean_time)))
|
||||
SDR_output.append(SDR_t_output)
|
||||
SDR_diff.append(SDR_t_output-SDR_t_noisy)
|
||||
|
||||
noisy_time=librosa.resample(np.transpose(noisy_time),self.fs, 48000) #P.Kabal PEAQ code is hardcoded at Fs=48000, so we have to resample
|
||||
wav_noisy_name=os.path.join(self.test_results_filepath, name_noisy+".wav")
|
||||
sf.write(wav_noisy_name, noisy_time, 48000) #overwrite without prenoise
|
||||
|
||||
clean_time=librosa.resample(np.transpose(clean_time),self.fs, 48000) #without prenoise please!!!
|
||||
wav_clean_name=os.path.join(self.test_results_filepath, name_clean+".wav")
|
||||
sf.write(wav_clean_name, clean_time, 48000)
|
||||
|
||||
pred_time=librosa.resample(np.transpose(pred_time),self.fs, 48000) #without prenoise please!!!
|
||||
wav_output_name=os.path.join(self.test_results_filepath, name_pred+".wav")
|
||||
sf.write(wav_output_name, pred_time, 48000)
|
||||
|
||||
#save pred at 48k
|
||||
#print("calculating PEMOQ")
|
||||
#odg_noisy,odg_output =self.calculate_PEMOQ(wav_clean_name,wav_noisy_name,wav_output_name)
|
||||
#PEMOQ_odg_noisy.append(odg_noisy)
|
||||
#PEMOQ_odg_output.append(odg_output)
|
||||
#PEMOQ_odg_diff.append(odg_output-odg_noisy)
|
||||
|
||||
#print("calculating PEAQ")
|
||||
#odg_noisy,odg_output =self.calculate_PEAQ(wav_clean_name,wav_noisy_name,wav_output_name)
|
||||
#PEAQ_odg_noisy.append(odg_noisy)
|
||||
#PEAQ_odg_output.append(odg_output)
|
||||
#PEAQ_odg_diff.append(odg_output-odg_noisy)
|
||||
|
||||
print("calculating VGGish")
|
||||
VGGish_clean_embeddings=process_wav(wav_clean_name)
|
||||
VGGish_noisy_embeddings=process_wav(wav_noisy_name)
|
||||
VGGish_output_embeddings=process_wav(wav_output_name)
|
||||
dist_noisy = np.linalg.norm(VGGish_noisy_embeddings-VGGish_clean_embeddings)
|
||||
dist_output = np.linalg.norm(VGGish_output_embeddings-VGGish_clean_embeddings)
|
||||
VGGish_noisy.append(dist_noisy)
|
||||
VGGish_output.append(dist_output)
|
||||
VGGish_diff.append(-(dist_output-dist_noisy))
|
||||
os.remove(wav_clean_name)
|
||||
os.remove(wav_noisy_name)
|
||||
os.remove(wav_noisy_name_pre)
|
||||
os.remove(wav_output_name)
|
||||
os.remove(wav_output_name_proc)
|
||||
|
||||
num=num+1
|
||||
|
||||
frame = { 'num':nums,'PEAQ(ODG)_noisy': PEAQ_odg_noisy, 'PEAQ(ODG)_output': PEAQ_odg_output, 'PEAQ(ODG)_diff': PEAQ_odg_diff, 'PEMOQ(ODG)_noisy': PEMOQ_odg_noisy, 'PEMOQ(ODG)_output': PEMOQ_odg_output, 'PEMOQ(ODG)_diff': PEMOQ_odg_diff,'SDR_noisy': SDR_noisy, 'SDR_output': SDR_output, 'SDR_diff': SDR_diff, 'VGGish_noisy': VGGish_noisy, 'VGGish_output': VGGish_output,'VGGish_diff': VGGish_diff }
|
||||
|
||||
metrics=pd.DataFrame(frame)
|
||||
metrics.to_csv(os.path.join(self.test_results_filepath,"metrics.csv"),index=False)
|
||||
metrics=metrics.set_index('num')
|
||||
|
||||
return metrics
|
||||
def inference_inner(self, folder_name):
|
||||
nums=[]
|
||||
|
||||
PEAQ_odg_noisy=[]
|
||||
PEAQ_odg_output=[]
|
||||
PEAQ_odg_diff=[]
|
||||
|
||||
PEMOQ_odg_noisy=[]
|
||||
PEMOQ_odg_output=[]
|
||||
PEMOQ_odg_diff=[]
|
||||
|
||||
SDR_noisy=[]
|
||||
SDR_output=[]
|
||||
SDR_diff=[]
|
||||
|
||||
VGGish_noisy=[]
|
||||
VGGish_output=[]
|
||||
VGGish_diff=[]
|
||||
|
||||
self.test_results_filepath = os.path.join(self.path_experiment,folder_name)
|
||||
if not os.path.exists(self.test_results_filepath):
|
||||
os.makedirs(self.test_results_filepath)
|
||||
num=0
|
||||
for element in tqdm(self.dataset_test.take(self.num_test_segments)):
|
||||
test_element=tf.data.Dataset.from_tensors(element)
|
||||
noisy=element[0].numpy()
|
||||
noisy_time=self.do_istft(noisy)
|
||||
name_noisy=str(num)+'_noisy'
|
||||
clean=element[1].numpy()
|
||||
clean_time=self.do_istft(clean)
|
||||
name_clean=str(num)+'_clean'
|
||||
print("inferencing")
|
||||
pred = self.model.predict(test_element.batch(1))
|
||||
if self.args.unet.num_stages==2:
|
||||
pred=pred[0]
|
||||
pred_time=self.do_istft(pred[0])
|
||||
name_pred=str(num)+'_output'
|
||||
|
||||
nums.append(num)
|
||||
pred_time=pred_time.numpy().astype(np.float32)
|
||||
clean_time=clean_time.numpy().astype(np.float32)
|
||||
SDR_t_noisy=10*np.log10(np.mean(np.square(clean_time))/np.mean(np.square(noisy_time-clean_time)))
|
||||
SDR_t_output=10*np.log10(np.mean(np.square(clean_time))/np.mean(np.square(pred_time-clean_time)))
|
||||
SDR_noisy.append(SDR_t_noisy)
|
||||
SDR_output.append(SDR_t_output)
|
||||
SDR_diff.append(SDR_t_output-SDR_t_noisy)
|
||||
|
||||
print("generating wavs")
|
||||
noisy_time=librosa.resample(np.transpose(noisy_time),self.fs, 48000) #P.Kabal PEAQ code is hardcoded at Fs=48000, so we have to resample
|
||||
clean_time=librosa.resample(np.transpose(clean_time),self.fs, 48000)
|
||||
pred_time=librosa.resample(np.transpose(pred_time),self.fs, 48000)
|
||||
|
||||
wav_noisy_name=os.path.join(self.test_results_filepath, name_noisy+".wav")
|
||||
sf.write(wav_noisy_name, noisy_time, 48000)
|
||||
wav_clean_name=os.path.join(self.test_results_filepath, name_clean+".wav")
|
||||
sf.write(wav_clean_name, clean_time, 48000)
|
||||
wav_output_name=os.path.join(self.test_results_filepath, name_pred+".wav")
|
||||
sf.write(wav_output_name, pred_time, 48000)
|
||||
|
||||
print("calculating PEMOQ")
|
||||
odg_noisy,odg_output =self.calculate_PEMOQ(wav_clean_name,wav_noisy_name,wav_output_name)
|
||||
PEMOQ_odg_noisy.append(odg_noisy)
|
||||
PEMOQ_odg_output.append(odg_output)
|
||||
PEMOQ_odg_diff.append(odg_output-odg_noisy)
|
||||
print("calculating PEAQ")
|
||||
odg_noisy,odg_output =self.calculate_PEAQ(wav_clean_name,wav_noisy_name,wav_output_name)
|
||||
PEAQ_odg_noisy.append(odg_noisy)
|
||||
PEAQ_odg_output.append(odg_output)
|
||||
PEAQ_odg_diff.append(odg_output-odg_noisy)
|
||||
|
||||
print("calculating VGGish")
|
||||
VGGish_clean_embeddings=process_wav(wav_clean_name)
|
||||
VGGish_noisy_embeddings=process_wav(wav_noisy_name)
|
||||
VGGish_output_embeddings=process_wav(wav_output_name)
|
||||
dist_noisy = np.linalg.norm(VGGish_noisy_embeddings-VGGish_clean_embeddings)
|
||||
dist_output = np.linalg.norm(VGGish_output_embeddings-VGGish_clean_embeddings)
|
||||
VGGish_noisy.append(dist_noisy)
|
||||
VGGish_output.append(dist_output)
|
||||
VGGish_diff.append(-(dist_output-dist_noisy))
|
||||
os.remove(wav_clean_name)
|
||||
os.remove(wav_noisy_name)
|
||||
os.remove(wav_output_name)
|
||||
|
||||
num=num+1
|
||||
|
||||
frame = { 'num':nums,'PEAQ(ODG)_noisy': PEAQ_odg_noisy, 'PEAQ(ODG)_output': PEAQ_odg_output, 'PEAQ(ODG)_diff': PEAQ_odg_diff, 'PEMOQ(ODG)_noisy': PEMOQ_odg_noisy, 'PEMOQ(ODG)_output': PEMOQ_odg_output, 'PEMOQ(ODG)_diff': PEMOQ_odg_diff,'SDR_noisy': SDR_noisy, 'SDR_output': SDR_output, 'SDR_diff': SDR_diff, 'VGGish_noisy': VGGish_noisy, 'VGGish_output': VGGish_output,'VGGish_diff': VGGish_diff }
|
||||
|
||||
metrics=pd.DataFrame(frame)
|
||||
metrics.to_csv(os.path.join(self.test_results_filepath,"metrics.csv"),index=False)
|
||||
metrics=metrics.set_index('num')
|
||||
|
||||
return metrics
|
||||
|
||||
|
||||
def inference_real(self, folder_name):
|
||||
self.test_results_filepath = os.path.join(self.path_experiment,folder_name)
|
||||
if not os.path.exists(self.test_results_filepath):
|
||||
os.makedirs(self.test_results_filepath)
|
||||
num=0
|
||||
for element in tqdm(self.dataset_real.take(self.num_real_test_segments)):
|
||||
test_element=tf.data.Dataset.from_tensors(element)
|
||||
noisy=element.numpy()
|
||||
noisy_time=self.do_istft(noisy)
|
||||
name_noisy="recording_"+str(num)+'_noisy.wav'
|
||||
pred = self.model.predict(test_element.batch(1))
|
||||
if self.args.unet.num_stages==2:
|
||||
pred=pred[0]
|
||||
pred_time=self.do_istft(pred[0])
|
||||
name_pred="recording_"+str(num)+'_output.wav'
|
||||
sf.write(os.path.join(self.test_results_filepath, name_noisy), noisy_time, self.fs)
|
||||
sf.write(os.path.join(self.test_results_filepath, name_pred), pred_time, self.fs)
|
||||
self.generate_images(noisy,name_noisy)
|
||||
self.generate_images(pred[0],name_pred)
|
||||
num=num+1
|
||||
|
||||
|
||||
def process_in_matlab(self,wav_noisy_name,wav_output_name,mode): #Opening and closing matlab to calculate PEAQ, rudimentary way to do it but easier. Make sure to have matlab installed
|
||||
addpath=self.alg_dir
|
||||
#odgmatfile_noisy=os.path.join(self.test_results_filepath, "odg_noisy.mat")
|
||||
#odgmatfile_pred=os.path.join(self.test_results_filepath, "odg_pred.mat")
|
||||
#bashCommand = "matlab -nodesktop -r 'addpath(\"PQevalAudio\", \"PQevalAudio/CB\",\"PQevalAudio/Misc\",\"PQevalAudio/MOV\", \"PQevalAudio/Patt\"), [odg, MOV]=PQevalAudio(\"0_clean_48.wav\",\"0_noise_48.wav\"), save(\"odg_noisy.mat\",\"odg\"), save(\"mov.mat\",\"MOV\") , exit'"
|
||||
bashCommand = "matlab -nodesktop -r 'addpath(genpath(\""+addpath+"\")), declick_and_denoise(\""+wav_noisy_name+"\",\""+wav_output_name+"\",\""+mode+"\") , exit'"
|
||||
print(bashCommand)
|
||||
p1 = subprocess.Popen(bashCommand, stdout=subprocess.PIPE, shell=True)
|
||||
(output, err) = p1.communicate()
|
||||
|
||||
print(output)
|
||||
|
||||
p1.wait()
|
||||
|
||||
def calculate_PEMOQ(self,wav_clean_name,wav_noisy_name,wav_output_name): #Opening and closing matlab to calculate PEAQ, rudimentary way to do it but easier. Make sure to have matlab installed
|
||||
addpath=self.PEMOQ_dir
|
||||
odgmatfile_noisy=os.path.join(self.test_results_filepath, "odg_pemo_noisy.mat")
|
||||
odgmatfile_pred=os.path.join(self.test_results_filepath, "odg_pemo_pred.mat")
|
||||
#bashCommand = "matlab -nodesktop -r 'addpath(\"PQevalAudio\", \"PQevalAudio/CB\",\"PQevalAudio/Misc\",\"PQevalAudio/MOV\", \"PQevalAudio/Patt\"), [odg, MOV]=PQevalAudio(\"0_clean_48.wav\",\"0_noise_48.wav\"), save(\"odg_noisy.mat\",\"odg\"), save(\"mov.mat\",\"MOV\") , exit'"
|
||||
bashCommand = "matlab -nodesktop -r 'addpath(genpath(\""+addpath+"\")), [ ODG]=PEMOQ(\""+wav_clean_name+"\",\""+wav_noisy_name+"\"), save(\""+odgmatfile_noisy+"\",\"ODG\"), exit'"
|
||||
print(bashCommand)
|
||||
|
||||
p1 = subprocess.Popen(bashCommand, stdout=subprocess.PIPE, shell=True)
|
||||
(output, err) = p1.communicate()
|
||||
|
||||
print(output)
|
||||
|
||||
bashCommand = "matlab -nodesktop -r 'addpath(genpath(\""+addpath+"\")), [ ODG]=PEMOQ(\""+wav_clean_name+"\",\""+wav_output_name+"\"), save(\""+odgmatfile_pred+"\",\"ODG\"), exit'"
|
||||
|
||||
p2 = subprocess.Popen(bashCommand, stdout=subprocess.PIPE, shell=True)
|
||||
(output, err) = p2.communicate()
|
||||
|
||||
print(output)
|
||||
p1.wait()
|
||||
p2.wait()
|
||||
#I save the odg results in a .mat file, which I load here. Not the most optimal method, sorry :/
|
||||
annots_noise = loadmat(odgmatfile_noisy)
|
||||
annots_pred = loadmat(odgmatfile_pred)
|
||||
#Consider loading also the movs!!
|
||||
return annots_noise["ODG"][0][0], annots_pred["ODG"][0][0]
|
||||
|
||||
def calculate_PEAQ(self,wav_clean_name,wav_noisy_name,wav_output_name): #Opening and closing matlab to calculate PEAQ, rudimentary way to do it but easier. Make sure to have matlab installed
|
||||
addpath=self.PEAQ_dir
|
||||
odgmatfile_noisy=os.path.join(self.test_results_filepath, "odg_noisy.mat")
|
||||
odgmatfile_pred=os.path.join(self.test_results_filepath, "odg_pred.mat")
|
||||
#bashCommand = "matlab -nodesktop -r 'addpath(\"PQevalAudio\", \"PQevalAudio/CB\",\"PQevalAudio/Misc\",\"PQevalAudio/MOV\", \"PQevalAudio/Patt\"), [odg, MOV]=PQevalAudio(\"0_clean_48.wav\",\"0_noise_48.wav\"), save(\"odg_noisy.mat\",\"odg\"), save(\"mov.mat\",\"MOV\") , exit'"
|
||||
bashCommand = "matlab -nodesktop -r 'addpath(genpath(\""+addpath+"\")), [odg, MOV]=PQevalAudio(\""+wav_clean_name+"\",\""+wav_noisy_name+"\"), save(\""+odgmatfile_noisy+"\",\"odg\"), save(\"mov.mat\",\"MOV\") , exit'"
|
||||
p1 = subprocess.Popen(bashCommand, stdout=subprocess.PIPE, shell=True)
|
||||
(output, err) = p1.communicate()
|
||||
|
||||
print(output)
|
||||
|
||||
bashCommand = "matlab -nodesktop -r 'addpath(genpath(\""+addpath+"\")), [odg, MOV]=PQevalAudio(\""+wav_clean_name+"\",\""+wav_output_name+"\"), save(\""+odgmatfile_pred+"\",\"odg\"), save(\"mov.mat\",\"MOV\") , exit'"
|
||||
p2 = subprocess.Popen(bashCommand, stdout=subprocess.PIPE, shell=True)
|
||||
(output, err) = p2.communicate()
|
||||
|
||||
print(output)
|
||||
p1.wait()
|
||||
p2.wait()
|
||||
#I save the odg results in a .mat file, which I load here. Not the most optimal method, sorry :/
|
||||
annots_noise = loadmat(odgmatfile_noisy)
|
||||
annots_pred = loadmat(odgmatfile_pred)
|
||||
#Consider loading also the movs!!
|
||||
return annots_noise["odg"][0][0], annots_pred["odg"][0][0]
|
||||
|
||||
def inference(self, name, method=None):
|
||||
print("Inferencing :",name)
|
||||
if self.dataset_test!=None:
|
||||
if method=="EM":
|
||||
return self.inference_inner_classical(name, "EM")
|
||||
elif method=="wiener":
|
||||
return self.inference_inner_classical(name, "wiener")
|
||||
elif method=="wiener_declick":
|
||||
return self.inference_inner_classical(name, "wiener_declick")
|
||||
elif method=="EM_declick":
|
||||
return self.inference_inner_classical(name, "EM_declick")
|
||||
else:
|
||||
return self.inference_inner(name)
|
||||
Loading…
Reference in New Issue
Block a user