designIIR(phi, D, m0, h, M)

function [num, den, gamma] = designIIR(phi, D, m0, h, M)

% designIIR: design filters F_i(z) (as in Section III.B.)
%
% Usage:    [num, den, gamma] = designIIR(phi, D, m0, h, M)
%
% INPUTS: 
%   phi: input rational transfer functions
%   D:   vector of fractional delays
%   m0:  system delay tolerance
%   h:   fast sampling interval
%   M:   superresolution factor (integer)
%
% OUTPUTS: 
%   gamma:  the H infinity norm of the induced error system K
%   num: numerator vectors
%   den: denominator vectors
%
% filter F_i(z) will be an IIR filter with num{i} is the coefficients of
% the numerator and den{i} is the coefficients of the denominator.
%
% See also: getF, demo

% Get the ingeter and residual of the delays
m = floor(D(:)/h);
d = D(:) - m * h;

% Get the digital system as in Prop. 1 and 2
Ad = getAd(phi, h, d);
Bd = getBd(phi, h, d);
Cd = getCd(phi, h, d);
Dd = zeros(size(Cd,1), size(Bd,2));

% The Integer Delay Operator
sysd = IntDelayOp([m0; m]);

% Get the digital system as in Prop. 3 by taking into account the integer
% delay operators
sys = sminreal( sysd * ss(Ad, Bd, Cd, Dd, -1) );    % sminreal to reduce the system's dimension on the fly 
% sys = sysd * ss(Ad, Bd, Cd, Dd, -1);

% Get the system P (see Fig. 6)
[Ap, Bp, Cp, Dp, p, q] = getP(sys, M);

% Convert P into analog, since hinfsyn works in analog domain
[Ac, Bc, Cc, Dc] = Digital2Analog(Ap, Bp, Cp, Dp);

% Design in continuous time
P = pck(Ac, Bc, Cc, Dc);

ubd = 1;
lbd = 0;
inc = .0001;

[F, g, gamma] = hinfsyn(P,p,q,lbd,ubd,inc);

% Get the synthesis system F
[Af, Bf, Cf, Df] = unpck(F);

% Convert F back to digital domain
[Af, Bf, Cf, Df] = Analog2Digital(Af, Bf, Cf, Df);

% Convert to filter coefficients
[num, den] = getF(Af, Bf, Cf, Df, M);