function Dpoly = idpoly(Dtf, varargin)
% Conversion to polydata.
%
%  Type: 'measured', or 'split'.
%     'measured': Convert to measured model component (A=C=D=1).
%     'split':    Split model inputs so that that the last Ny inputs are
%                 treated as noise inputs. Original model must have Ny or
%                 more inputs.

%   Author(s): Rajiv Singh
%   Copyright 2009-2012 The MathWorks, Inc.

Warn = ctrlMsgUtils.SuspendWarnings('Control:ltiobject:UseSSforInternalDelay'); %#ok<NASGU>

Split = nargin>1 && strncmpi(varargin{1},'s',1);
[Ny, Nu] = iosize(Dtf);
if Split && localPredictorStabCorrectionRequired(Dtf,Ny,Nu)
   % Go via SS to "correct" the predictor, if required.
   Dpoly = idpoly(ss(Dtf),'split');
   return
end

% Note: Ny>1 && Split might be TRUE. But then noisenum, noiseden are
% diagonal.
a = []; c = []; d = []; NV = eye(Ny);
if Split, Nu = Nu-Ny; end

[b, f, negDelay] = localMakeMonic(Dtf.num(:,1:Nu), Dtf.den(:,1:Nu), Dtf.Ts~=0, false);
ioDelay = Dtf.Delay.IO(:,1:Nu)+negDelay;
id = Dtf.Delay.Input(1:Nu,:);
if Split
   [c, d, ~, NV] = localMakeMonic(Dtf.num(:,Nu+1:end), Dtf.den(:,Nu+1:end), Dtf.Ts~=0, true);
end

if any(ioDelay(:)<0)
   ctrlMsgUtils.warning('Ident:transformation:improperDiscreteLTI','idpoly')
   md = min(ioDelay(:));
   id = id + md;
   ioDelay = ioDelay - md;
end

% Absorb output delays into IO delay matrix. REVISIT: error out if output delay?
ioDelay = ioDelay+repmat(Dtf.Delay.Output,[1 Nu]);

% Migrate D and F into A if they are equal
if Split
   a = num2cell(eye(Ny));
   for ky = 1:Ny
      if Nu==0 || idpack.isequalwithreltol([f(ky,:),d(ky)],1e-6)
         a{ky,ky} = d{ky};
         f(ky,:) = num2cell(ones(1,Nu));
         d{ky} = 1;
      end
   end
   
   if isequal(a,num2cell(eye(Ny))), a = []; end
   if isequal(f,num2cell(ones(Ny,Nu))), f = []; end
   if Nu==0 || isequal(d,num2cell(ones(Ny,1))), d = []; end
end

PolyStruct = pmodel.polynomial(a,b,c,d,f,ioDelay);
% Fix leading zeros of B polynomial in discrete-time case.
PolyStruct = pmodel.polynomial.fixNkInB(PolyStruct,Dtf.Ts);

Dpoly = idpack.polydata(PolyStruct,Dtf.Ts);
Dpoly.Delay.Input = id;
Dpoly.NoiseVariance = NV;
Dpoly = setLastOperation(Dpoly,'tf');

%--------------------------------------------------------------------------
function [num, den, negDelay, nv] = localMakeMonic(num, den, isDT, isNoise)
% Make denominator polynomial monic in the transfer function num/den.
% Return zero access over pole as negative delay for discrete-time data.

if ~iscell(num)
   num = {num};
   den = {den};
end

[ny,nu] = size(num);
negDelay = zeros(ny,nu);
nv = eye(ny);
if isNoise
   NoiseInd = 1:ny+1:ny^2;
end

for ct = 1:ny*nu
   if isNoise && ~any(NoiseInd==ct)
      continue;
   end
   numct = num{ct};
   denct = den{ct};
   nr = find(numct,1,'first');
   dr = find(denct,1,'first');
   if isNoise && nr~=dr
      ctrlMsgUtils.error('Ident:transformation:NoiseTF2idpoly')
   end
   if isempty(nr)
      nr = min(length(numct),dr);
   end
   % Remove leading zeros and normalize 
   numct   = numct(min(nr,dr):end)/denct(dr);
   den{ct} = denct(dr:end)/denct(dr);
   if isDT      
      negDelay(ct) = min(nr-dr,0);
      % remove trailing zeros from numerator
      nr = find(numct,1,'last');
      if ~isempty(nr)
         numct = numct(1:nr);
      end
   else
      % remove leading zeros from numerator
      % (assume them to be fillers)
      nr = find(numct,1,'first');
      if ~isempty(nr)
         numct = numct(nr:end);
      end
   end
   
   if isNoise
      noisestd = 1;
      if numct(1)~=1
         noisestd = numct(1);
         numct = numct/noisestd;
      end
      nv(ct) = noisestd^2;
   end
   num{ct} = numct;
end

if isNoise
   num = num(NoiseInd); num = num(:);
   den = den(NoiseInd); den = den(:);
end
   
%--------------------------------------------------------------------------
function boo = localPredictorStabCorrectionRequired(Dtf,Ny,Nu)
% Check if idss predictor stabilization needs to be used.
Nu = Nu-Ny;
nnum = Dtf.num(:,Nu+1:end); Ts = Dtf.Ts;
boo = false;
for ky = 1:Ny
   for ky2 = 1:Ny
      if ky==ky2
          r = roots(nnum{ky,ky});
          if ~isempty(r)
            boo = (Ts~=0 && any(abs(r)>1)) || (Ts==0 && any(real(r)>0));
          end
      else
         % for nondiagonal C/D, require predictor correction
         boo = norm(nnum{ky,ky2})>0;
      end
      if boo, return, end
   end
end
