%tento reverberator bol navrhnuty Vanaanenovou a kolegami. Vstup je stereo alebo mono
%signal a vzorkovaci kmitocet.Vystup je stereo signal a celkova impulzna
%odozva reverberatoru

function [y,IRtotal]=vanaanen(x,Fvz)
close all;
clc;

Fvz=44100;          %vzorkovaci kmitocet
implength=1.7*Fvz;   %dlzka vektoru impulznej odozvy vo vzorkoch
beta=0.2;
wet=0.5;

%oneskorenia v jednotlivych linkach
k1=2957;
k2=4019;
k3=4391;
k4=4463;

%parametre jednotlivých fazovacich clankov
ga1=0.5;            % zosilenie 
k11=619;              % oneskorenie
ga2=0.5;            % zosilenie 
k22=853;              % oneskorenie  
ga3=0.5;            % zosilenie 
k33=983;              % oneskorenie
ga4=0.5;            % zosilenie 
k44=1009;              % oneskorenie

%celkove oneskorenie jednotlivych liniek, fazovacie clanky+oneskorenia
delay1=k1+k11;
delay2=k2+k22;
delay3=k3+k33;
delay4=k4+k44;
%% prvotne odrazy

%lavy kanal
bleft(709)=0.841;               
bleft(859)=0.504;
bleft(1031)=0.490;
bleft(1361)=0.379;
bleft(1699)=0.307;
bleft(1031)=0.490;
bleft(1361)=0.379;
bleft(1699)=0.307;
%pravy kanal
bright(811)=0.811;
bright(919)=0.484;
bright(1049)=0.453;
bright(1423)=0.351;
bright(1753)=0.282;
bright(1049)=0.453;
bright(1423)=0.351;
bright(1753)=0.282;

%impulzna odozva
FIR_left=impz(bleft,1);
FIR_right=impz(bright,1);
%% definicia prvkov modulu neskorej reverberacie

%%delay1
b=zeros(1,k1+1);
a=zeros(1,k1+1);

b(k1+1)=1;
a(1)=1;

D1=tf(b,a);

%%delay2
b=zeros(1,k2+1);
a=zeros(1,k2+1);

b(k2+1)=1;
a(1)=1;

D2=tf(b,a);

%%delay3
b=zeros(1,k3+1);
a=zeros(1,k3+1);

b(k3+1)=1;
a(1)=1;

D3=tf(b,a);

%%delay4
b=zeros(1,k4+1);
a=zeros(1,k4+1);

b(k4+1)=1;
a(1)=1;

D4=tf(b,a);

%%lowpass filtre

%%lowpass1
k=0.34;
g=0.74;                           
b=[0 (g-(g*k))];
a=[1 -k];

LP1=tf(b,a);

%%lowpass2
k=0.45;
g=0.66;                           
b=[0 (g-(g*k))];
a=[1 -k];


LP2=tf(b,a);

%%lowpass3
k=0.48;
g=0.64;                           
b=[0 (g-(g*k))];
a=[1 -k];

LP3=tf(b,a);

%%lowpass4
k=0.49;
g=0.63;                           
b=[0 (g-(g*k))];
a=[1 -k];

LP4=tf(b,a);

%%allpass filtre

%%alpass1
   
b=zeros(1,k11+1);
a=zeros(1,k11+1);
b(1)=ga1;
b(k11+1)=1;
a(1)=1;
a(k11+1)=ga1;

AP1=tf(b,a);

%%alpass2
    
b=zeros(1,k22+1);
a=zeros(1,k22+1);
b(1)=ga2;
b(k22+1)=1;
a(1)=1;
a(k22+1)=ga2;

AP2=tf(b,a);

%%alpass3
      
b=zeros(1,k33+1);
a=zeros(1,k33+1);
b(1)=ga3;
b(k33+1)=1;
a(1)=1;
a(k33+1)=ga3;

AP3=tf(b,a);

%%alpass4
      
b=zeros(1,k44+1);
a=zeros(1,k44+1);
b(1)=ga4;
b(k44+1)=1;
a(1)=1;
a(k44+1)=ga4;

AP4=tf(b,a);
%% seriove zapojenie prvkov

%%vetva1
vetva1direct=D1*LP1*AP1;
vetva1=feedback(vetva1direct,1,+1);

[b1,a1]=tfdata(vetva1,'v');
IR1=impz(b1,a1,implength);

[FIR_left,IR1]=eqtflength(FIR_left,IR1);
IR1=conv(FIR_left,IR1);
IR1=IR1(1:implength);

%%vetva2
vetva2direct=D2*LP2*AP2;
vetva2=feedback(vetva2direct,1,+1);

[b2,a2]=tfdata(vetva2,'v');
IR2=impz(b2,a2,implength);

[FIR_left,IR2]=eqtflength(FIR_left,IR2);
IR2=conv(FIR_left,IR2);
IR2=IR2(1:implength);

%%vetva3
vetva3direct=D3*LP3*AP3;
vetva3=feedback(vetva3direct,1,+1);

[b3,a3]=tfdata(vetva3,'v');
IR3=impz(b3,a3,implength);

[FIR_left,IR3]=eqtflength(FIR_left,IR3);
IR3=conv(FIR_left,IR3);
IR3=IR3(1:implength);

%%vetva4
vetva4direct=D4*LP4*AP4;
vetva4=feedback(vetva4direct,1,+1);

[b4,a4]=tfdata(vetva4,'v');
IR4=impz(b4,a4,implength);

[FIR_left,IR4]=eqtflength(FIR_left,IR4);
IR4=conv(FIR_left,IR4);
IR4=IR4(1:implength);

IR=IR1+IR2+IR3+IR4;

IR1buffer(implength+delay1)=0;
IR2buffer(implength+delay2)=0;
IR3buffer(implength+delay3)=0;
IR4buffer(implength+delay4)=0;  

%spetna väzba
for i=1:2
IR1buffer=delay(IR,delay1)*beta;
IR2buffer=delay(IR,delay2)*beta;
IR3buffer=delay(IR,delay3)*beta;
IR4buffer=delay(IR,delay4)*beta;

IR1=vektorplus(IR1,IR1buffer(1:implength));
IR2=vektorplus(IR2,IR2buffer(1:implength));
IR3=vektorplus(IR3,IR3buffer(1:implength));
IR4=vektorplus(IR4,IR4buffer(1:implength));

IR=IR1+IR2+IR3+IR4;
end

IRtotal=-IR1-IR2+IR3+IR4;

[FIR_left,IRtotal]=eqtflength(FIR_left,IRtotal);
IRleft=conv(FIR_left,IRtotal);

[FIR_right,IRtotal]=eqtflength(FIR_right,IRtotal);
IRright=conv(FIR_right,IRtotal);

IRfinal(:,1)=IRleft';
IRfinal(:,2)=IRright';
%% filtrovanie

if size(x,2)==2
xl=x(:,1);
xr=x(:,2);

[xl,IRleft]=eqtflength(xl,IRleft);
yl=conv(xl,IRleft);
yr=conv(xr,IRright);
size(yl)
size(xl)
size(yr)
size(xr)
pause
yleft=vektorplus((wet*yl),((1-wet)*xl));
yright=vektorplus((wet*yr),((1-wet)*xr'));
[yleft,yright]=eqtflength(yleft,yright);

y(:,1)=yleft;
y(:,2)=yright;
elseif size(x,2)==1
[x,IRleft]=eqtflength(x,IRleft);
yl=conv(x,IRleft);
[x,IRright]=eqtflength(x,IRright);
yr=conv(x,IRright);

yleft=vektorplus(wet*yl,(1-wet)*x);
yright=vektorplus(wet*yr,(1-wet)*x);
[yleft,yright]=eqtflength(yleft,yright);

y(:,1)=yleft;
y(:,2)=yright;
% y=y';
end

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

%% IMPULZNA ODOZVA

figure('Position',[0 20 675 800])  % "Position" nastavi velkost a poziciu okna
                                  % figure('Position',[left bottom width heigth])  
set(gcf, 'color', [1 1 1]) % nastavenie farby pozadia
n=1:length(IRfinal);

subplot(3,1,1)
plot(n/Fvz,IRfinal);
h=title('\bfReverberátor');
set(h,'FontName','MS Sans Serif')
xlabel('{{\itt}} (s) \rightarrow')
ylabel('\ith(-) \rightarrow')
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('\ith(-) \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('\ith(-) \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(IRfinal);
I=fft(IRfinal)/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 suboru 
% % wavwrite(y,Fvz,16,'reverberant.wav');
