close all;
clear all;
clc

% vyber Barkovho pasma pre zobrazenie v grafoch
Nbi = 6;

% pridanie cesty pre priecinok s odozvami
addpath('Databáza impulzných odoziev')

% vyber suboru s impulznou odozvou
[Subor,Cesta,Index] = uigetfile('Databáza impulzných odoziev\ *.wav','Výber impulznej odozvy');

% h - impulzna odozva
% fvz - vzorkovaci kmitocet
[h,fvz] = audioread(Subor);

% celkove trvanie odozvy
Tcelk = length(h)/fvz;

% casovy vektor vzoriek odozvy
Th_vekt = linspace(0,Tcelk,length(h));

% vstupne parametre pre spektrogram
DlzkaOkna = ones(512,1); % vektor poctu okien 
nOkno = 512; % pocet casti rozdelenia signalu
Prekrytie = 0.5; % prekrytie susediadich casti rozdelenia
TypOkna = 'hann';
nFFT = fvz/50; % pocet bodov diskretnej fourierovej transformacie

% vypocet spektrogramu
% S - matica casovo-kmitoctovej distribucie (spektrogramu)
% F - vektor kmitoctov odvodenych z vypoctu spektrogramu
% T - casovy vektor bodov spektrogramu
[S,F,T] = spectrogram(h,DlzkaOkna,nOkno*Prekrytie,2*nFFT,fvz);

% EDR v zavislosti na kmitocte

% zmensenie maximalnej hodnoty kmitoctov na 20000 Hz
[~,F_h_20k] = min(abs(F-20000)); % hodnota indexu kmitoctu 20 kHz
F_20k = F(1:F_h_20k); % vektor kmitoctov zo spektrogramu do 20 kHz

% spektrogram pre kmitocty do 20 kHz
S_20k = S(1:F_h_20k,1:size(S,2));

% vypocet echogramu
S2 = S_20k.*conj(S_20k);

% vypocet EDR
% F_rozd - pocet rozdeleni v kmitoctovej osi
% T_rozd - pocet rozdeleni v casovej osi
[F_rozd,T_rozd] = size(S2);
S2_edr = zeros(F_rozd,T_rozd); % vytvorenie matice EDR
for i=1:F_rozd
    S2_edr(i,:) = fliplr(cumsum(fliplr(S2(i,:)))); % EDR
end

% premena relativnej velkosti EDR na dB
S2_edr_db = 10*log10(abs(S2_edr/((1*10^-12))));

% normalizacia EDR na 0 dB
offset = max(max(S2_edr_db)); % hodnota odchylky najvacsej hodnoty EDR od 0 dB
S2_edr_db_norm = S2_edr_db-offset; % EDR s normalizovanymi hodnotami


% vykreslenie grafu
figure('Name',Subor);clf; % vytvorenie a pomenovanie okna s grafom
surf(T,F_20k,S2_edr_db_norm); % typ grafu a zadane veliciny ktore sa zobrazia
view(0,90); % zmena zobrazenia
title('Reliéf poklesu energie'); % nazov grafu
xlabel('{\itt} (s) \rightarrow') 
ylabel('{\itf} (Hz)  \rightarrow') % nazov osi
zlabel('{\itL} (dB) \rightarrow') 
colormap jet % zmena farebnej skaly grafu
axis tight; % zmena rozpatia osi
shading flat % nastavenie tienovania grafu
c=colorbar; % nazov osi farebnej skaly
c.Label.String = '{\itL} (dB) \rightarrow';
set (gca,'zlim',[-120 0]); % zmena rozpatia osi Z
set (gcf,'pos',[820 500 550 500]) % urcenie pozicie a rozmerov grafu pre vykreslenie na obrazovke

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% EDR podla barkovych pasiem

% spektrogram pre kmitocty do 15500 Hz
[~,F_h_15500] = min(abs(F-15500)); % hodnota indexu kmitoctu 15,5 kHz
S_15500 = S(1:F_h_15500,1:size(S,2)); % spektrogram pre kmitocty do 15,5 kHz

% vypocet echogramu pre kmitocty do 15500 Hz
S2_15500 = S_15500.*conj(S_15500);

% zmensenie maximalnej hodnoty kmitoctov na 15500 Hz
F_15500 = F(1:F_h_15500);

% vektor strednych kmitoctov barkovych pasiem
fS = [0 50 150 250 350 450 570 700 840 1000 1170 1370 1600 1850 2150 2500 2900 3400 4000 4800 5800 7000 8500 10500 13500];

% vektor hranicnych kmitoctov barkovych pasiem
fH = [0 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];
% vektor hranicnych kmitoctov barkovych pasiem pre pouzitie pri filtrovani pasiem
fH_filter = [20 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];

% vektor barkovych pasiem
Bark_vekt = linspace(1,length(fS)-1,length(fS)-1);

% casovy vektor barkovych pasiem
T_vekt = linspace(1,size(S2_15500,2),size(S2_15500,2));

% vektor obalky
Obalka_vekt = ones(1,length(F_15500));
for y = round(fH/50+1)
    Obalka_vekt(y) = 0.5; % vysledna obalka
end

% vypocet EDR
% F_rozd - pocet rozdeleni v kmitoctovej osi
% T_rozd - pocet rozdeleni v casovej osi
[F_rozd_bark,T_rozd_bark] = size(S2_15500);

S2_15500_e = zeros(F_rozd_bark,T_rozd_bark); % vytvorenie matice EDR
for i2=1:F_rozd_bark
    S2_15500_e(i2,:) = fliplr(cumsum(fliplr(S2_15500(i2,:)))); % EDR bez obalky
end

for cas = T_vekt
    for bark = Bark_vekt
            z = Obalka_vekt((round(fH(bark)/50+1)):(round(fH(bark+1)/50+1))); % obalka v jednotlivych pasmach
            S2_15500_ed(bark,cas) = z*S2_15500_e(round(fH(bark)/50+1):(round(fH(bark+1)/50+1)),cas);
            S2_15500_edr(bark,cas) = S2_15500_ed(bark,cas)/(length(z)-1); % EDR pre barkove pasma
    end
end

% premena relativnej velkosti EDR na dB
S2_15500_edr_db = 10*log10(abs(S2_15500_edr/((1*10^-12))));

% normalizacia EDR na 0 dB
offset_bark = max(max(S2_15500_edr_db)); % hodnota odchylky najvacsej hodnoty EDR od 0 dB
S2_15500_edr_db_norm (Bark_vekt,:) = S2_15500_edr_db-offset_bark; % normalizovane EDR

% vykreslenie grafu
figure('Name',Subor) % vytvorenie a pomenovanie okna s grafom
surf(T,Bark_vekt,S2_15500_edr_db_norm); % typ grafu a zadane veliciny ktore sa zobrazia
view(0,90); % zmena zobrazenia
title('Reliéf poklesu energie'); % nazov grafu
xlabel('{\itt} (s) \rightarrow') 
ylabel('{\itf} (Hz)  \rightarrow') % nazov osi
zlabel('{\itL} (dB) \rightarrow') 
colormap jet % zmena farebnej skaly grafu
axis tight; % zmena rozpatia osi
shading flat % nastavenie tienovania grafu
c=colorbar; % nazov osi farebnej skaly
c.Label.String = '{\itL} (dB) \rightarrow';
set (gca,'zlim',[-120 0]); % zmena rozpatia osi Z
set (gcf,'pos',[1370 500 550 500]) % urcenie pozicie a rozmerov grafu pre vykreslenie na obrazovke

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% krivky dozvuku v pasmach podla Barkovej stupnice

% zobrazenie Schroederovych integralov pre vsetky Barkove pasma
figure('Name',Subor) % vytvorenie a pomenovanie okna s grafom
title(['Schröderove integrály pre Barkove pásma']); % nazov grafu
hold on
for bark = 1:24
    plot(T,S2_15500_edr_db_norm(bark,:)) % typ grafu a zadane veliciny ktore sa zobrazia
end
set(gca,'ylim',[-90,0]) % zmena zobrazenia
xlabel('{\itt} (s) \rightarrow ') % nazov osi
ylabel('{\itL} (dB) \rightarrow')
grid on % zobrazenie mriezky
legend show % zobrazenie legendy

% prepocitanie kriviek dozvuku podla casovej osi povodnej impulznej odozvy
for bark = 1:24
% prva cast krivky dozvuku

% najdenie hodnot z ktorych sa vypocita prva cast krivky (vsetky pasma)
    S2_RT = find (S2_15500_edr_db_norm(bark,:)<S2_15500_edr_db_norm(bark,1) & S2_15500_edr_db_norm(bark,:)>S2_15500_edr_db_norm(bark,round(0.7/10*length(S2_15500_edr_db_norm))));

    Y = S2_15500_edr_db_norm(bark,S2_RT); % vektor hodnot z ktorych sa vypocita prva cast krivky

% casovy vektor tychto hodnot        
    T_vekt_RT = T(S2_RT);

% vypocet prvej casti krivky
    P = polyfit(T_vekt_RT,Y,1); % vytvorenie linearneho poklesu z vektoru hodnot z ktorych sa vypocita prva cast krivky
    L_RT = polyval(P,T_vekt_RT); % vektor hodnot linearneho poklesu

    A(bark) = (L_RT(2)-L_RT(1))/(T_vekt_RT(2)-T_vekt_RT(1)); % smernica poklesu hladiny akustickeho tlaku
    RT60(bark,:) = -60/A(bark); % doba dozvuku
    B(bark) = S2_15500_edr_db_norm(bark,1); % pociatocna hladina akustickeho tlaku prvej casti krivky
    
% predlzenie prvej casti krivky    
    T_vekt_RT60_h = Th_vekt; % casovy vektor pre predlzenie priamky

    T_vekt_RT60_index_h = linspace(1,length(T_vekt_RT60_h),length(T_vekt_RT60_h)); % vektor vzoriek priamky

% prva cast krivky dozvuku pre vsetky Barkove pasma
    L_vekt_RT60_h(bark,T_vekt_RT60_index_h) = A(bark)*T_vekt_RT60_h+B(bark);
%%%%%%
% druha cast krivky dozvuku
    
% najdenie hodnot z ktorych sa vypocita druha cast krivky (vsetky pasma)
    S2_RT_2 = find (S2_15500_edr_db_norm(bark,:)<S2_15500_edr_db_norm(bark,round(4.5/10*length(S2_15500_edr_db_norm))) & S2_15500_edr_db_norm(bark,:)>S2_15500_edr_db_norm(bark,round(7/10*length(S2_15500_edr_db_norm))));
    
    Y_2 = S2_15500_edr_db_norm(bark,S2_RT_2); % vektor hodnot z ktorych sa vypocita druha cast krivky
    
% casovy vektor tychto hodnot    
    T_vekt_RT_2 = T(S2_RT_2);

% vypocet druhej casti krivky 
    P_2 = polyfit(T_vekt_RT_2,Y_2,1); % vytvorenie linearneho poklesu z vektoru hodnot z ktorych sa vypocita druha cast krivky
    L_RT_2 = polyval(P_2,T_vekt_RT_2); % vektor hodnot linearneho poklesu

    A_2(bark) = (L_RT_2(2)-L_RT_2(1))/(T_vekt_RT_2(2)-T_vekt_RT_2(1)); % smernica poklesu hladiny akustickeho tlaku
    RT60_2 = -60/A_2(bark); % doba dozvuku
    B_2(bark) = (S2_15500_edr_db_norm(bark,(S2_RT_2(1)))-A_2(bark)*T(S2_RT_2(1))); % pociatocna hladina akustickeho tlaku druhej casti krivky

% predlzenie druhej casti krivky      
    T_vekt_RT60_2_h = Th_vekt; % casovy vektor pre predlzenie priamky

    T_vekt_RT60_index_2_h = linspace(1,length(T_vekt_RT60_2_h),length(T_vekt_RT60_2_h)); % vektor vzoriek priamky

% druha cast krivky dozvuku pre vsetky Barkove pasma
    L_vekt_RT60_2_h(bark,T_vekt_RT60_index_2_h) = A_2(bark)*T_vekt_RT60_2_h+B_2(bark);
%%%%%%
% vysledna krivka dozvuku

    L_vekt_krivka_h(bark,:) = L_vekt_RT60_2_h(bark,:) - L_vekt_RT60_h(bark,:); % definicia krivky ako rozdielu povodnych priamok
    C_h = [1,-A(bark);1,-A_2(bark)]; % cast sustavy dvoch rovnic o dvoch neznamych
    a_h = [B(bark);B_2(bark)]; % cast sustavy dvoch rovnic o dvoch neznamych
    yx_h(1:2,bark) = C_h\a_h; % najdenie priesecniku
    
    L_vekt_priesecnik_h(bark) = (B_2(bark)-B(bark))/(A(bark)-A_2(bark)); % cas v ktorom sa nachadza priesecnik
    [~,L_vekt_index_h] = min(abs(Th_vekt-L_vekt_priesecnik_h(bark))); % index postupnosti priesecnika
    p_h = L_vekt_index_h-1; % oprava indexu postupnosti priesecnika
    L_vekt_krivka_h(bark,1:p_h) = zeros(1,p_h); % vynulovanie casti pred priesecnikom
    L_vekt_krivka_vysl_db_h(bark,:) = L_vekt_krivka_h(bark,:)+L_vekt_RT60_h(bark,:); % vysledna krivka
end


% pridanie tabulky s udajmi z Excel suboru
[Subor_num,Subor_txt,Subor_raw] = xlsread('info.xlsx');

% ziskanie udaju o objeme miestnosti
Subor_nazov = strsplit(Subor,{'-','.'}); % ziskanie nazvu suboru pre vyhladanie v tabulke
Subor_info = Subor_num(find(strcmp(Subor_raw, Subor_nazov(1)))-2*length(Subor_raw)-1,:); % vyhladanie udajov o danej odozve z tabulky

% objem miestnosti uvedený v zdroji
V = Subor_info(4);

% Schroederov kmitocet podľa objemu uvedeneho v zdroji
f_sch = 2000*sqrt(RT60(9,:)/V);


% zobrazenie Schroederovho integralu a krivky dozvuku
figure('Name',Subor) % nazov okna s grafom
title(['Schröderov integrál a časti krivky dozvuku pre ',num2str(Nbi),'. Barkovo pásmo']); % nazov grafu
hold on

plot(T,S2_15500_edr_db_norm(Nbi,:)) % typ grafu a zobrazene veliciny
set(gca,'ylim',[-90,0]) % zmena zobrazenia
xlabel('{\itt} (s) \rightarrow ')
ylabel('{\itL} (dB) \rightarrow') % nazov osi
grid on % zobrazenie mriezky

plot(Th_vekt,L_vekt_RT60_h(Nbi,:),'-') % typ grafu a zobrazene veliciny
legend(['Schröderov integrál (',num2str(Nbi),'. Barkovo pásmo'],['Časti krivky dozvuku (',num2str(Nbi),'. Barkovo pásmo)']) % legenda

plot(Th_vekt,L_vekt_RT60_2_h(Nbi,:),'-') % typ grafu a zobrazene veliciny


% zobrazenie krivky dozvuku v casovej ose povodnej odozvy a Schroederovho integralu
figure('Name',Subor) % nazov okna s grafom
title(['Schröderov integrál a krivka dozvuku pre ',num2str(Nbi),'. Barkovo pásmo']); % nazov grafu
hold on
plot(T,S2_15500_edr_db_norm(Nbi,:)) % typ grafu a zobrazene veliciny
set(gca,'ylim',[-90,0]) % zmena zobrazenia
xlabel('{\itt} (s) \rightarrow ')
ylabel('{\itL} (dB) \rightarrow') % nazov osi
grid on % zobrazenie mriezky

plot(Th_vekt,L_vekt_krivka_vysl_db_h(Nbi,:),'-') % typ grafu a zobrazene veliciny
legend(['Schroederov integrál (',num2str(Nbi),'. Barkovo pásmo'],['Krivka dozvuku (',num2str(Nbi),'. Barkovo pásmo)']) % legenda


for bark = 1:24
% vytvorenie priebehu vyslednych kriviek dozvuku
    n2 = wgn(1,length(h),0); % excitacny signal
    e2(bark,:) = 10.^(L_vekt_krivka_vysl_db_h(bark,:)/20)*((1*10^-12)); % odlogaritmovana obalka
    
    out2(bark,:) = n2.*e2(bark,:); % signal s postupnostou podla obalky

% filtrovanie priebehov kriviek dozvuku pasmovymi priepustami pre Barkove pasma
    [filtA(bark,:),filtB(bark,:)] = butter(2,[fH_filter(bark)/fvz fH_filter(bark+1)/fvz]);

% vysledna namodelovana impulzna odozva    
    h_nove(bark,:) = filter(filtA(bark,:),filtB(bark,:),out2(bark,:));
end

% spojenie filtrovanych casti modelovanej odozvy v Barkovych pasmach do vyslednej impulznej odozvy
for Th = 1:length(h)
    h_nove_2(1,Th) = sum(h_nove(:,Th));
end

% denormaliziacia hlasitosti vytvorenej odozvy
h_max = max(max(h)); % ziskanie najvyssej hodnoty vzorky povodnej odozvy
h_nove_2_max = max(max(h_nove_2)); % ziskanie najvyssej hodnoty vzorky namodelovanej odozvy
h_nove_vysl = h_nove_2/(h_nove_2_max/h_max); % zmena najvyssej hodnoty vzorky modelovanej odozvy na najvyssiu hodnotu vzorky povodnej odozvy


% zobrazenie priebehov povodnej a namodelovanej impulznej odozvy 
figure
subplot(2,1,1) % rozdelenie grafu na 2 casti
plot(Th_vekt,h,'k') % typ casti grafu a vykreslene veliciny
title('Priebeh pôvodnej impulznej odozvy v čase'); % nazov casti grafu
xlabel('{\itt} (s) \rightarrow ')
ylabel('{\ith} (-) \rightarrow') % nazov osi
grid on % mriezka

subplot(2,1,2)
plot(Th_vekt,h_nove_vysl,'k')
title('Priebeh namodelovanej impulznej odozvy v čase');
xlabel('{\itt} (s) \rightarrow ')
ylabel('{\ith} (-) \rightarrow') % nazov osi
grid on % mriezka
set (gcf,'pos',[1370 150 550 500]) % pozicia a rozmery okna s grafom

% cesta a nazov suboru vytvorenej impulznej odozvy
Subor_nazov_n = char(strcat(Cesta,Subor_nazov(1),'-',Subor_nazov(2),'-new.wav'));

% zapis namodelovanej impulznej odozvy na disk
audiowrite(Subor_nazov_n,h_nove_vysl,fvz)