function [yhat,thpred]=predict(z,theta,m)
%PREDICT Computes the m-step ahead prediction
%	
%	YP = predict(Z,TH,M)
%	
%	Z : The output - input data for which the prediction is computed.
%	    Z = [Y U], where Y and U are column vectors (one column for each 
%	    output and input)
%	TH: The model in the THETA format (see help theta)
%	M : The prediction horizon. Old outputs up to time t-M are used to
%	    predict the output at time t. All relevant inputs are used.
%	    M = inf gives a pure simulation of the system.(Default M=1).
%	YP: The resulting predicted output.
% 
% 	With [YP,THPRED] = predict(Z,TH,M) the predictor is returned in the
%	THETA format
%	See also COMPARE

%	L. Ljung 10-1-89
%	Copyright (c) 1989-90 by the MathWorks, Inc.
%	All Rights Reserved.

if nargin<3, m=1;end
if isempty(m), m=1;end
if isthss(theta), eval('yhat=predicts(z,theta,m);')
if nargout==2,disp('The predictor will not be returned for a theta-model that is defined via a state space model'),end,return,end
[a,b,c,d,f]=th2poly(theta);
nu=theta(1,3);
T=gett(theta);
if T<0,error('This is a continuous time model of input-output type. Please sample it (using thc2thd) before applying predict!'),end

if nu>0,ff=1;for ku=1:nu, 
bt=b(ku,:);for kku=1:nu,if kku~=ku,bt=conv(bt,f(kku,:));end,end
bb(ku,:)=bt;
ff=conv(ff,f(ku,:));end
a=conv(conv(a,ff),d);c=conv(c,ff);else a=conv(a,d);end


na=length(a);nc=length(c);nn=max(na,nc);
a=[a,zeros(1,nn-na)];c=[c,zeros(1,nn-nc)];
[f,g]=deconv(conv([1 zeros(1,m-1)],c),a);
ng=length(g);

if nu>0,
df=conv(d,f);

for ku=1:nu
bf(ku,:)=conv(bb(ku,:),df);
end
nbf=length(bf(1,:));nn=max(ng,nbf);
gg=[[g,zeros(1,nn-ng)];[bf,zeros(nu,nn-nbf)]];
else gg=g;end
thpred=poly2th(c,gg);
yhat=idsim(z,thpred);

