ГОСТ Р 56047-2014 Системы охранные телевизионные. Компрессия оцифрованных аудиоданных. Классификация. Общие требования и методы оценки алгоритмов

Метрика качества

Диапазон значений метрик качества по классам алгоритмов компрессии оцифрованных аудиоданных

Класс III

Класс II

Класс I

Пиковое отношение сигнал/шум (PSNR), дБ

Менее 30

[30; 40]

Свыше 40

Коэффициент различия форм сигналов

Более 10

[10; 10]

Менее 10

Объективная оценка аудиоданных с точки зрения восприятия (PEAQ)

[-3,98; -2,3)*

[-2,3; -0,62]

(-0,62; 0,22]*

Примечание - Метрики качества отражают изменения оцифрованных аудиоданных (после их обработки алгоритмами компрессии и декомпрессии), которые могут оказать критическое влияние на возможность использования восстановленных аудиоданных для установления наличия звуковых сигналов, дифференциации звуков и речи.

Наименования выходной переменной модели

Описание

BandwidthRefB

Ширина полосы исходного аудиосигнала

BandwidthTestB

Ширина полосы восстановленного аудиосигнала

Total NMRB

Отношение уровня шума к порогу маскирования

WinModDiff1B

Оконная разница модуляций

ADBB

Среднее поблочное искажение

EHSB

Гармоническая структура ошибки

AvgModDiff1B

Средняя разница модуляции 1

AvgModDiff2B

Средняя разница модуляции 2

RmsNoiseLoudB

Громкость шума

MFPDB

Максимальная вероятность обнаружения искажения

RelDistFramesB

Относительное искажение фреймов

где =-0,198719;

=0,0550197;

=-0,00102438;

(А.95)

=5,05622·10;

=9,01033·10;

=5,95072.

=6,39468;

=1,71332;

Индекс переменой, i

Наименование переменной

Минимальное значение [i]

Максимальное значение [i]

0

BandwidthRefB

393,916656

921,0

1

BandwidthTestB

361,965332

881,131226

2

TotalNMRB

-24,045116

16,212030

3

WinModDiff1B

1,110661

107,137772

4

ADBB

-0,206623

2,886017

5

EHSB

0,074318

13,933351

6

AvgModDiff1B

1,113683

63,257874

7

AvgModDiff2B

0,950345

1145,018555

8

RmsNoiseLoudB

0,029985

14,819740

9

MFPDB

0,000101

1,0

10

RelDistFramesB

0,0

1,0

Индекс i

Вес [i, 0]

Вес [i, 1]

Вес [i, 2]

0

-0,502657

0,436333

1,219602

1

4,307481

3,246017

1,123743

2

4,984241

2,211189

-0,192096

3

0,0511056

-1,762424

4,331315

4

2,321580

1,789971

-0,754560

5

-5,303901

-3,452257

-10,814982

6

2,730991

-6,111805

1,519223

7

0,624950

-1,331523

-5,955151

8

3,10288

0,871260

-5,922878

9

-1,051468

-0,939882

-0,142913

10

-1,804679

-0,503610

0,620456

Bias

Bias [0]

Bias [1]

Bias [2]

-0,307594

-2,518254

0,654841

-2,207228

Индекс j

Вес [j]

0

-3,817048

1

4,017138

2

4,629582

function ODG = PQevalAudio (Fref, Ftest, StartS, EndS)

% Оценка качества аудиоданных с точки зрения восприятия (Perceptual evaluation of audio quality)

% - StartS - индекс значения, соответствующего началу сигнала.

% - EndS - индекс значения, соответствующего концу сигнала.

% глобальные переменные

global MOVC PQopt

% параметры

NF = 2048;

Nadv = NF/2;

Version = 'Basic';

% настройки

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

% вычислить количество значений и каналов для каждого входного файла

WAV(1) = PQwavFilePar (Fref);

WAV(2) = PQwavFilePar (Ftest);

% согласовать размеры файлов

PQ_CheckWAV (WAV);

if (WAV(1).Nframe ~= WAV(2).Nframe)

disp ('>>>Number of samples differ: using the minimum');

end

% границы данных

Nchan = WAV(1).Nchan;

[Starts, Fstart, Fend] = PQ_Bounds (WAV, Nchan, StartS, EndS, PQopt);

% фреймов PEAQ

Np = Fend - Fstart + 1;

if (PQopt.Ni < 0)

PQopt.Ni=ceil (Np/abs(PQopt.Ni));

end

% инициализация структуры MOV

MOVC = PQ_InitMOVC (Nchan, Np);

Nc = PQCB (Version);

for (j = 0:Nchan-1)

Fmem(j+1) = PQinitFMem (Nc, PQopt.PCinit);

end

is = 0;

for (i = -Fstart:Np-1)

% считать фрейм данных

xR = PQgetData (WAV(1), StartS(1) + is, NF); % Reference file

xT = PQgetData (WAV(2), StartS(2) + is, NF); % Test file

is = is + Nadv;

% обработка фрейма

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)

% вывести MOV в новую структуру

PQframeMOV (i, MOVI); % выходные значения теперь в глобальной переменной MOVC

% вывод значений

if (PQopt.Ni ~= 0 & mod (i, PQopt.Ni) == 0)

%

PQprtMOVCi (Nchan, i, MOVC);

end

end

end

% усреднение по времени значений MOV

if (PQopt.DelayOverlap)

Nwup = Fstart;

else

Nwup = 0;

end

MOVB = PQavgMOVB (MOVC, Nchan, Nwup);

% запуск нейронной сети

ODG = PQnNet (MOVB);

% совокупный вывод значений

% PQprtMOV (MOVB, ODG);

%----------

function PQ_CheckWAV (WAV)

% проверка файлов

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;

% границы данных

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

% номер первого фрейма

Fstart = floor ((Lim(1) - StartS(1)) / PQ_NADV);

% номер последнего фрейма

Fend = floor ((Lim(2) - StartS(1) + 1 - PQopt.EndMin) / PQ_NADV);

%----------

function MOVC = PQ_lnitMOVC (Nchan, Np)

MOVC.MDiff.Mt1 В = 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);

function ODG = PQnNetB (MOV) '

% нейронная сеть для получения конечного значения метрики

persistent amin amax wx wxb wy wyb bmin bmax I J CLIPMOV

global PQopt

if (isempty (amin))

I = length (MOV);

if (l == 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

% нейтронная сеть

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

Dl = Dl + wy(j+1) * sigmoid (arg);

end

ODG = bmin + (bmax - bmin) * sigmoid (Dl);

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 = ...

[ 2.166.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

function [Nc, fc, fl, fu, dz] = PQCB (Version)

% параметры критической полосы пропускания

B = inline (7 * asinh (f / 650)');

Bl = inline ('650 * sinh (z / 7)');

fL = 80;

fU = 18000;

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 = Bl (zl);

fc = Bl (zc);

fu = Bl (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,

578434

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

function Es = PQspreadCB (E. Ver)

% распространение возбуждений

% E и Es - энергии

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

e = 0.4;

if (~ strcmp (Ver, Version))

Version = Ver;

[Nc, fc, fl, fu, dz] = PQCB (Version);

end

% вычисление памяти

aUCEe = zeros (1, Nc);

Ene = zeros (1, Nc);

Es = zeros (1, Nc);

% вычисление термов, зависящих от энергии

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);

glL = (1 -aL^ (m+1))/(1 - aL);

glU = (1 - aUCE^ (Nc-m)) / (1 - aUCE);

En = E(m+1) / (glL + glU - 1);

aUCEe(m+1) = aUCE^e;

Ene(m+1) = En^e;

end

% распространение вниз

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

% распространение вверх 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/е) / Bs(i+1);

end

function Eb = PQgroupCB (X2, Ver)

% группировка вектора энергии ДПФ и критической полосы пропускания

% Х2 - вектор значений, возведенных в степень 2

% Eb - вектор возбуждений

persistent Nc kl ku Ul Uu Version

Emin = 1e-12;

if (~ strcmp (Ver, Version))

Version = Ver;

NF = 2048;

Fs = 48000;

[Nc, kl, ku, Ul, Uu] = PQ_CBMapping (NF, Fs, Version);

end

% выделение памяти

Eb = zeros (1, Nc);

% вычисление возбуждений в каждой полосе

for (i = 0:Nc-1)

Ea = UI(i+1) * X2(kl(i+1)+1);

for (k = (kl(i+1)+1):(ku(i+1)-1))

Ea = Ea + X2(k+1);

end

Ea = Ea + Uu(i+1) * X2(ku(i+1)+1);

Eb(i+1) = max(Ea, Emin);

end

%---------------------------------------

function [Nc, kl, ku, Ul, Uu] = PQ_CBMapping (NF, Fs, Version)

[Nc, fc, fl, fu] = PQCB (Version);

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;

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;

if(kl(i+1) == ku(i+1))

Uu(i+1) = 0;

else

Uu(i+1) = (min(fui, (k+0.5)*df)...

- max(fli, (k-0.5)*df)) / df;

end

break;

end

end

end

function [MOVI, Fmem] = PQeval (xR, xT, Fmem)

% PEAQ - обработка единичного фрейма

NF = 2048;

Version = 'Basic';

% оконное ДПФ

X2(1,:) = PQDFTFrame (xR);

X2(2,:) = PQDFTFrame (xT);

[EbN, Es] = PQ_excitCB (X2);

[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,:));

% адаптация паттернов возбуждения

[ЕР, Fmem.Adap] = PQadapt (Ehs, Fmem.Adap, Version, 'FFT');

% паттерны модуляции

[M, ERavg, Fmem.Env] = PQmodPatt (Es, Fmem.Env);

% громкость

MOVI.Loud.NRef = PQIoud (Ehs(1,:), Version, 'FFT');

MOVI.Loud.NTest = PQIoud (Ehs(2,:), Version, 'FFT');

% разница модуляций

MOVI.MDiff = PQmovModDiffB (M, ERavg);

% громкость шума

MOVI.NLoud.NL = PQmovNLoudB (M, EP);

% полоса пропускания

MOVI.BW = PQmovBW (X2);

% отношений шума к маскированию

MOVI.NMR = PQmovNMRB (EbN, Ehs(1,:));

% вероятность обнаружения

MOVI.PD = PQmovPD (Ehs);

% ошибка гармонической структуры

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

% выделение памяти

XwN2 = zeros (1, NF/2+1);

% фильтрация на основе модели внешнего и среднего уха

Xw2(1,:) = W2.* X2(1,1:NF/2+1);

Xw2(2,:) = W2.* X2(2,1:NF/2+1);

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

Eb(1,:) = PQgroupCB (Xw2(1,:), Version);

Eb(2,:) = PQgroupCB (Xw2(2,:), Version);

EbN = PQgroupCB (XwN2, Version);

E(1,:) = Eb(1,:) + EIN;

E(2,:) = Eb(2,:) + EIN;

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

Ehs = zeros (1, Nc);

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

function X2 = PQDFTFrame (x)

persistent hw

NF = 2048;

if (isempty (hw))

Amax = 32768;

fc= 1019.5;

Fs = 48000;

Lp = 92;

GL = PQ_GL (NF, Amax, fc/Fs, Lp);

hw = GL * PQHannWin (NF);

end

xw = hw.* x;

X = PQRFFT (xw, NF, 1);

X2 = PQRFFTMSq (X, NF);

%----------------------------------------

function GL = PQ_GL (NF, Amax, fcN, Lp)

W = NF - 1;

gp = PQ_gp (fcN, NF, W);

GL = 10^ (Lp / 20) / (gp * Amax/4 * W);

%----------

function gp = PQ_gp (fcN, NF, W)

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));

function Lim = PQdataBoundary (WAV, Nchan, StartS, Ns)

PQ_L = 5;

Amax = 32768;

NBUFF = 2048;

PQ_ATHR = 200 * (Amax / 32768);

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

Lim(2) = -1;

is = StartS;

while (is <= EndS)

Nf = min (EndS - is + 1, NBUFF);

ie = is + Nf - 1;

js = EndS - (ie - StartS + 1) + 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

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)

s = s + (abs (x(i+L-1 + 1)) - abs (x(i-1 + 1)));

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)

s = s + (abs (x(i-L+1+1)) - abs (x(i+1+1)));

if (s > Thr)

ie = i;

return

end

end

function x = PQgetData (WAV, i, N)

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

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

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)

Amax = 32768;

Nchan = WAV.Nchan;

х = 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 * (wavread (WAV.Fname, [i+1 i+Ns-1+1]))';

end

function hw = PQHannWin (NF)

hw = zeros (1, NF);

for(n = 0:NF-1)

hw(n+1) = 0.5 * (1 - cos(2 * pi * n / (NF-1)));

end

function Fmem = PQinitFMem (Nc, PCinit)

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;

function EIN = PQIntNoise (f)

N = length (f);

for (m = 0:N-1)

INdB = 1.456 * (f(m+1)/ 1000)^(-0.8);

EIN(m+1) = 10^(INdB/10);

End

function X = PQRFFT (x, N, ifn)

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)];

xl = [0x(N/2+1+1:N-1+1) 0];

x = complex ([xR xR(N/2-1+1:-1:1+1)], [xl - xl(N/2-1+1:-1:1+1)]);

X = real (ifft (x, N));

end

function X2 = PQRFFTMSq (X, N)

X2 = zeros (1, N/2+1);

Х2(0+1) = Х(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;

function [a, b] = PQtConst (t100, tmin, f, Fs)

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

function W2 = PQWOME (f)

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

function WAV = PQwavFilePar (File)

persistent iB

if (isempty (iB))

iB = 0;

else

iB = mod (iB + 1, 2);

end

[size WAV.Fs Nbit] = wavread (File, 'size');

WAV.Fname = File;

WAV.Nframe = size(1);

WAV.Nchan = size(2);

WAV.iB = iB;

PQgetData (WAV, 0, 0);

function MOV = PQavgMOVB (MOVC, Nchan, Nwup)

Fs = 48000;

NF = 2048;

Nadv=NF / 2;

Fss = Fs / Nadv;

tdel = 0.5;

tex = 0.050;

[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 = PQIoudTest (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 = Iog10 (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 * Iog10 (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);

L = floor (tavg * Fss);

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;

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

AvgModDiff1В = s / Nchan;

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);

s = 0;

for (j = 0:Nchan-1)

s = s + PQ_RMSAvg (NLoud.NL(j+1,Ndel+1:Np-1+1));

end

RmsNoiseLoudB = s / Nchan;

%-----------------------------------

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

%----------

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

%----------

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

%----------

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

%----------

function LinAvg = PQ_LinAvg (x)

N = length (x);

s = 0;

for (i = 0:N-1)

s = s + x(i+1);

end

LinAvg = s / N;

%----------

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

function PQframeMOV (i, MOVI)

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;

function Ndel = PQIoudTest (Loud)

[Nchan, Np] = size (Loud.NRef);

Thr = 0.1;

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)

Np = length (NRef);

it = Np;

for (i = 0:Np-1)

if (NRef(i+1) > Thr & NTest(i+1) > Thr)

it = i;

break;

end

end

function BW = PQmovBW (X2)

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

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

function MDiff - PQmovModDiffB (M, ERavg)

persistent Nc Ete

if (isempty (Nc))

e = 0,3;

[Nc, fc] = PQCB ('Basic');

Et = PQIntNoise (fc);

for (m = 0:Nc-1)

Ete(m+1) = E(m+1)^e;

end

end

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

MD1В = num1В / (offset1В + 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;

function NL = PQmovNLoudB (M, EP)

persistent Nc Et

if (isempty (Nc))

[Nc, fc] = PQCB ('Basic');

Et = PQIntNoise (fc);

end

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) - serf * 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)е - 1);

end

NL = (24 / Nc) * s;

if (NL < NLmin)

NL = 0;

end

function NMR = PQmovNMRB (EbN, Ehs)

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

function PD = PQmovPD (Ehs)

Nc = length (Ehs);

PD.p = zeros (1, Nc);

PD.q = zeros (1, Nc);

persistent с g d1 d2 bP bM

if (isempty (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 * Iog10 (Ehs(1,m+1));

EdBT = 10 * Iog10 (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);

PD.q(m+1) = abs (fix(edB)) / s;

end

function PQprtMOV (MOV, ODG)

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;

function PQprtMOVCi (Nchan, i, MOVC)

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.Mt1 B(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

function EHS = PQmovEHS (xR, xT, X2)

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)';

if (EnRef < EnThr & EnTest < EnThr)

EHS = -1;

return;

end

D = zeros (1, kmax);

for (k = 0:kmax-1)

D(k+1) = log (X2(2,k+1) / X2(1,k+1));

end

C = PQ_Corr(D, NL, M);

Cn = PQ_NCorr (C, D, NL, M);

Cnm = (1 / NL) * sum (Cn(1:NL));

Cw = Hw. * (Cn - Cnm);

% DFT

cp = PQRFFT (Cw, NL, 1);

c2 = PQRFFTMSq (cp, NL);

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 С = PQ_Corr (D, NL, M)

NFFT = 2 * NL;

D0 = [D(1:M)zeros(1,NFFT-M)];

D1 = [D(1:M+NL-1) zeros(1 ,NFFT-(M+NL-1))];

d0 = PQRFFT (DO, NFFT, 1);

d1 = PQRFFT (D1, NFFT, 1);

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)

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)

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;

function [EP, Fmem] = PQadapt (Ehs, Fmem, Ver, Mod)

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

EP = zeros (2, Nc);

R = zeros (2, Nc);

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

CL = (sn/sd)^2;

for (m = 0:Nc-1)

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

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

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);

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)

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

function Ntot = PQIoud (Ehs, Ver, Mod)

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);

с =1.07664;

else

[Nc, fc] = PQFB;

с =1.26539;

end

E0 = 1e4;

Et = PQ_enThresh (fc);

s = PQ_exlndex (fc);

for (m = 0:Nc-1)

Ets(m+1) = с * (Et(m+1) / (s(m+1) * Е0))^е;

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))^е -1);

sN = sN + max(Nm, 0);

end

Ntot = (24 / Nc) * sN;

%====================

function s = PQ_exlndex (f)

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)

N = length (f);

for (m = 0:N-1)

EtdB = 3.64 * (f(m+1) / 1000)^(-0.8);

Et(m+1) = 10^(EtdB/10);

end

function [M, ERavg, Fmem] = PQmodPatt (Es, Fmem)

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

M = zeros (2, Nc);

е = 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 (Ее - 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,:);

Файл: common.h

#define DEBUG

#define HANN 2048

#define BARK 109

#define DOUBLE

#if defined(DOUBLE)

#define module(x) fabs((double) x)

#define p(x,y) pow((double)x, (double)y)

#elif defined(LDOUBLE)

#define module(x) fabsl((long double) x)

#define p(x,y) powl((long double)x, (long double)y)

#endif

/************************ detprob ***********************/

#defineC1 1.0

/************************** end **************************/

/*********************** harmstruct **********************/

#defineAVGHANN

#define SKIPFRAME

#define GETMAX

#define Fup 18000.0

#define Flow 80.0

#define PATCH 1

/***************************** end **************************/

/**************************** peaqb ************************/

#define LOGVARIABLE

#ifdef LOGVARIABLE

#define LOGALLFRAMES

#endif

/****************************** end **************************/

struct processing {

double fftref[HANN/2];

double ffttest[HANN/2];

double ffteref[HANN/2];

double ffletest[HANN/2];

double fnoise[HANN/2];

double pptest[BARK];

double ppref[BARK];

double ppnoise[BARK];

double E2test[BARK];

double E2ref[BARK];

double Etest[BARK];

double Eref[BARK];

double Mref[BARK];

double Modtest[BARK];

double Modref[BARK];

};

Файл: peaqb.h

#define LOGRESULT "analized"

#ifdef DEBUG

#define LOGFILE "debugged.txt"

#endif

#define OPT_REF

0x01

#define OPT_TEST

0x02

#define THRESHOLDDELAY 0.050

#define AVERAGINGDEALAY 0.5

#define B(f) 7 * asinh((double)f /650)

#define Bl(z) 650 * sinh((double)z /7)

/* Function prototypes */

void fatalerr(char *,...);

void usage(char*);

void logvariable(const char*, double*, int);

void ProcessFrame(signed int*, signed int*, int, signed int *,

signed int*, int, int, int, int);

/* Prototypes end */

Файл: peaqb.c

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <stdarg.h>

#include <getopt.h>

#include <assert.h>

#include <math.h>

#include <fftw.h>

#include <common.h>

#include <wavedump.h>

#include <getframe.h>

#include <bandwidth.h>

#include <levpatadapt.h>

#include <moddiff.h>

#include <modulation.h>

#include <loudness.h>

#include <neural.h>

#include <nmr.h>

#include <detprob.h>

#include <energyth.h>

#include <harmstruct.h>

#include <boundary.h>

#include <critbandgroup.h>

#include <earmodelfft.h>

#include <noiseloudness.h>

#include <reldistframes.h>

#include <spreading.h>

#include <timespreading.h>

#include threshold.h>

#include <peaqb.h>

extern int errno;

char *fileref, *filetest;

double hannwindow[HANN];

double Etesttmpch1[BARK], Etesttmpch2[BARK], Ereftmpch1[BARK],

Ereftmpch2[BARK], Cffttmpch1[HANN/2], Cffttmpch2[HANN/2];

int delaytime1, delaytime2;

int count = 0;

int harmsamples = 1;

fftw_plan plan, plan2;

/* Bark Tables */

double *fL, *fC, *fU;

int bark;

struct levpatadaptin levinch1, Ievinch2;

struct modulationin modintestch1, modintestch2, modinrefch1, modinrefch2;

struct moddiffin moddiffinch1, moddiffinch2;

struct bandwidthout bandwidthch1, bandwidthch2;

struct outframes processed;

int

main(int argc, char *argv[])

{

signed int ch1ref[HANN];

signed int ch2ref[HANN];

signed int ch1test[HANN];

signed int ch2test[HANN];

int opt_line = 0;

int rateref, numchref, bitsampleref, Ipref;

int ratetest, numchtest, bitsampletest, Iptest;

int boundflag, totalframes = 0;

FILE *fpref, *fptest;

struct boundaryflag boundbe = {0, 0};

struct out oveRet;

/* Parse command line */

if (argc < 3)

usage(argv[0]);

{

int с = 0;

while ((c = getopt(argc, argv, "r:t:h")) != EOF)

switch (c) {

case 'h':

usage(argv[0]);

break;

case 'r':

opt_line |= OPT_REF;

fileref = optarg;

break;

case 't':

opt_line |= OPT_TEST;

filetest = optarg;

break;

}

}

/* Input control */

if (!(opt_line & OPT_REF) || !*fileref)

fatalerr ("err: -r/--reference <arg> required");

if (!(opt_line & OPT_TEST) || !*filetest)

fatalerr ("err: -t/--test <arg> required");

Ipref = LevelPression(fileref);

Iptest = LevelPression(filetest);

/* Init routines */

// make Hann Window (2.1.3)

{

int k;

for(k=0;k<HANN;k++)

hannwindow[k] = 0.5*sqrt((double)8/3)*

(1 - cos((double)2*M_PI*k/(HANN -1)));

}

// make Bark tables (2.1.5)

{

int k;

double zL, zU;

double *zl, *zc, *zu;

zL = B(Flow);

zU = B(Fup);

bark = ceil((zU-zL)/dz);

fL = (double *)malloc(bark * sizeof(double));

fC = (double *)malloc(bark * sizeof(double));

fU = (double *)malloc(bark * sizeof(double));

zl = (double *)malloc(bark * sizeof(double));

zc = (double *)malloc(bark * sizeof(double));

zu = (double *)malloc(bark * sizeof(double));

assert(fL != NULL && fC != NULL && fU != NULL && zl != NULL

&& zc != NULL && zu != NULL);

for(k=0;k<bark;k++) {

zl[k] = zL + k*dz;

zu[k] = zL + (k+1)*dz;

zc[k] = 0.5*(zl[k] + zu[k]);

}

zu[bark-1] = zU;

zc[bark - 1] = 0.5 * (zl[bark-1] + zu[bark-1]);

for(k=0;k<bark;k++) {

fL[k] = Bl(zl[k]);

fU[k] = Bl(zu[k]);

fC[k] = Bl(zc[k]);

}

free(zl);

free(zu);

free(zc);

}

// Initialize temp var

memset(&levinch1, 0x00, sizeof(struct levpatadaptin));

memset(&levinch2, 0x00, sizeof(struct levpatadaptin));

memset(&modintestch1, 0x00, sizeof(struct modulationin));

memset(&modintestch2, 0x00, sizeof(struct modulationin));

memset(&modinrefch1, 0x00, sizeof(struct modulationin));

memset(&modinrefch2, 0x00, sizeof(struct modulationin));

memset(Etesttmpch1, 0x00, BARK * sizeof(double));

memset(Etesttmpch2, 0x00, BARK * sizeof(double));

memset(Ereftmpch1, 0x00, BARK * sizeof(double));

memset(Ereftmpch2, 0x00, BARK * sizeof(double));

memset(Cffttmpch1, 0x00, (HANN/2) * sizeof(double));

memset(Cffttmpch2, 0x00, (HANN/2) * sizeof(double));

memset(&moddiffinch1, 0x00, sizeof(struct moddiffin));

memset(&moddiffinch2, 0x00, sizeof(struct moddiffin));

memset(&bandwidthch1, 0x00, sizeof(struct bandwidthout));

memset(&bandwidthch2, 0x00, sizeof(struct bandwidthout));

// ref file

if ((fpref = fopen(fileref,"r")) == NULL)

fatalerr("err: %s", strerror(errno));

if ((rateref = SampleRate(fpref)) == -1)

fatalerr("err: error in WaveHeader");

if ((numchref = NumOfChan(fpref)) == -1)

fatalerr("err: error in WaveHeader");

if ((bitsampleref = BitForSample(fpref)) == -1)

fatalerr("err: error in WaveHeader");

if (FindData(fpref) == -1)

fatalerr("err: can't find Data Field");

// test file

if ((fptest = fopen(filetest,"r")) == NULL)

fatalerr("err: %s", strerror(errno));

if ((ratetest = SampleRate(fptest)) == -1)

fatalerr("err: error in WaveHeader");

if ((numchtest = NumOfChan(fptest)) == -1)

fatalerr("err: error in WaveHeader");

if ((bitsampletest = BitForSample(fptest)) == -1)

fatalerr("err: error in WaveHeader");

if(FindData(fptest) == -1)

fatalerr("err: can't find Data Field");

fprintf(stdout,"\n PEAQb Algorithm. Author Giuseppe Gottardi 'oveRet'"

" <gottardi@ailinux.org>\n");

fprintf(stdout,"\nRef File %s"

"\n - Sample Rate: %d"

"\n - Number Of Channel: %d"

"\n - Bits for Sample: %d"

"\n - Level Playback: %d\n\n", fileref, rateref,

numchref, bitsampleref, Ipref);

fprintf(stdout,"\nTest File %s"

"\n - Sample Rate: %d"

"\n - Number Of Channel: %d"

"\n - Bits for Sample: %d"

"\n - Level Playback: %d\n\n", filetest, ratetest,

numchtest, bitsampletest, Iptest);

// Processing

if(ratetest != rateref)

fatalerr("err: Can't process Wave Files with different Sample Rate");

if(numchref != numchtest)

fatalerr("err: Can't process Mono Wave with Stereo Wave");

// Find delaytimel for Loudness Threshold

delaytime1 = ceilf((float)THRESHOLDDELAY*ratetest*2/HANN);

// Find delaytime2 for Delayed Averaging

delaytime2 = ceilf((float)AVERAGINGDEALAY*ratetest*2/HANN);

// make fft plan

plan = fftw_create_plan(HANN, FFTW_FORWARD, FFTW_MEASURE);

while(harmsamples < (Fup/ratetest)*(HANN/2.0)/2.0)

harmsamples *= 2;

plan2 = fftw_create_plan(harmsamples, FFTW_FORWARD, FFTW_MEASURE);

if(numchref == 1) {

if (fseek(fpref, (HANN/2)*bitsampleref/8, SEEK_CUR) == -1)

fatalerr("err: %s", strerror(errno));

if (fseek(fptest, (HANN/2)*bitsampletest/8, SEEK_CUR) == -1)

fatalerr("err: %s", strerror(errno));

#ifdef DATABOUND_BE

#undef DATABOUND_ONE

{

int i = 0, flag = 0, f1, f2;

long dataref, datatest, br1, br2;

dataref = ftell(fpref);

datatest = ftell(fptest);

while(1) {

br1 = ftell(fpref);

br2 = ftell(fptest);

f1 = GetMonoFrame(fpref, (signed int *)ch1ref,

bitsampleref/8, HANN);

f2 = GetMonoFrame(fptest, (signed int*)ch1test,

bitsampletest/8, HANN);

if(f1 && f2) {

totalframes++;

if(boundary(ch1ref, ch1test, NULL, NULL, HANN) && !flag) {

boundbe. begin = totalframes;

flag = 1;

}

}

else {

fseek(fptest, br1, SEEK_SET);

fseek(fpref, br2, SEEK_SET);

break;

}

}

fseek(fptest, -(HANN/2)*bitsampletest/8, SEEK_CUR);

fseek(fpref, -(HANN/2)*bitsampleref/8, SEEK_CUR);

while(i<totalframes) {

GetMonoFrame(fpref, (signed int *)ch1ref, bitsampleref/8, HANN);

GetMonoFrame(fptest, (signed int *)ch1 test, bitsampletest/8, HANN);

fseek(fptest, -2*(HANN/2)*bitsampletest/8, SEEK_CUR);

fseek(fpref, -2*(HANN/2)*bitsampleref/8, SEEK_CUR);

i++;

if(boundary(ch1ref, ch1test, NULL, NULL, HANN)) {

boundbe.end = totalframes-i;

break;

}

}

fseek(fptest, datatest, SEEK_SET);

fseek(fpref, dataref, SEEK_SET);

}

#endif

while (GetMonoFrame(fpref, (signed int *)ch1ref,

bitsampleref/8, HANN)

&& GetMonoFrame(fptest, (signed int *)ch1test,

bitsampletest/8, HANN)) {

count++;

#ifdefDATABOUND_BE

if(count >= boundbe.begin && count <= boundbe.end)

boundflag = 1;

else

boundflag = 0;

#else

boundflag = boundary(ch1ref, ch1test, NULL, NULL, HANN);

#ifdef DATABOUND_ONE

{

static int flag 1 =0, flag2 = 0;

if(boundflag && !flag1)

flag1 = 1;

if(!boundflag && flag1)

flag2 = 1;

if(flag2)

boundflag = 0;

}

#endif

#endif

ProcessFrame((signed int *)ch1ref,

(signed int*)NULL, Ipref,

(signed int *)ch1test,

(signed int *)NULL,

Iptest, rateref, boundflag, HANN);

oveRet = neural(processed);

fprintf(stdout,"\nframe: %d"

#ifdef DATABOUND_BE

"/%d"

"\ndata boundary: %d -> %d"

#endif

"\nBandwidthRefb: %g"

"\nBandwidthTestb: %g"

"\nTotalNMRb %g"

"\nWinModDiff1b: %g"

"\nADBb: %g"

"\nEHSb: %g"

"\nAvgModDiff1b: %g"

"\nAvgModDiff2b %g"

"\nRmsNoiseLoudb: %g"

"\nMFPDb: %g"

"\nRelDistFramesb: %g"

"\nDI: %g"

"\nODG: %g\n",

count,

#ifdef DATABOUND_BE

totalframes, boundbe.begin, boundbe.end,

#endif

processed. BandwidthRefb,

processed.BandwidthTestb, processed.TotalNMRb,

processed.WinModDiff1b, processed.ADBb,

processed.EHSb, processed.AvgModDiff1b,

processed.AvgModDiff2b,

processed.RmsNoiseLoudb, processed.MFPDb,

processed.RelDistFramesb,

oveRet.Dl, oveRet.ODG);

}

{

FILE *res;

res = fopen(LOGRESULT,"a+");

fprintf(res,"\nFile: %s\n"

"\nframe: %d"

"\nBandwidthRefb: %g"

"\nBandwidthTestb: %g"

"\nTotalNMRb %g"

"\nWinModDiff1b: %g"

"\nADBb: %g"

"\nEHSb: %g"

"\nAvgModDiff1b: %g"

"\nAvgModDiff2b %g"

"\nRmsNoiseLoudb: %g"

"\nMFPDb: %g"

"\nRelDistFramesb: %g"

"\nDI: %g"

"\nODG: %g\n",

filetest, count, processed.BandwidthRefb,

processed.BandwidthTestb, processed.TotalNMRb,

processed.WinModDiff1b, processed.ADBb,

processed.EHSb, processed.AvgModDiff1b,

processed.AvgModDiff2b,

processed.RmsNoiseLoudb, processed.MFPDb,

processed.RelDistFramesb,

oveRet.DI, oveRet.ODG);

fclose(res);

}

}

if(numchref == 2) {

if (fseek(fpref, HANN*bitsampleref/8, SEEK_CUR) == -1)

fatalerr("err: %s", strerror(errno));

if (fseek(fptest, HANN*bitsampletest/8, SEEK_CUR) == -1)

fatalerr("err: %s", strerror(errno));

#ifdefDATABOUND_BE

#undef DATABOUND_ONE

{

int i = 0, flag = 0, f1, f2;

long dataref, datatest, br1, br2;

dataref = ftell(fpref);

datatest = ftell(fptest);

while(1) {

br1 = ftell(fpref);

br2 = ftell(fptest);

f1 = GetStereoFrame(fpref, (signed int *)ch1ref,

(signed int *)ch2ref, bitsampleref/8, HANN);

f2 = GetStereoFrame(fptest, (signed int *)ch1test,

(signed int *)ch2test, bitsampletest/8, HANN);

if(f1 && f2) {

totalframes++;

if(boundary(ch1ref, ch1test, ch2ref, ch2test, HANN) && Iflag) {

boundbe.begin = totalframes;

flag = 1;

}

}

else {

fseek(fptest, br1, SEEK_SET);

fseek(fpref, br2, SEEK_SET);

break;

}

}

fseek(fptest, -HANN*bitsampletest/8, SEEK_CUR);

fseek(fpref, -HANN*bitsampleref/8, SEEK_CUR);

while(i<totalframes) {

GetStereoFrame(fpref, (signed int *)ch1ref,

(signed int *)ch2ref, bitsampleref/8, HANN);

GetStereoFrame(fptest, (signed int *)ch1test,

(signed int *)ch2test, bitsampletest/8, HANN);

fseek(fptest, -2*HANN*bitsampletest/8, SEEK_CUR);

fseek(fpref, -2*HANN*bitsampleref/8, SEEK_CUR);

i++;

if(boundary(ch1ref, ch1test, ch2ref, ch2test, HANN)) {

boundbe.end = totalframes-i+1;

break;

}

}

fseek(fptest, datatest, SEEK_SET);

fseek(fpref, dataref, SEEK_SET);

}

#endif

while (GetStereoFrame(fpref, (signed int *)ch1ref,

(signed int *)ch2ref, bitsampleref/8, HANN)

&& GetStereoFrame(fptest, (signed int *)ch1test,

(signed int *)ch2test, bitsampletest/8, HANN)) {

count++;

#ifdef DATABOUND_BE

if(count >= boundbe.begin && count <= boundbe.end)

boundflag = 1;

else

boundflag = 0;

#else

boundflag = boundary(ch1ref, ch1test, ch2ref, ch2test, HANN);

#ifdef DATABOUND_ONE

{

static int flag1 =0, flag2 = 0;

if(boundflag && !flag1)

flag1 = 1;

if(!boundflag && flag1)

flag2 = 1;

if(flag2)

boundflag = 0;

}

#endif

#endif

ProcessFrame((signed int *)ch1ref,

(signed int*)ch2ref, Ipref,

(signed int *)ch1test,

(signed int *)ch2test,

Iptest, rateref, boundflag, HANN);

oveRet = neural(processed);

fprintf(stdout,"\nframe: %d"

#ifdef DATABOUND_BE

"/%d"

"\ndata boundary: %d -> %d"

#endif

"\nBandwidthRefb: %g"

"\nBandwidthTestb: %g"

"\nTotalNMRb %g"

"\nWinModDiff1b: %g"

"\nADBb: %g"

"\nEHSb: %g"

"\nAvgModDiff1b: %g"

"\nAvgModDiff2b %g"

"\nRmsNoiseLoudb: %g"

"\nMFPDb: %g"

"\nRelDistFramesb: %g"

"\nDI: %g"

"\nODG: %g\n",

count,

#ifdefDATABOUND_BE

totalframes, boundbe.begin, boundbe.end,

#endif

processed. BandwidthRefb,

processed.BandwidthTestb, processed.TotalNMRb,

processed.WinModDiff1b, processed.ADBb,

processed.EHSb, processed.AvgModDiff1b,

processed. AvgModDiff2b,

processed.RmsNoiseLoudb, processed.MFPDb,

processed. RelDistFramesb,

oveRet.DI, oveRet.ODG);

}

{

FILE *res;

res = fopen(LOGRESULT,"a+");

fprintf(res,"\nFile: %s\n"

"\nframe: %d"

"\nBandwidthRefb: %g"

"\nBandwidthTestb: %g"

"\nTotalNMRb %g"

"\nWinModDiff1b: %g"

"\nADBb: %g"

"\nEHSb: %g"

"\nAvgModDiff1b: %g"

"\nAvgModDiff2b %g"

"\nRmsNoiseLoudb: %g"

"\nMFPDb: %g"

"\nRelDistFramesb: %g"

"\nDI: %g"

"\nODG: %g\n",

filetest, count, processed.BandwidthRefb,

processed.BandwidthTestb, processed.TotalNMRb,

processed.WinModDiff1b, processed.ADBb,

processed.EHSb, processed.AvgModDiff1b,

processed. AvgModDiff2b,

processed.RmsNoiseLoudb, processed.MFPDb,

processed. RelDistFramesb,

oveRet.DI, oveRet.ODG);

fclose(res);

}

}

fftw_destroy_plan(plan);

fclose(fpref);

fclose(fptest);

return 0;

}

void

ProcessFrame(signed int *ch1ref, signed int*ch2ref, int Ipref,

signed int *ch1test, signed int*ch2test, int Iptest,

int rate, int boundflag, int hann)

{

int k;

static int ch = 1;

double Ntotaltest, Ntotalref;

struct levpatadaptout lev;

struct moddiffout mod;

struct processing processch1, processch2;

earmodelfft(ch1ref, Ipref, hann, processch1.ffteref,

processch1.fftref);

earmodelfft(ch1test, Iptest, hann, processch1.fltetest,

processch1.ffttest);

critbandgroup(processch1.ffteref, rate, hann, processch1. ppref);

AddlntNoise(processch1.ppref);

critbandgroup(processch1.fltetest, rate, hann, processch1.pptest);

AddlntNoise(processch1.pptest);

for(k=0;k<hann/2;k++)

processch1.fnoise[k] = module(processch1.flteref[k])

- module(processch1.fftetest[k]);

critbandgroup(processch1.fnoise, rate, hann, processch1.ppnoise);

spreading(processch1.pptest, processch1.E2test);

spreading(processch1.ppref, processch1.E2ref);

timespreading(processch1.E2test, Etesttmpch1, rate, processch1.Etest);

timespreading(processch1.E2ref, Ereftmpch1, rate, processch1.Eref);

threshold(processch1.Eref, processch1.Mref);

modulation(processch1.E2test, rate, &modintestch1, processch1.Modtest);

modulation(processch1.E2ref, rate, &modinrefch1, processch1.Modref);

// Data boundary

if(boundflag) {

static int countboundary = 1;

static double RelDistFramesb = 0, nmrtmp = 0;

bandwidth(processch1.ffttest, processch1.fftref, hann,

&bandwidthch1);

processed.BandwidthRefb = bandwidthch1.BandwidthRefb;

processed.BandwidthTestb = bandwidthch1.BandwidthTestb;

processed.TotalNMRb = nmr(processch1.ppnoise, processch1.Mref,

&nmrtmp, countboundary);

processed.RelDistFramesb = reldistframes(processch1.ppnoise,

processch1.Mref,

&RelDistFramesb,

countboundary);

countboundary++;

// Data boundary + Energy threshold

if(energyth(ch1test, ch1ref, hann)) {

static int countenergy = 1;

static double EHStmp = 0;

processed.EHSb = harmstruct(processch1. ffttest,

processch1.fftref,

&EHStmp, rate, Cffttmpch1,

harmsamples, &countenergy);

countenergy++;

}

}

// Delayed Averaging

if(count > delaytime2) {

static double nltmp = 0;

static int noise = 0, internal_count = 0, loudcounter = 0;

mod = moddiff(processch1.Modtest, processch1.Modref,

(double *)&(modinrefch1.Etildetmp));

processed.WinModDiff1b = ModDiff1(mod, &moddiffinch1,

count - delaytime2);

processed.AvgModDiff1b = ModDiff2(mod, &moddiffinch1);

processed.AvgModDiff2b = ModDiff3(mod, &moddiffinch1);

Ntotaltest = loudness(processch1.Etest);

Ntotalref = loudness(processch1.Eref);

if(Ntotaltest > 0.1 || Ntotalref > 0.1) {

noise = 1;

#ifdefined(LOUDMODO2)

internal_count = 0;

#endif

}

// Delayed Averaging + loudness threshold

if(noise && internal_count <= delaytime1) {

// skip 0.05 sec (about 3 frames)

internal_count++;

loudcounter++;

}

else {

lev= levpatadapt(processch1.Etest, processch1.Eref, rate,

&levinch1, hann);

processed.RmsNoiseLoudb = noiseloudness(processch1.Modtest,

processch1.Modref,

lev, &nltmp,

count - delaytime2

- loudcounter);

}

}

/*{

extern double Cfft[];

extern int maxk;

FILE *fp;

logvariable("Cfftsx.txt", Cfft, 128);

fp = fopen("Cfftsxmaxpos.txt", "a+");

fprintf(fp,"%d\n", maxk);

fclose(fp);

}*/

if(*ch2ref && *ch2test) {

ch=2;

earmodelfft(ch2ref, Ipref, hann, processch2.ffteref,

processch2.fftref);

earmodelfft(ch2test, Iptest, hann, processch2.fftetest,

processch2.ffttest);

critbandgroup(processch2.ffteref, rate, hann, processch2.ppref);

AddlntNoise(processch2.ppref);

critbandgroup(processch2.fftetest, rate, hann, processch2.pptest);

AddlntNoise(processch2.pptest);

for(k=0;k<hann/2;k++)

processch2.fnoise[k] = module(processch2.ffteref[k])

- module(processch2.fftetest[k]);

critbandgroup(processch2.fnoise, rate, hann, processch2.ppnoise);

spreading(processch2.pptest, processch2.E2test);

spreading(processch2.ppref, processch2.E2ref);

timespreading(processch2.E2test, Etesttmpch2, rate,

processch2.Etest);

timespreading(processch2.E2ref, Ereftmpch2, rate,

processch2.Eref);

threshold(processch2.Eref, processch2.Mref);

modulation(processch2.E2test, rate, &modintestch2,

processch2.Modtest);

modulation(processch2.E2ref, rate, &modinrefch2,

processch2.Modref);

// Data boundary

if(boundflag) {

static int countboundary = 1;

static double RelDistFramesb = 0, nmrtmp = 0;

bandwidth(processch2.ffttest, processch2.fftref, hann,

&bandwidthch2);

processed.BandwidthRefb += bandwidthch2.BandwidthRefb;

processed.BandwidthTestb += bandwidthch2.BandwidthTestb;

processed.BandwidthRefb /= 2.0;

processed.BandwidthTestb /= 2.0;

processed.TotalNMRb += nmr(processch2.ppnoise,

processch2.Mref, &nmrtmp,

countboundary);

processed. RelDistFramesb += reldistframes(processch2.ppnoise,

processch2.Mref,

&RelDistFramesb,

countboundary);

processed.TotalNMRb /= 2.0;

processed.RelDistFramesb /= 2.0;

countboundary++;

// Data boundary + Energy threshold

if(energyth(ch2test, ch2ref, hann)) {

static int countenergy = 1;

static double EHStmp = 0;

processed.EHSb += harmstruct(processch2.ffttest,

processch2.fftref,

&EHStmp, rate, Cffttmpch2,

harmsamples, &countenergy);

processed.EHSb/=2.0;

countenergy++;

}

}

// Delayed Averaging

if(count > delaytime2) {

static double nltmp = 0;

static int noise = 0, internal_count = 0, loudcounter = 0;

mod = moddiff(processch2.Modtest, processch2.Modref,

(double *)&(modinrefch2.Etildetmp));

processed.WinModDiff1b += ModDiff1(mod, &moddiffinch2,

count - delaytime2);

processed.AvgModDiff1b += ModDiff2(mod, &moddiffinch2);

processed.AvgModDiff2b += ModDiff3(mod, &moddiffinch2);

processed.WinModDiff1b /= 2.0;

processed.AvgModDiff1b /= 2.0;

processed.AvgModDiff2b /= 2.0;

Ntotaltest = loudness(processch2.Etest);

Ntotalref = loudness(processch2.Eref);

if(Ntotaltest> 0.1 || Ntotalref > 0.1) {

noise = 1;

#ifdefined(LOUDMODO2)

internal_count = 0;

#endif

}

// Delayed Averaging + loudness threshold

if(noise && internal_count <= delaytime1) {

// skip 0.05 sec (about 3 frames)

internal_count++;

loudcounter++;

}

else {

lev= levpatadapt(processch2.Etest, processch2.Eref, rate,

&levinch2, hann);

processed.RmsNoiseLoudb += noiseloudness(processch2.Modtest,

processch2.Modref,

lev, &nltmp,

count - delaytime2

- loudcounter);

processed.RmsNoiseLoudb /= 2.0;

}

}

}

{

static int ndistorcedtmp = 0;

static double Ptildetmp = 0, PMtmp = 0, Qsum = 0;

if(ch == 2)

processed.ADBb = detprob(processch1.Etest, processch2.Etest,

processch1.Eref, processch2.Eref,

&Ptildetmp, &PMtmp, &Qsum,

&ndistorcedtmp, hann);

else

processed.ADBb = detprob(processch1.Etest, NULL,

Processch1.Eref, NULL,

&Ptildetmp, &PMtmp, &Qsum,

&ndistorcedtmp, hann);

processed. MFPDb = PMtmp;

}

/*

#ifdef LOGVARIABLE

logvariable("fftetestsx.txt", processch1.fftetest, hann/2);

logvariable("ffterefsx.txt", processch1.ffteref, hann/2);

logvariable("ffttestsx.txt", processch1.ffttest, hann/2);

logvariable("fftrefsx.txt", processch1.fftref, hann/2);

logvariable("Etestsx.txt", processch1.Etest, bark);

logvariable("Erefsx.txt", processch1.Eref, bark);

logvariable("E2testsx.txt", processch1.E2test, bark);

logvariable("E2refsx.txt", processch1.E2ref, bark);

logvariable("pptestsx.txt", processch1.pptest, bark);

logvariable("pprefsx.txt", processch1.ppref, bark);

logvariable("ppnoisesx.txt", processch1.ppnoise, bark);

logvariable("Mrefsx.txt", processch1.Mref, bark);

logvariable("Modtestsx.txt", processch1.Modtest, bark);

logvariable("Modrefsx.txt", processch1.Modref, bark);

logvariable("fftetestdx.txt", processch2.fftetest, hann/2);

logvariable("ffterefdx.txt", processch2.ffteref, hann/2);

logvariable("ffttestdx.txt", processch2.ffttest, hann/2);

logvariable("fftrefdx.txt", processch2.fftref, hann/2);

logvariable("Etestdx.txt", processch2.Etest, bark);

logvariable("Erefdx.txt", processch2.Eref, bark);

logvariable("E2testdx.txt", processch2.E2test, bark);

logvariable("E2refdx.txt", processch2.E2ref, bark);

logvariable("pptestdx.txt", processch2.pptest, bark);

logvariable("pprefdx.txt", processch2.ppref, bark);

logvariable("ppnoisedx.txt", processch2.ppnoise, bark);

logvariable("Mrefdx.txt", processch2.Mref, bark);

logvariable("Modtestdx.txt", processch2.Modtest, bark);

logvariable("Modrefdx.txt", processch2.Modref, bark);

#endif

*/

/*{

extern double Cfft[];

extern int maxk;

FILE *fp;

logvariable("Cfftdx.txt", Cfft, 128);

fp = fopen("Cfftdxmaxpos.txt", "a+");

fprintf(fp,"%d\n", maxk);

fclose(fp);

}*/

return;

}

void

fatalerr(char * pattern,...)/* Error handling routine */

{

va_list ap;

va_start(ap, pattern);

fprintf(stderr,"PEAQ-");

vfprintf(stderr, pattern, ap);

fprintf(stderr," (exit forced).\n");

va_end(ap);

exit(-1);

}

#ifdef DEBUG

void

debug(char * pattern,...) /* Debug handling routine */

{

FILE log;

va_list ap;

va_start(ap, pattern);

log=fopen(LOGFILE,"a+");

vfprintf(log, pattern, ap);

va_end(ap);

fclose(log);

return;

}

#endif

#ifdef LOGVARIABLE

void

logvariable(const char *filename, double *var, int len)

{

FILE *fp;

int k;

#ifdef LOGALLFRAMES

fp = fopen(filename,"a+");

#else

fp = fopen(filename,"w");

#endif

for(k=0;k<len;k++)

fprintf(fp,"%g\n",var[k]);

fclose(fp);

return;

}

#endif

void

usage(char * name) /* Print usage */

{

fprintf(stderr, "PEAQ Algorithm. Giuseppe Gottardi 'oveRet'"

"<gottardi@ailinux.org>\n\n");

fprintf(stderr, "usage: %s <option>\n", name);

fprintf(stderr,"

-r reffile[:lp] (Ip default = 92)\n"

" -t testfile[:lp] (Ip default = 92)\n"

"-h print this help\n");

exit (0);

}

Файл: bandwidth.h

#define ZEROTHRESHOLD 921

#define BwMAX 346

struct bandwidthout {

double sumBandwidthRefb;

int countref;

double sumBandwidthTestb;

int counttest;

double BandwidthRefb;

double BandwidthTestb;

};

/* Function prototypes */

int bandwidth(double *, double *, int, struct bandwidthout *);

/* Prototypes end */

Файл: bandwidth.c

#include <stdlib.h>

#include <math.h>

#include <common.h>

#include <bandwidth.h>

int

bandwidth(double *ffttest, double *fftref, int hann,

struct bandwidthout *out)

{

int k, BwRef = 0, BwTest = 0;

double Flevtest, Flevref;

double ZeroThreshold;

ZeroThreshold = 20.0 * log10((double)ffttest[ZEROTHRESHOLD]);

for(k=ZEROTHRESHOLD;k<hann/2;k++) {

Flevtest = 20.0 * Iog10((double)ffttest[k]);

if(Flevtest > ZeroThreshold)

ZeroThreshold = Flevtest;

}

for(k=ZEROTHRESHOLD-1; k>=0;k--) {

Flevref = 20.0 * Iog10((double)fftref[k]);

if(Flevref >= 10.0+ZeroThreshold) {

BwRef = k + 1;

break;

}

}

for(k=BwRef-1;k>=0;k--) {

Flevtest = 20.0 * Iog10((double)ffttest[k]);

if(Flevtest >= 5.0+ZeroThreshold) {

BwTest = к + 1;

break;

}

}

if( BwRef >BwMAX) {

out->sumBandwidthRefb += (double)BwRef;

out->countref++;

}

if(BwTest > BwMAX) {

out->sumBandwidthTestb += (double)BwTest;

out->counttest++;

}

if(out->countref == 0)

out->BandwidthRefb = 0;

else

out->BandwidthRefb = out->sumBandwidthRefb/(double)out->countref;

if(out->counttest == 0)

out->BandwidthTestb = 0;

else

out->BandwidthTestb = out->sumBandwidthTestb/(double)out->counttest;

return 0;

}

Файл: boundary.h

#define BOUNDWIN 5

#defineBOUNDLIMIT 200

struct boundaryflag {

int begin;

int end;

};

/* Function prototypes */

int boundary(signed int *, signed int *, signed int *, signed int *, int);

/* Prototypes end */

Файл: boundary.c

#include <stdlib.h>

#include <math.h>

#include <common.h>

#include <boundary.h>

int

boundary(signed int *ch1ref, signed int *ch1test, signed int *ch2ref,

signed int *ch2test, int hann)

{

int k, i, sum;

int ch1t = 0, ch1r = 0, ch2t = 0, ch2r = 0;

for(k=0;k<hann-BOUNDWIN+1;k++) {

if(!ch1t) {

sum = 0;

for(i=0;i<BOUNDWIN;i++)

sum += abs(ch1test[k+i]);

if(sum > BOUNDLIMIT)

ch1t=1;