function [mp,m] = tfest(data,nn,varargin)
% TFEST: Estimate Continuous Time Transfer function
%
%   M = TFEST(Data,[nb nf]) or M = TFEST(Data,Mi)
%
%   M: Y = B/F U  MIMO
%
%   If Mi is specified as an IDPOLY object, it will be used at the initial
%   model.
%   If the orders nb of B and nf of F (ny-by-nu matrices) (nb<=nf) are
%   specified, an initial estimate will be formed using OE. By 
%   
%   With
%   [M,T] = TFEST(Data,[nb nf],P1,V1,...,Pn,Vn)
%   further Property/Value pairs for the IDTF object M can be set, including
%   algorithm options like SearchMethod, MaxIter, Tol. See
%   IDPROPS/... for more details

[dum,Ny,Nu]=size(data);
if isa(nn,'idmodel')     
    mi=nn;
else
    [nmy,nc]=size(nn);
    if nmy~=Ny
        error('Ident:tfest:wrongny',...
            'The orders must have the same number of rows as the number of outputs.')
    end
    if nc ~= 2*Nu
        error('Ident:tfest:wrongorder',...
            'The orders should be [nb nf] with nb and nf as Ny-by-Nu matrices')
    end
    nb = nn(:,1:Nu); 
    nf = nn(:,Nu+1:end);
    if any(any(nb>nf+1))
        error('Ident:tfest:toolargenb',...
            'The numerator orders may exceed the denominator orders by at most 1.')
    end
    ts = pvget(data,'Ts');
    
    mi = tfinivaloe(data,nn);
     
end
ms = idp2ids(mi);
m = pem(data,ms,'ts',0,'sea','gna',varargin{:});
mp = localoepoly(m,mi.nb,mi.nf);
mp = inherit(mp,m);
%************************************************************
function ms=idp2ids(mp);
nf = mp.nf;
nb = mp.nb;
n = sum(nf);
f=mp.f;
b=mp.b;
[ny,nu]=size(mp);
A=[]; B=[]; C=[]; As =[]; Bs=[]; 
for ku =1:nu
    ak = [-f(ku,2:nf(ku)+1)', eye(nf(ku),nf(ku)-1)];
    bk = [zeros(nf(ku)-nb(ku),1);b(ku,1:nb(ku))'];
    ask = ak;
    ask(1:nf(ku),1)= NaN;
    bsk = [zeros(nf(ku)-nb(ku),1);NaN*ones(nb(ku),1)];
    [na]=size(A,1); [nbr,nbc]=size(B);
    A = [[A,zeros(na,nf(ku))];[zeros(nf(ku),na),ak]];
     As = [[As,zeros(na,nf(ku))];[zeros(nf(ku),na),ask]];
     B = [[B,zeros(nbr,1)];[zeros(nf(ku),nbc),bk]];
             Bs = [[Bs,zeros(nbr,1)];[zeros(nf(ku),nbc),bsk]];
             C = [C,1,zeros(1,nf(ku)-1)];
end
Cs = C;
ms = idss(A,B,C,zeros(1,nu),zeros(n,1),zeros(n,1),'As',As,'Bs',Bs,...
    'Cs',C,'Ds',zeros(1,nu),'Ks',zeros(n,1),'X0s',zeros(n,1),'ts',0);
%**************************************************************************

        

function mp = tfinivaloe(data,nn,mode,lambda)
%CTTFINIVAL: Initial estimate of CTTF parameters
%
%   M = TFINIVALOE(DATA,[nNum nNum],MODE,Lambda)
%   DATA: Iddata object.
%   nNum: Numerator orders
%   dNum: Denominator orders
% This version for MISO and strictly proper CT TFs
 
[N,ny,nu]=size(data);
 
nb = nn(:,1:nu);
nf = nn(:,nu+1:end);
for ky = 1:ny
    mdn = oe(data(:,ky,:),[nf,nf,ones(size(nf))]);
end
 
mc = d2c(mdn,'cov','none');
b = mc.b;
for ku=1:nu
    b(ku,1:nf(ku)-nb(ku))=zeros(1,nf(ku)-nb(ku));
end
mc.b =b;
mp = mc;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mp = localoepoly(m,nb,nf)
[Ny,Nu]=size(m);
r=1;
as = m.a;
bs = m.b;
ds = m.d;
 cov = pvget(m,'CovarianceMatrix');
 if ~isempty(cov)&&~ischar(cov)
     covp = localcovpoly(m,cov,nb,nf);
 else
     covp = cov;
 end
for ky=1:Ny
    csu= cumsum([r,nf(ky,:)]);
    for ku =1:Nu
        f{ky,ku}=[1 -as(csu(ku):csu(ku+1)-1,csu(ku)).'];
        b{ky,ku}=bs(csu(ku):csu(ku+1)-1,ku).';
        if ds(ky,ku)~=0
            b{ky,ku} =[0 b{ky,ku}]+ds(ky,ku)*f{ky,ku};
        end
    end
    r = csu(Nu+1);
end
b1 =[]; f1 =[];
for ku=1:Nu
    b1=[b1;b{1,ku}];
    f1 = [f1;f{1,ku}];
end
mp = idpoly(1,b1,1,1,f1,'ts',0);     
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function covp = localcovpoly(m,cov,nb,nf)
d = pvget(m,'D');
p = pvget(m,'ParameterVector');

if norm(d,1)==0
par = [1:length(p)]';
m=pvset(m,'ParameterVector',par,'CovarianceMatrix',[]);
mp = localoepoly(m,nb,nf);
pp = pvget(mp,'ParameterVector');
 % so far no Dterms
covp = cov(abs(pp),abs(pp)).*(sign(pp)*sign(pp)');
else
    m.cov = [];
    mp = localoepoly(m,nb,nf);
    pp = pvget(mp,'ParameterVector');
    for kp = 1:length(p)
        p1 =p;
        p1(kp) = p1(kp)-1;
        m1=parset(m,p1);
        mp1=localoepoly(m1,nb,nf);
        pp1=pvget(mp1,'ParameterVector');
        L(:,kp) = pp1-pp;
    end
    covp = L*cov*L';
    
end




