function [ y ] = IIRzerophase( x, Bcasc, Acasc )
%[y] = IIRZEROPHASE( x, Bcasc, Acasc )
%IIRZEROPHASE vykoná filtráciu metódou forward-backward, pri ktorej nedochádza
%k časovému posunu výstupného signálu od vstupného. Filter je
%reprezentovaný čitateľom a menovateľom prenosu
% Vstupy:   x - vstupný signál filtru vo forme matice
%           Bcasc - jednotlivé sekcie čitateľa prenosu
%           Acasc - jednotlivé sekcie menovateľa prenosu
% Výstup:   y - výstupný signál filtru
%inšpirované matlabovskou funkciou filtfilt.m
%[1] https://www.dsprelated.com/freebooks/filters/State_Space_Realization.html
%[2] https://www.dsprelated.com/freebooks/filters/State_Space_Filters.html
%[3] https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/12886/versions/3/previews/Teqcspec/filtfilt.m/index.html

assert(isnumeric(x)==1,'Vstupný signál musí pozostávať z číselných hodnôt');
assert(isnumeric(Bcasc)==1,'Čitateľ prenosu musí pozostávať z číselných hodnôt');

b = size(Bcasc,1);

if nargin<3
    Acasc = [zeros(b,1) ones(b,1)];
end
assert(isnumeric(Acasc)==1,'Menovateľ prenosu musí pozostávať z číselných hodnôt');

a = size(Acasc,1);

assert(a==b,'Počet sekcií čitateľa a menovateľa musí byť rovnaký');

ord = max(numel(Acasc),numel(Bcasc))-b;
[K,N] = size(x);

Nin = 3*ord;

if N < Nin
    assert(1==0,'Vstupný signál je krátky');
end

n = max(size(Acasc,2)-1,size(Bcasc,2)-1);

%% URCENIE POCIATOCNEHO STAVU FILTRU [3]

zin = zeros(b,n);
for i=1:b
    Bh = IIRtruncate(Bcasc(i,:));
    Ah = IIRtruncate(Acasc(i,:));
    Bh(n+2) = 0;
    Ah(n+2) = 0;
    ps  = (Bh(2:n+1) - Bh(1)*Ah(2:n+1));
    res = ( eye(n) - [-Ah(2:n+1);eye(n-1,n)]' ) \ ps';
    zin(i,:) = res';
end

%% FILTRACIA S NULOVOU FAZOU

yn = zeros(K,N);

for j = 1:K
    y = x(j,:);    
    
    for i=1:b
        y = [2*y(1)-y(Nin+1:-1:2) y 2*y(end)-y(end-1:-1:end-Nin)];  %nutné rozšírenie, aby koncové hodnoty neutekali od vstupného signálu

        y = IIRfilt(y,Bcasc(i,:),Acasc(i,:),zin(i,:)*y(1));         %filtrácia prvým smerom
        y = y(end:-1:1);                                            %otočenie signálu
        y = IIRfilt(y,Bcasc(i,:),Acasc(i,:),zin(i,:)*y(1));         %filtrácia druhým smerom pre vynulovanie fázy

        y = y(end-Nin:-1:Nin+1);                                    %zúženie signálu na počet vzorkov vstupného
    end  
    yn(j,:) = y;
end

y = yn;

end

