% funkcia na vypocet objektivnych parametrov akustickeho priestoru7
% vyuziva suboru octdsgn.m z toolboxu octave na navrh pasmovych priepusti 
% doba dozvuku je extrapolovana z poklesu z -5 na -35 dB
% Inter aural cross corelation IACCE je pocitane iba v pripade ze vstupny
% signal je stereo
% Vstupom je impulzna odozva priestoru, pripadne reverberatoru a vzorkovaci
% kmitocet

function [RT,EDT,C7,C50,C80,D50,H,itdg,BR,HFR,Ts,IACCE,EKdelta] = objective(x,Fvz)

% stredne frekvencie oktav
fstredne = [125 250 500 1000 2000 4000 8000];

%integracny cas IACC
t2=3528;

% rad filtru
N = 3;

% pocet kanalov
kanaly = size(x,2);

t7 = round(0.007*Fvz);
t50 = round(0.05*Fvz);
t80 = round(0.08*Fvz);

%vytvori sadu filtrov na prislusiacich strednych kmitoctoch
a = zeros(length(fstredne),(2*N)+1);
b = zeros(length(fstredne),(2*N)+1);

for i = 1:length(fstredne)
    [b(i,:),a(i,:)] = octdsgn(fstredne(i),Fvz);
end

% vytvorenie matic a vektorov na naplnenie
z = zeros([length(fstredne) size(x)]);
rt = zeros([length(fstredne) kanaly]);
edt = zeros([length(fstredne) kanaly]);
itdg = zeros(1,kanaly);

ct7 = zeros(1,kanaly);
ct50 = zeros(1,kanaly);
ct80 = zeros(1,kanaly);

%filtrovanie
for i = 1:kanaly
    %vyhlada maximalnu spicku odpovedajucu priamej vlne
    itdg(i) = find(x(:,i)>0,1,'first'); 
    
    for j = 1:length(fstredne)
        y = filter(b(j,:),a(j,:),x(:,i)); % vyfiltruje signal
        temp = cumtrapz(y(end:-1:1).^2); % vytvori krivku dozvuku
        z(j,:,i) = temp(end:-1:1);
                
        [rt(j,i)] = dozvuk(z(j,itdg:end,i),[-5,-35],Fvz); % doba dozvuku
        [edt(j,i)] = dozvuk(z(j,itdg:end,i),[0,-10],Fvz); % pocatecni doba dozvuku
        
        disp(j);
        ct7(j,i) = 10.*log(trapz(y(itdg(i):itdg(i)+t7).^2)/trapz(y(itdg(i)+t7+1:end).^2));
        ct50(j,i) = 10.*log(trapz(y(itdg(i):itdg(i)+t50).^2)/trapz(y(itdg(i)+t50+1:end).^2));
        ct80(j,i) = 10.*log(trapz(y(itdg(i):itdg(i)+t80).^2)/trapz(y(itdg(i)+t80+1:end).^2));
    end
        % hodnoty Ts,C7,C50,C80
        p=(itdg(i):length(x(:,i)));
        p=p/Fvz;
        ts(i)=  trapz((x(itdg(i):end,i).^2)'.*p(1:end))/trapz(x(itdg(i):end,i).^2);
       
end
    
%Inter aural cross corelation coeficient IACC
if size(x,2)==1
    IACCE='NaN'
    
else
    
xleft=x(min(itdg):t2+min(itdg),1);
xright=x(min(itdg):min(itdg)+t2,2);    
    
xleftsquare=x(min(itdg):t2+min(itdg),1).^2;
xrightsquare=x(min(itdg):t2+min(itdg),2).^2;

xleftsum = (sum(xleftsquare))/Fvz;
xrightsum =(sum(xrightsquare))/Fvz;

tau = round(0.001*Fvz);
xright = [zeros(tau,1); xright; zeros(tau,1)];

for i = -tau:tau
    num = xleft.*(xright(i+tau+1:i+tau+t2+1));
    IACF(i+tau+1) = (sum(num))/Fvz;
end

IACCE = max(IACF/sqrt(xleftsum*xrightsum));

end

%rt60 pre jednotlive oktavy
if kanaly~=1;
rt60 = (rt(:,1)+rt(:,2))/2;
else
rt60=rt;    
end

RT=sum(rt60)/length(fstredne);

%EDT pre jednotlive oktavy
if kanaly~=1;
EDT = (edt(:,1)+edt(:,2))/2;
else
EDT=edt;    
end

EDT=sum(EDT)/length(fstredne);

%initial time delay gap v ms
itdg=min(itdg)/Fvz*1000;    

%bass ratio
BR=(rt60(1)+rt60(2))/(rt60(3)+rt60(4));
%high frequency ratio
HFR=(rt60(3)+rt60(4))/(rt60(5)+rt60(6));

%c50
if kanaly~=1;
C50 = (ct50(3,1)+ct50(4,1)+ct50(3,2)+ct50(4,2))/4;
else
C50=(ct50(3,1)+ct50(4,1))/2;    
end

%c7
if kanaly~=1;
C7 = (ct7(3,1)+ct7(4,1)+ct7(3,2)+ct7(4,2))/4;
else
C7=(ct7(3,1)+ct7(4,1))/2;    
end

%c80
if kanaly~=1;
C80 = (ct80(3,1)+ct80(4,1)+ct80(3,2)+ct80(4,2))/4;
else
C80=(ct80(3,1)+ct80(4,1))/2;    
end


%d50
D50=(10^(0.1*C50))/(1+10^(0.1*C50));

%H
H=-C50;

%ts
if kanaly~=1;
Ts = ((ts(:,1)+ts(:,2))/2)*1000;
else
Ts=ts*1000;    
end

% Echo Criterion
deltate=0.014;   %koeficient - te=0.009 pro řeč, te=0.014 pro hudbu

IRek=(x.^2)';

IRek=IRek(1:length(IRek));

EKnum=IRek;
EKden=IRek;

sumnom(length(IRek))=0;
sumden(length(IRek))=0;

for i=1:length(IRek)
    EKnum(1,i)=EKnum(1,i).*(i/Fvz);
end

sumnom(1,1) = 0;
sumden(1,1) = IRek(1,1);

for i = 2:length(IRek)
sumnom(1,i)=sumnom(1,i-1)+EKnum(1,i);
sumden(1,i)=sumden(1,i-1)+EKden(1,i);
end

tau = round(deltate*Fvz);

EK=sumnom./sumden;
EKdelta(length(IRek))=0;

for i = 1:length(IRek)-tau
EKdelta(1,i) = EK(1,tau+i)- EK(1,1+i);
end

EKdelta = EKdelta./deltate;