function [ahat, bhat, L] = batlrpi(y, U, N, D, T) 
%
% BATLRPI Gives estimated parametersof a model of the type
% Ay=Bu+Der.
% Function Call:
%   [ahat, bhat, L] = batlrpi(y, U, N, D, T)
% Inputs :
% y and U are vectors of sequential output and input
% respectively assumed t o be of the same length.
%   N(1) = na number of estimated A parameters
%   N(2) = nb number of estimated B parameters
%   N(3) = k the assumed time-delay in samples(>=1) 
%   N(4) = N1 Initial Costing Horizon
%   N(5) = N2 Final Costing Horizon
% D, T are D and T in the above model
% Outputs:
% ahat, bhat are A and B in the above model
% L is the converged filter for LRPI.
%
%
na=N(l);nb=N(2); k=N(3);N1 =N(4); N2=N(5); y=y(:); u=u(:);
to1 = l e - 3 ;
err = 1;
Nmax = length(y);
theta = zeros(na + nb, 1); 
while err > tol,
  if na > 0,
    L-eqivl([l theta(nb+l:na+nb)‘], D,T, N1, N2); 
  else
    L=eqivl([l 01, D, N1, N2); 
  end;
  yf= filter(conv(L, D), T, y);
  uf = filter(conv(L, D), T, U); 
  n=max(na, nb+k-1); 
  X=uf(n-k+l:Nmax-1 -k+l); 
  for i= 2: nb,
    X=[Xuf(n-i+1 -k+l:Nmax-i-k+l)]; 
  end
  for i = 1: na,
    X=[X-yf(n -i+1:Nmax-i)];
  end
  Y=yf(n+l:Nmax);
  thetanew = X\Y;
  err = norm(theta - thetanew, 1)/max(norm(theta, 1), 1e-6); 
  theta = thetanew;
end
bhat = [zeros(l, k) theta(1 : nb)’];
ahat= [ltheta(nb+1:na+nb)’];
return

function [L] = eqivl(a, D, T, N1, N2)
%
% EQlVL Computes the equivalent L prefilter as
%
% L(q)L(q^-1) =sum over j {Ni -> N2) of Ej(q)Ej(q^-1) 
% Function Call:
%   L=eqivl(A,D,T,N1,N2)
% Inputs:
% A, D, T are polynomials in the equation:
%   T = EjAD + q^jF_j
%   N1 and N2 are the limiting values of j.
% Output:
%   L finite pulse responsefilter with the above property 
%   with unity steady state gain
%
E=filter(T, conv(a, D), [1,zeros(1, N2)]); 
m=conv(E(1:N1), E(N1:-1:1));
for i= N1+1:N2
  mi=conv(E(1:i),E(i:-1:1));
  m= [ 0 m0 ] + mi; 
end;
L=spctr(m); 
L=L/sum(L); 
return


function [x] = spctr(m, tol)
%SPCTR Performs Spectral factorisation using Peterka's method. Function Call:
% X = spctr(M) 
% Input:
%   M : A polynomial which is going to be split into stable and unstable
%   bits. It is assumed that the polynomial of the form:
%
%   M(q) = N(q)N(q^-1)
%   tol' Tolerance limit for accuracy of the convergence.
% Output:
%   X: Polynomial which is the spectral factor of M.
%
if nargin==1,
  tol=1e-3,
end; 
Im = length(m);
x = [ 1 zeros(1, fix(lm/2))];
r = 1;
while r > to1
  [x, rem] = deconv(m, x); 
  lx = length(x);
  x=x(lx: -1: 1)/x(lx);
  r = norm(rem, 1) ;
end
return

function [r, s, pc] = igpc(a, b, n1 , n2, nu, lambda, t, pnz)
%  A PC-MATLAB macro to calculate the transfer function form of the discrete
%  time GPC Controller. The computations are done via the state-space as it
%  is generally easier to do so. The function call is:
%
%  [r, s, pc] = igpc(a, b, nl, n2, nu, lambda, t, pnz) 
%
% Inputs:
%
%  a and b are the A(q^-1) and B(q^-1) polynomials
%  n1 is the initial costing horizon
%  n2 is the final costing horizon
%  nu is the control horizon lambdaisthecontrolcostingparameteronwhat's left t is the perenial observer polynomial!

pnz is the P of Py auxiliary output outputs:
rand s are the coefficients of w(t) - S(q-')y(t) = R(q-')u(t) pc is the closed-loop pole polynomial
setupthematricesfortheRiccatiequationupdate
iflambda==0,lambda=eps;end
adelta =conv(a, [ l-11);
I= length(b); bpol= b(2:/);
np = max(size(pnz)); pnz = pnz/sum(pnz); pt = conv(pnz, t); npt = max(size(pt));
na = max(size(ade1ta));
ifna<npt+1,adelta= [adeltazeros(1,npt-na+1)];end;
na = max(size(ade1ta));
nb = max(size(bpo1)); if nb < na. bpol= bpol zeros(1, na - nb)]; end ;
