% reverberator navrhli Zafar Rafii a Bryan Pardo. Ide o modifikovany
% Moorerov reverberator, ktory je celkovo ovladany piatimi parametrami -
% oneskorenim a zosilenim prveho filtra zo sekcie paralelne zapojenych hrebenovych
% filtrov, rozdielom v oneskorenie fazovacich clankov, zosilenim riadicim
% balanc dry/wet signalu a medznym kmitoctom dolnej priepustnosti.
% 
% vstupne parametre su signal a jeho vzorkovaci kmitocet, vystup je signal
% po reverberacii a impulzna odozva reverberatoru

function [y,IR] = rafii(x,Fvz)


%% dital reverberator - Pardo, Rafii
clc
close all;

%dlzka impulznej odozvy
implength=1.2*Fvz;

%%definicia parametrov fazovacich clankov
%delay a rozdiel medzi jednotlivymi kanalmi m
da=1327;
m=163;

%delay od ktoreho su odvodene delaya comb filtrov
%zosilenie prveho comb filtra
origindelay=5147;
g=0.5;

%medzny kmitocet dolnej priepustnosti v hz
fc=3000;
% v pripade ze chceme pre jednotlive vetvy rozdielne medzne kmitocty
% fc1=5000;

%%dry/wet
wet=0.5;
dry=1-wet;
%% fir filter, prvotne odrazy
a(1)=1;

b(419)=0.68;               
b(631)=0.644;
b(757)=0.587;
b(907)=0.540;
b(1117)=0.498;
b(1327)=0.465;
b(1549)=0.413;
b(1871)=0.371;
b(2161)=0.324;
b(2459)=0.281;
b(2659)=0.217;
b(2857)=0.181;
b(3037)=0.180;

early=impz(b,a);
%% definicia prvkov neskorej reverberacie

%%comb1
k=round(origindelay*1.5);

T60=k*((log10(0.001))/log10(g));

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

b(1)=0;
b(k+1)=1;
a(1)=1;
a(k+1)=-g;

comb1=tf(b,a);

%%comb2
k=round(origindelay*1.42);
g=10^((k*log10(0.001))/T60);

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

b(1)=0;
b(k+1)=1;
a(1)=1;
a(k+1)=-g;

comb2=tf(b,a);

%%comb3
k=round(origindelay*1.34);
g=10^((k*log10(0.001))/T60);

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

b(1)=0;
b(k+1)=1;
a(1)=1;
a(k+1)=-g;

comb3=tf(b,a);

%%comb4
k=round(origindelay*1.26);
g=10^((k*log10(0.001))/T60);

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

b(1)=0;
b(k+1)=1;
a(1)=1;
a(k+1)=-g;

comb4=tf(b,a);

%%comb5
k=round(origindelay*1.16);
g=10^((k*log10(0.001))/T60);

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

b(1)=0;
b(k+1)=1;
a(1)=1;
a(k+1)=-g;

comb5=tf(b,a);

%%comb6
k=round(origindelay*1.07);
g=10^((k*log10(0.001))/T60);

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

b(1)=0;
b(k+1)=1;
a(1)=1;
a(k+1)=-g;

comb6=tf(b,a);

%%alpassleft
g=1/sqrt(2);             
k=da+round(m/2);              
      
b=zeros(1,k+1);
a=zeros(1,k+1);

b(1)=g;
b(k+1)=1;
a(1)=1;
a(k+1)=g;

APleft=tf(b,a);

%%alpassright
g=1/sqrt(2);            
k=da-round(m/2);              
      
b=zeros(1,k+1);
a=zeros(1,k+1);

b(1)=g;
b(k+1)=1;
a(1)=1;
a(k+1)=g;

APright=tf(b,a);

%%Lowpassleft
fc=fc/Fvz;

g=2-cos(2*pi*fc)-sqrt(((cos(2*pi*fc)-2)^2)-1);             

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

LPleft=tf(b,a);

%%Lowpassright

% v pripade ze chceme pre jednotlive vetvy rozdielne medzne kmitocty
% fc=fc1/Fvz;   

g=2-cos(2*pi*fc)-sqrt(((cos(2*pi*fc)-2)^2)-1);            

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

LPright=tf(b,a);

%% zjednodušenie štruktúry

% paralelna kombinacia hrebenovych filtrov
combparallel=comb1-comb2+comb3-comb4+comb5-comb6;
disp('paralelna kombinacia hrebenovych filtrov hotova');

% seriova kombinacia v oboch kanaloch
allpassleft=APleft*LPleft*wet;
allpassright=APright*LPright*wet;
disp('seriova kombinacia alpass filtrov hotova');

% seriove zapojeni zjednodusenych hrebenovych filtrov a lavej, respektive pravej vetvy s
% fazovacim clanok a dolnou priepustnostou
wholeleft=combparallel*allpassleft;
disp('celkova prenosova funcia laveho kanala');
[IRb,IRa]=tfdata(wholeleft);
IRb=IRb{1,1};
IRa=IRa{1,1};

wholeright=combparallel*allpassright;
disp('celkova prenosova funcia praveho kanala');
[IRbb,IRaa]=tfdata(wholeright);
IRbb=IRbb{1,1};
IRaa=IRaa{1,1};

%celkove impulzne odozvy laveho a praveho kanalu
IRl=impz(IRb,IRa,implength);
IRr=impz(IRbb,IRaa,implength);

%konvolucia s prvotnych odrazov s neskorou reverberaciou
[early,IRl]=eqtflength(early,IRl);
IRleft=conv(early,IRl);

[early,IRr]=eqtflength(early,IRr);
IRright=conv(early,IRr);

[IRleft,IRright]=eqtflength(IRleft,IRright);

% spojenie laveho a praveho kanalu do stereo impulznej odozvy
IR(:,1)=IRleft;
IR(:,2)=IRright;
%% filtrovanie

% ak je signal stereo upravi sa na mono
if size(x,2)~=1 
xcombin=0.5*x(:,1)+0.5*x(:,2);
elseif size(x,2)==1
    xcombin=x;
end

%ziskanie laveho a praveho kanalu sterea
yl=conv(xcombin,IRleft);
yr=conv(xcombin,IRright);

% kombinacia dry a wet signalu
if dry==1
    y=x;
elseif size(x,2)==2
y(:,1)=vektorplus(yl,(dry*x(:,1))');
y(:,2)=vektorplus(yr,(dry*x(:,2))');
elseif size(x,2)==1
y(:,1)=vektorplus(yl',dry*x'); 
y(:,2)=vektorplus(yr',dry*x');
end

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

%odstrani nuly z konca impulznej odozvy
I=find(abs(IR(:,1))>10^(-5),1,'last');
IR=IR(1:I,:);

%% zobrazenie impulznej odozvy praveho kanalu
figure('Position',[10 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 reverberátoru');
set(h,'FontName','MS Sans Serif')
xlabel('{{\itt}} (s) \rightarrow')
ylabel('\ith(-) \rightarrow')

grid on
%% zobrazenie vstupneho signalu
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
%% zobrazenie vystupneho signalu
n=(1:length(y));

subplot(3,1,3)
plot(n/Fvz,y);
h=title('\bfVýstupny 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(IRleft);
I=fft(IRleft)/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 signalu do suboru
% wavwrite(y,Fvz,16,'reverberant.wav');





