%% FDN reverberator postaveny na zaklade jotovho reverberatoru
% ide o reverberator postaveny na oneskorovacej spetnovezobnej sieti prvykrat
% navrhnuty Jotom, v ktorom su prevedene drobne modifikacie - zvecseny
% pocetoneskorovacich liniek, tomu odpovedajuca ina matica, pridany modul
% prvotnych odrazov atd. Ide o funkciu ktorej vstupom je stereo alebo mono
% signal a vystup je vystupny mono signal a impulzna odozva
% reverberatoru.

function [y,IRfinal]=fdnjot(x,Fvz)

close all;
clc;

%dry wet a dlzka impulznej odozvy
wet=0.5;
Td=0.7;
implength=Td*Fvz;   %dlzka vektoru impulznej odozvy vo vzorkoch

%% prvotne odrazy FIR filter
earlya=1;

earlyb(419)=0.511;  
earlyb(759)=-0.422;  
earlyb(1307)=0.352;
earlyb(2017)=-0.312; 
earlyb(2143)=0.372;
earlyb(2539)=-0.23;
earlyb(3319)=0.198; 
earlyb(3571)=-0.159;
earlyb(3797)=0.118;
earlyb(4373)=0.078;

FIR=impz(earlyb,earlya);
%% definica oneskorenia
oneskorenie=[1579 1949 2293 2647 2903 3217 3647 4111]; %vektor oneskoreni v jednotlivych linkach vektoru 
                                                       % vo vzorkoch
                                                       % -prvocisla
impzasob(1,oneskorenie(1))=1;   
impzasob(2,oneskorenie(2))=1;
impzasob(3,oneskorenie(3))=1;
impzasob(4,oneskorenie(4))=1;
impzasob(5,oneskorenie(5))=1;   
impzasob(6,oneskorenie(6))=1;
impzasob(7,oneskorenie(7))=1;
impzasob(8,oneskorenie(8))=1;
%% vektor zosileni g
g=10.^(-3*oneskorenie/Td/Fvz);    %definicia zosilenia, vztah [4] v kapitole 3
g=g';                       
%% matica zosileni,ortogonalna
a=-0.25;
A=[0.75 a a a a a a a;
   a a 0.75 a a a a a
   a 0.75 a a a a a a
   a a a a 0.75 a a a
   a a a 0.75 a a a a
   a a a a a a 0.75 a
   a a a a a 0.75 a a
   a a a a a a a 0.75]; 
%% spatna vazba, generovanie impulznej odozvy

for i=1:implength
    
    impzasob(:,i)=impzasob(:,i).*g;             %na patricnu poziciu kde je zpozdenie, vynasobi prislusnym oneskorenieom
                                        
    imp(:,i)=impzasob(:,i);                     % vystup
     
    feedback=impzasob(:,i)'*A;                  % spatna vazba
    impzasob(1,i+oneskorenie(1))=feedback(1);   % nacitanie do impzasobu,mnozenie ozvien v impzasobi, modelovanie echa
    impzasob(2,i+oneskorenie(2))=feedback(2);
    impzasob(3,i+oneskorenie(3))=feedback(3);
    impzasob(4,i+oneskorenie(4))=feedback(4);
    impzasob(5,i+oneskorenie(5))=feedback(5);  
    impzasob(6,i+oneskorenie(6))=feedback(6);
    impzasob(7,i+oneskorenie(7))=feedback(7);
    impzasob(8,i+oneskorenie(8))=feedback(8);
    
end 

% sucet impulznych odoziev z jednotlivych liniek
IR=sum(imp);
%% filtrovanie lowpassom

%%definicia dolnej priepustnosti
fc=3000/Fvz;

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

bb=[1-gg 0];
aa=[1 -gg];

iLP=impz(bb,aa,implength); % impulzna odozva low pass filtru

IRfinal=conv(IR,iLP);      % filtrovanie impulznej odozvy lowpassom

%seriove zapojenie FIR filtru
[FIR,IRfinal]=eqtflength(FIR,IRfinal);
IRfinal=conv(FIR,IRfinal);

%paralelne zapojenie FIR filtru, odkomentovat ak ho chce uzivatel pouzit
% a zakomentovat seriove
% q=delay(IRfinal,1571);     % oneskorenie difuzneho dozvuku od prvotnych odrazov
% IRfinal=vektorplus(FIR,q);
% IRfinal=IRfinal(1:implength);
%% filtrovanie signalu
if size(x,2)==1
[x,IRfinal]=eqtflength(x,IRfinal);
y=conv(xmono,IRfinal);
elseif size(x,2)~=1
xl=x(:,1);
xr=x(:,2);
[xl,IRfinal]=eqtflength(xl,IRfinal);
y(:,1)=conv(xl,IRfinal);
[xr,IRfinal]=eqtflength(xr,IRfinal);
y(:,2)=conv(xr,IRfinal);
end

% dry/wet
y=vektorplus(wet*y,(1-wet)*x);
y=y';

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

%odstrani nuly
I=find(abs(IRfinal)>10^(-5),1,'last');
IRfinal=IRfinal(:,1:I);
%% vykreslenie impluznej odozvy
n=(1:length(IRfinal));
figure('Position',[10 20 755 900])  % "Position" nastavi velkost a poziciu okna
                                 
subplot(3,1,1)
plot(n/Fvz,IRfinal);
h=title('\bfImpulzná charakteristika reverberátoru');
set(h,'FontName','MS Sans Serif')
xlabel('{{\itt}} (s) \rightarrow')
ylabel('\ith(-) \rightarrow')
% set(gca,'ylim',[-0.15 0.5])
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ý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
                                   
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 signalu do suboru
% wavwrite(y,Fvz,16,'paralel.wav');

