% Implementacia Dattorovho reverberatoru. Tento reverberator vytvara iba
% difuzny dozvuk. Vstup funkcie je stereo alebo mono signal a vzorkovaci
% kmitocet. Vystup je vstupny signal po reverberacii a Impulzna odozva
% reverberatoru.Na vytvorenie impulznej odozvy reverberatoru je pouzity
% jednotkovy impulz. 

function [y,IR] = dattorro(x,Fvz)
clc
close all

%jendotkovy impulz na vytvorenie impulznej odozvy
implength=1.5*Fvz;
xmono=1;
xmono(implength)=0;

%% definicia parametrov

%parametre reverberárotoru
decay = 0.5;
decay_diffusion_1 = 0.7;
decay_diffusion_2 = 0.5;
input_diffusion_1 = 0.75;
input_diffusion_2 = 0.625;
bandwidth = 0.0005;
damping=0.90;
wet=0.5;

%oneskorenia jednotlivych prvkov
delay1 = 210;
delay2 = 159;
delay3 = 561;
delay4 = 410;
delay5 = 996;
delay6 = 6598;
delay7 = 2667;
delay8 = 6249;
delay55 = 1345;
delay66 = 6249;
delay77 = 3935;
delay88 = 4687;

%oneskorenia v output tapoch tanku pre lavy a pravy kanal
leftdelay1 = 394;
leftdelay2 = 4407;
leftdelay3 = 2834;
leftdelay4 = 2957;
leftdelay5 = 2949;
leftdelay6 = 277;
leftdelay7 = 1580;

rightdelay1 = 523;
rightdelay2 = 5375;
rightdelay3 = 1820;
rightdelay4 = 3961;
rightdelay5 = 3128;
rightdelay6 = 496;
rightdelay7 = 179;

%% definicia prvkov

%lowpass

b=[1-bandwidth 0];
a=[1 -bandwidth];

LP1=tf(b,a);

%allpass1
k=delay1;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=input_diffusion_1;
b(k+1)=1;
a(1)=1;
a(k+1)=input_diffusion_1;

AP1=tf(b,a);  

%allpass2
k=delay2;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=input_diffusion_1;
b(k+1)=1;
a(1)=1;
a(k+1)=input_diffusion_1;

AP2=tf(b,a);

%allpass3
k=delay3;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=input_diffusion_2;
b(k+1)=1;
a(1)=1;
a(k+1)=input_diffusion_2;

AP3=tf(b,a);

%allpass4
k=delay4;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=input_diffusion_2;
b(k+1)=1;
a(1)=1;
a(k+1)=input_diffusion_2;

AP4=tf(b,a);

%tank vetva 1

%allpass5
k=delay5;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=decay_diffusion_1;
b(k+1)=1;
a(1)=1;
a(k+1)=decay_diffusion_1;

AP5=tf(b,a);

%delay
k=delay6;

b=zeros(1,k+1);
b(k+1)=1;
a=1;

D1=tf(b,a);

%lowpass

b=[1-damping 0];
a=[1 -damping];

LP2=tf(b,a);

%allpass6
k=delay7;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=decay_diffusion_2;
b(k+1)=1;
a(1)=1;
a(k+1)=decay_diffusion_2;

AP6=tf(b,a);

%delay
k=delay8;

b=zeros(1,k+1);
b(k+1)=1;
a=1;

D2=tf(b,a);

%tank vetva 2

%allpass55
k=delay55;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=decay_diffusion_1;
b(k+1)=1;
a(1)=1;
a(k+1)=decay_diffusion_1;

AP55=tf(b,a);

%delay
k=delay66;

b=zeros(1,k+1);
b(k+1)=1;
a=1;

D11=tf(b,a);

%lowpass

b=[1-damping 0];
a=[1 -damping];

LP22=tf(b,a);

%allpass66
k=delay77;

b=zeros(1,k+1);
a=zeros(1,k+1);
b(1)=decay_diffusion_2;
b(k+1)=1;
a(1)=1;
a(k+1)=decay_diffusion_2;

AP66=tf(b,a);

%delay
k=delay88;

b=zeros(1,k+1);
b(k+1)=1;
a=1;

D22=tf(b,a);

%% prenosove funcie kombinacii prvkov

%dekorelacna faza
first=LP1*AP1*AP2*AP3*AP4;
[b1,a1]=tfdata(first);
b1=b1{1,1};
a1=a1{1,1};

%prenosova funkcia tanku
tank1=AP5*D1*LP2*decay*AP6*D2*decay;
tank2=AP55*D11*LP22*decay*AP66*D22*decay;

tank=(1+tank2)/(1-tank1*tank2);
[b2,a2]=tfdata(tank);
b2=b2{1,1};
a2=a2{1,1};

%% reverberacia

%dekorelacia signalu a ziskanie referencneho signalu v tanku 
beforetank=filter(b1,a1,xmono);
tank_ref=filter(b2,a2,beforetank);

%ziskavanie oneskorenych vystupnych signalov z roznych bodov tanku
[node1] = allpass(tank_ref, delay5,-decay_diffusion_1);
[leftchannel] = delay(-node1, leftdelay5);
[rightchannel] = delay(node1, rightdelay1);
[rightbuffer] = delay(node1, rightdelay2);
rightchannel = vektorplus(rightchannel, rightbuffer);

[node2] = delay(node1, delay6);

[node3] = lowpass(decay*node2, damping);
[leftbuffer] = allpasszasobnik(node3, delay7,decay_diffusion_2, leftdelay6);
leftchannel = vektorplus(leftchannel, leftbuffer);
[rightbuffer] = allpasszasobnik(node3, delay7, decay_diffusion_2, rightdelay3);
rightchannel = vektorplus(rightchannel, rightbuffer);

[node4] = allpass(node3, delay7, decay_diffusion_2);
[leftbuffer] = delay(-node4, leftdelay7);
leftchannel = vektorplus(leftchannel, leftbuffer);
[rightbuffer] = delay(node4, rightdelay4);
rightchannel = vektorplus(rightchannel, rightbuffer);

[node5] = delay(node4, delay8);

[node6] = allpass((beforetank+decay*node5), delay55, -decay_diffusion_1);

[leftbuffer] = delay(node6, leftdelay1);
leftchannel = vektorplus(leftchannel, leftbuffer);
[leftbuffer] = delay(node6, leftdelay2);
leftchannel = vektorplus(leftchannel, leftbuffer);
[rightbuffer] = delay(-node6, rightdelay5);
rightchannel = vektorplus(rightchannel, rightbuffer);

[node7] = delay(node6, delay66);

[node8] = lowpass(decay*node7, damping);
[leftbuffer] =  allpasszasobnik(node8, delay77, decay_diffusion_2, leftdelay3);
leftchannel = vektorplus(leftchannel, leftbuffer);
[rightbuffer] = allpasszasobnik(node8, delay77, decay_diffusion_2, rightdelay6);
rightchannel = vektorplus(rightchannel, rightbuffer);

[node9] = allpass(node8, delay88, decay_diffusion_2);
[leftbuffer] = delay(node9, leftdelay4);
leftchannel = vektorplus(leftchannel, leftbuffer);
[rightbuffer] = delay(-node9, rightdelay7);
rightchannel = vektorplus(rightchannel, rightbuffer);

%vytvorenie stereo signalu
IR(:,1)=leftchannel;
IR(:,2)=rightchannel;

% filtrovanie vstupneho signalu
if size(x,2)==2
    xl=x(:,1);
    xr=x(:,2);
    
    [xl,leftchannel]=eqtflength(xl,leftchannel);
    [xr,rightchannel]=eqtflength(xr,rightchannel);
    
    y(:,1)=conv(xl, leftchannel);
    y(:,2)=conv(xr, rightchannel);
    
    y(:,1)=vektorplus(wet*y(:,1),((1-wet)*xl)');
    y(:,2)=vektorplus(wet*y(:,2),((1-wet)*xr)');
    
elseif size(x,2)==1
     [x,leftchannel]=eqtflength(x,leftchannel);
     y=conv(x,leftchannel);
     
     %dry wet
     y=vektorplus(wet*y,(1-wet)*x);
     y=y';    
     
end

%%odstranenie nul z konca signalu
I=find(abs(y(:,1))>10^(-5),1,'last');
y=y(1:I,:);

%%odstranenie nul z konca signalu
I=find(abs(IR(:,1))>10^(-5),1,'last');
IR=IR(1:I,:);
%% impulzna odozva
figure('Position',[20 20 755 900])  % "Position" nastavi velkost a poziciu okna
                             
n=(1:length(IR));

subplot(3,1,1)
plot(n/Fvz,IR);
h=title('\bfImpulzná odozva');
set(h,'FontName','MS Sans Serif')
xlabel('{{\itt}} (s) \rightarrow')
ylabel('\its(-) \rightarrow')
set(gca,'ylim',[-0.2 0.2])
grid on

%% vstupny signal

n=(1:length(x));

subplot(3,1,2)
plot(n/Fvz,x);
h=title('\bfVstupný signál');
set(h,'FontName','MS Sans Serif')
xlabel('{{\itt}} (s) \rightarrow')
ylabel('\its(-) \rightarrow')
grid on

%% vystupny signal
n=(1:length(y));

subplot(3,1,3)
plot(n/Fvz,y);
h=title('\bfVýstupný signál');
set(h,'FontName','MS Sans Serif')
xlabel('{{\itt}} (s) \rightarrow')
ylabel('\its(-) \rightarrow')
grid on

%% spektra vstupneho a vystupneho signalu + reverberatoru

figure('Position',[800 20 755 900])  % "Position" nastavi velkost a poziciu okna
                                   
set(gcf, 'color', [1 1 1]) % nastavenie farby pozadia

Nd=length(IR);
I=fft(IR)/Nd;                                        %výpočet DFT signálu y o délce Nd pomocí FFT, výsledek dělen Nd
I=I(1:(length(I)/2));
n=(1:length(I));

%odkomentovat v pripade ze uzivatel chce aj kmitoctove charakteristiky
%vstupneho a vystupneho signalu, pripadne zapisat vystupny signal do suboru
% subplot(3,1,1)

stem(n/Nd*Fvz,abs(I),'marker','none');
h=title('\bfModulová kmitočtová charakteristika reverberátoru');
set(h,'FontName','MS Sans Serif')
xlabel('{{\itf}} (Hz) \rightarrow')
ylabel('Modul (-) \rightarrow')
set(gca,'xlim',[20 20000])
grid on

%odkomentovat v pripade ze uzivatel chce aj kmitoctove charakteristiky
%vstupneho a vystupneho signalu, pripadne zapisat vystupny signal do suboru

% Nd=length(x);
% X=fft(x)/Nd; %výpočet DFT signálu y o délce Nd pomocí FFT, výsledek dělen Nd
% X=X(1:(length(X)/2));
% n=(1:length(X));
% subplot(3,1,2)
% stem(n/Nd*Fvz,abs(X),'marker','none');
% h=title('\bfSpektrum vstupného signálu');
% set(h,'FontName','MS Sans Serif')
% xlabel('{{\itf}} (Hz) \rightarrow')
% ylabel('Modul (-) \rightarrow')
% set(gca,'xlim',[20 20000])
% grid on
% 
% Nd=length(y);
% Y=fft(y)/Nd;                                        %výpočet DFT signálu y o délce Nd pomocí FFT, výsledek dělen Nd
% Y=Y(1:(length(Y)/2));
% n=(1:length(Y));
% subplot(3,1,3)
% stem(n/Nd*Fvz,abs(Y),'marker','none');
% h=title('\bfspektrum výstupného signálu');
% set(h,'FontName','MS Sans Serif')
% xlabel('{{\itf}} (Hz) \rightarrow')
% ylabel('Modul (-) \rightarrow')
% set(gca,'xlim',[20 20000])
% grid on


% zapis vystupneho signalu do suboru
% wavwrite(y,Fvz,16,'reverberant.wav');
