function [B,A,Bcasc,Acasc] = butterworth( Fs, Fpass, Fstop, Gp, Gs, typ, orderconstraint )
%[B,A,Bcasc,Acasc] = BUTTERWORTH( Fs, Fpass, Fstop, Gp, Gs, typ, orderconstraint )
%BUTTERWORTH z daných parametrov amplitúdovej odozvy vytvorí čitateľ a menovateľ prenosu frekvenčne priepustného
%           filtru s možnosťou LOWPASS a HIGHPASS, pre BANDPASS a BANDSTOP
%           viď BUTTERWORTHBAND
%Vstupy:    Fs - vzorkovacia frekvencia v Hertzoch
%           Fpass - hranica priepustnej oblasti (passband) - frekvencia v
%           Herzoch
%           Fstop - hranica blokujúcej oblasti (stopband) - frekvencia v
%           Hertzoch
%           Gp - amplitúdový prenos na hranici priepustnej oblasti (záporné decibely)
%           Gs - amplitúdový prenos na hranici blokujúcej oblasti (záporné decibely)
%           typ - 'low'/'high' na rozlíšenie návrhu lowpassu a highpassu
%           orderconstraint - v prípade požiadavku na párny stupeň filtru,
%           treba zadať 'even'
%Výstupy    A - vektor koeficientov menovateľa prenosu
%           B - vektor koeficientov čitateľa prenosu
%           Acasc - jednotlivé sekcie menovateľa prenosu
%           Bcasc - jednotlivé sekcie čitateľa prenosu
%[1] http://www.ee.ic.ac.uk/pcheung/teaching/ee2_signals/Lecture%209%20-%20Poles%20Zeros%20&%20Filters.pdf
%[2] https://www.youtube.com/watch?v=-EHRFrDujhc&t=5s
%[3] https://www.youtube.com/watch?v=hg3wMmTI_Oc
%[4] https://spinlab.wpi.edu/courses/ece503_2014/10-3bilinear_transform.pdf
%[5] PARKS, T. W. a C. S. BURRUS. Digital filter design. New York: Wiley, c1987. ISBN 9780471828969. 

% IIR frekvenčne priepustný filter - butterworth návrh [5]

if nargin<7
    orderconstraint = 'none';
end

Ts = 1;
assert(Gp>Gs, 'Navrhované decibely nesplňujú predpoklady klesajúcej amplitúdovej odozvy');

switch typ
    
    case 'low'
        fpass = Fpass/Fs;
        fstop = Fstop/Fs;
        
    case 'high'
        fpass = 0.5 - Fpass/Fs;
        fstop = 0.5 - Fstop/Fs;
        
    otherwise
        assert(1==0,'Neplatná voľba typu filtru');
end

if fpass > fstop
   hlp = fstop;
   fstop = fpass;
   fpass = hlp;
end

omp = 2*pi*fpass;
oms = 2*pi*fstop;

%% NAVRH ANALOGOVEHO EKVIVALENTU [1],[2],[3]

OMP = 2/Ts * tan(omp/2);
OMS = 2/Ts * tan(oms/2);

n = log10((10^(-Gs/10)-1)/((10^(-Gp/10)-1)))/(2*log10(OMS/OMP));
n = ceil(n);
if strcmp(orderconstraint,'even')==1
   if mod(n,2)==1
       n=n+1;
   end
end
OMC = OMP/((10^(-Gp/10)-1)^(1/2/n));

pk = zeros(1,floor(n/2));
for k = 1:floor(n/2)
    pk(k) = OMC*exp((1i*pi/2/n)*(2*k+n-1));
end

%% BILINEARNA TRANSFORMACIA NA DIGITALNY FILTER [4]

Pk = (1+pk*(Ts/2))./(1-pk*(Ts/2));
Zk = -ones(1,floor(n/2));
[A,B,Acasc,Bcasc] = gettf(Pk,Zk);

for i = 1:size(Acasc,1) %normalizácia priepustnej oblasti (passband)
    K = sum(Bcasc(i,:))/sum(Acasc(i,:));
    Bcasc(i,:) = Bcasc(i,:)./K;
    B = B./K;
end

if mod(n,2)==1 %v prípade nepárneho stupňa filtru leží jeden pól na reálnej ose
   ppos = -OMC;
   Ppos = (1+ppos*(Ts/2))./(1-ppos*(Ts/2));
   Acasc(end+1,:) = [0 1 -Ppos];
   Bcasc(end+1,:) = [0 1 1];
   K = sum(Bcasc(end,:))/sum(Acasc(end,:));
   Bcasc(end,:) = Bcasc(end,:)./K;
   A = convolution(A,[1 -Ppos]);
   B = convolution(B,[1 1]);
   B = B./K;
end

if strcmp(typ,'high') == 1
    [A,B] = spectralrev(A,B);
    for i = 1:size(Acasc,1)
       Acasc(i,:) = spectralrev(Acasc(i,:));
       Bcasc(i,:) = spectralrev(Bcasc(i,:));
    end
end


end

