function [Bd,Ad,n] = iirdig(fp, fs, Ap, As, fm, type)
%
% function [B,A] = iirdig(fp, fs, Ap As, fm, type)
%
% This function calculates an IIR digital filter of the type specified by variable type
%   type can be 'butter', 'cheby1' or 'cheby2'
%
% fp is a vector of passband frequencies (Hz)
% fs is a vector of stopband frequencies (Hz)
% Ap is the attenuation in db at the passband frequencies
% As is the attenuation in db at the stopband frequencies
% fm is the sampling frequency (Hz)
%


% Check the parameters fp and fs

np = size(fp);
ns = size(fs);

if (~isequal(np,ns))
    fprintf(1,'Error: fp and fs are of different sizes\n'); return;
elseif (min(np) == 1)  % fp and fs are vectors of the same size
    if (max(np) == 1)  % fp and fs are scalars (filter must be either a lowpass or highpass)
        if fp < fs
            filter = 'lp';
        elseif fp > fs
            filter = 'hp';
        else
            fprintf(1,'Error: fs and fp are equal\n'); return;
        end
    elseif (max(np)==2) % fp and fs are vectors (filter must be a bandpass or bandstop)
        if (fs(1) < min(fp)) & (fs(2) > max(fp))
            filter = 'bp';
        elseif (fp(1) < min(fs)) & (fp(2) > max(fs))
            filter = 'bs';
        else 
            fprintf(1,'Error: Values of fp and fs does not correspond to a known kind of filter\n'); return;
        end
    else
        fprintf(1,'Error: The length of fs or fp is bigger than 2\n'); return;
    end 
else
    fprintf(1,'Error: fp, fs or both are matrices\n'); return;
end

% Ordering fp and fs from smaller to bigger frequencies

if filter=='bp' | filter=='bs'
    x = min(fp); y = min(fs);
    if (x==2) tmp = fp(1); fp(1) = fp(2); fp(2) = tmp; end
    if (y==2) tmp = fs(1); fs(1) = fs(2); fs(2) = tmp; end
end


% Check the Attenuation parameters

ap = size(Ap);
as = size(As);

if (~isequal(ap,as))
    fprintf(1,'Error: Ap and As are of different sizes\n'); return;
elseif (ap ~= [1 1])
    fprintf(1,'Error: Ap and As are not scalars\n');return;
elseif (Ap <= 0 | As <= 0)
    fprintf(1,'Error: Ap or As or both are less or equal than 0\n'); return;
end

%% Check the sampling frequency
maxfreq = max(max([fp fs]));

if fm < 2*maxfreq
    fprintf(1,'Warning: fm < twice the maximum frequency specified in in fp and fs\n');
    fprintf(1,'         I suggest to increase the sampling frequency\n');
end


%% Check the filter type

ltype = lower(type);

switch ltype
case {'butter','cheby1','cheby2'}
otherwise
    fprintf(1,'Error: Unknown filter type'); return;
end



%% Prewarp specs so I can apply the bilinear transformation
wp_pw = 2*tan(2*pi*fp/(2*fm));
ws_pw = 2*tan(2*pi*fs/(2*fm));

%% Convert specs to low pass prototype specs

switch filter
case {'lp'}
    wp_lpp = wp_pw;
    ws_lpp = ws_pw;
case {'hp'}
    wp_lpp = ws_pw;
    ws_lpp = wp_pw;
case {'bp'}
    wx2 = prod(wp_pw); wx = sqrt(wx2);
    if prod(ws_pw) < wx2 
        ws_pw(1) = wx2/ws_pw(2);
    else
        ws_pw(2) = wx2/ws_pw(1);
    end
    wp_lpp = wp_pw(2)-wp_pw(1);
    ws_lpp = ws_pw(2)-ws_pw(1);
case {'bs'}
    wx2 = prod(wp_pw); wx = sqrt(wx2);
    if prod(ws_pw) < wx2 
        ws_pw(2) = wx2/ws_pw(1);
    else
        ws_pw(1) = wx2/ws_pw(2);
    end
    wp_lpp = ws_pw(2)-ws_pw(1);
    ws_lpp = wp_pw(2)-wp_pw(1);
end


%% Design the Lowpass prototype

switch ltype
case {'butter'}
    e2 = 10^(0.1*Ap)-1;
    n = ceil(log10(sqrt((10^(0.1*As)-1)/e2))/log10(ws_lpp/wp_lpp));
    wc = (1/e2)^(1/(2*n));
    [Z,P,K] = buttap(n);
    [B,A] = zp2tf(Z,P,K);
case {'cheby1'}
    e2 = 10^(0.1*Ap)-1;
    n = ceil(acosh(sqrt((10^(0.1*As)-1)/e2))/acosh(ws_lpp/wp_lpp));
    [Z,P,K] = cheb1ap(n,Ap);
    [B,A] = zp2tf(Z,P,K);
    wc = 1;
case {'cheby2'}
    e2 = 1/(10^(0.1*As)-1);
    n = ceil(acosh(sqrt(1/(10^(0.1*Ap)-1)/e2))/acosh(ws_lpp/wp_lpp));
    [Z,P,K] = cheb2ap(n,As);
    [B,A] = zp2tf(Z,P,K);
    wc = 1;
end

[H,W] = freqs(B,A,500);
Hdb = 20*log10(abs(H));
figure; subplot(3,1,1), plot(W,Hdb);

switch filter
case {'lp'}
    if (ltype == 'butter') | (ltype == 'cheby1')
        [Ban, Aan] = lp2lp(B, A, wp_lpp*wc);
    else
        [Ban, Aan] = lp2lp(B, A, ws_lpp);
    end
case {'hp'}
    if (ltype == 'butter') | (ltype == 'cheby1')
        [Ban, Aan] = lp2hp(B, A, ws_lpp/wc);
    else
        [Ban, Aan] = lp2hp(B, A, wp_lpp);
    end
case {'bp'}
    if (ltype == 'butter') | (ltype == 'cheby1')
        [Ban, Aan] = lp2bp(B, A, wx, wp_lpp*wc);
    else
        [Ban, Aan] = lp2bp(B, A, wx, ws_lpp);
    end
    Ap = [Ap Ap];
    As = [As As];
case {'bs'}
    if (ltype == 'butter') | (ltype == 'cheby1')
        [Ban, Aan] = lp2bs(B, A, wx, ws_lpp/wc);
    else
        [Ban, Aan] = lp2bs(B, A, wx, wp_lpp);
    end
    Ap = [Ap Ap];
    As = [As As];
end

[Han,W] = freqs(Ban,Aan,500);
Handb = 20*log10(abs(Han));
subplot(3,1,2), plot(W,Handb,[wp_pw ws_pw],-[Ap As], '*r');

%% Bilinear Transformation

[Bd Ad] = bilinear(Ban,Aan,1);
[Hd,F] = freqz(Bd,Ad,500,fm);
Hddb = 20*log10(abs(Hd));
subplot(3,1,3), plot(F,Hddb,[fp fs],-[Ap As],'*r');

[gd,F] = grpdelay(Bd,Ad,500,fm);
figure, plot(F,gd);

fprintf(1,'The filter order is n = %d\n',n);

return;