estorderiv.m

function[N0,nA0,nB0,d0]=estorderiv(y,u,nmax,sel)
% Function to estimate the orders of a discrete time system.
% [N0,nA0,nB0,d0]=estorderiv(y,u,nmax,sel)
% 
% Inputs:
% y = output data vector
% u = input data vector
% nmax = maximum system order (default 20)
% If sel is 1, the user will be asked to select the order N of the system 
% based on the order estimation curve.
% By default, sel=0. The value of N, minimizing the estimation criterion
% will be used.
% 
% Outputs:
% N0 = global order of the system
% nA0 = order of the denominator of the system
% nB0 = order of the numerator of the system
% d0 = identified delay
%
% Written by:  T.B. Airimitoaie, I.D. Landau and H. Prochazka
% 17th december 2009
% Last updated 15th february 2011

if nargin<2, error('"estorderivg" needs a minimum of 2 parameters. Write "help estorderivg" to the console for more informations.'); end
if nargin==2, nmax=20; sel=0; end
if nargin==3, sel=0; end

[V,S,VS]=estorderN(y,u,nmax);
[VS V];
figure
plot(0:nmax,VS,'r',0:nmax,V,'b')
title('Instrumental Variables estimation of system order (N)')
if sel
    N0=input('Please input here the order of the system N = ');
else
    [vmin,imin]=min(VS);
    N0 = imin-1;
    disp(['N = ' int2str(N0)]);
end

[V,S,VS]=estdiv(y,u,N0);
[VS V];
figure
plot(0:length(VS)-1,VS,'r',0:length(V)-1,V,'b')
title('IV estimation of N-d')
[vmin2,imin2]=min(VS);
d0 = N0-imin2+1;
disp(['d = ' int2str(d0)]);

[V,S,VS]=estnBiv(y,u,N0,d0);
[VS V];
figure
plot(0:length(VS)-1,VS,'r',0:length(V)-1,V,'b')
title('IV estimation of n_B+d')
[vmin2,imin2]=min(VS);
nB0 = imin2-1-d0;
disp(['nB = ' int2str(nB0)]);

[Va,Sa,VSa]=estnAiv(y,u,N0, d0, nB0);
[VSa Va];
figure
plot(0:length(VSa)-1,VSa,'r',0:length(Va)-1,Va,'b')
title('IV estimation of n_A')
[vmin2,imin2]=min(VSa);
nA0 = imin2-1;
disp(['nA = ' int2str(nA0)]);

function [VLs,S,VS]=estorderN(y,u,nmax)
if nargin==3,
elseif nargin==4, 
else return;
end;

sz=size(y);
if(sz(2)~=1), y=y';end;%column
sz=size(u);
if(sz(2)~=1), u=u';end;

VLs=[];
S=[];
VS=[];
L=2*nmax;
N=length(y)+1-L;
%construction of Z and y0
Z=[];
for k=1:2*nmax
    Z=[Z u(2*nmax-k+1:2*nmax-k+N-1)];
end
y0=y(L+1:length(y));

%QR-triangularization
[Q,R]=qr(Z);
Z1=R(1:L,:);
Zm=Z1'\Z';
Y0=Zm*y0;
for nn=0:nmax
    %construction of R(nn)
   Rnny=[];
   Rnnu=[];
   for k=1:nn
      yp=y(L-k+1:L-k+N-1);
      up=u(L-k+1:L-k+N-1);
      Rnny=[Rnny yp];
      Rnnu=[Rnnu up];
   end;
   if nn==0
       R1nn = 0;
       %theta computation
       theta_min=0;%inv(R1nn)*Y0;
   else
       R1nn=Zm*[Rnny Rnnu];
       %theta computation
       theta_min=(R1nn'*R1nn)\R1nn'*Y0;%inv(R1nn)*Y0;
   end
   criter=sum((Y0-R1nn*theta_min).^2)/N;
   VLs=[VLs;criter];
   penal=0.4*nn*log10(N)/N;
   S=[S;penal];
end
VLs=VLs/VLs(1);
VS=VLs+S;
VS=VS/VS(1);

function [VLs,S,VS]=estdiv(y,u,nmax)
if nargin==3,
elseif nargin==4, 
else return;
end;

sz=size(y);
if(sz(2)~=1), y=y';end;%column
sz=size(u);
if(sz(2)~=1), u=u';end;

VLs=[];
S=[];
VS=[];
L=2*nmax;
N=length(y)+1-L;   
%construction of R(nn) and Z
Z=[];
uz=u;
for k=1:2*nmax
    Z=[Z uz(2*nmax-k+1:2*nmax-k+N-1)];
end;
y0=y(L+1:length(y));

%QR-triangularization
[Q,R]=qr(Z);
Z1=R(1:L,:);
Zm=Z1'\Z';
Y0=Zm*y0;
for nn=0:nmax
   Rnny=[];
   Rnnu=[];
   for k=1:nmax
      yp=y(L-k+1:L-k+N-1);
      Rnny=[Rnny yp];
   end
   for k=1:nn
      up=u(L+(-nmax+k-1)+1:L+(-nmax+k-1)+N-1);
      Rnnu=[Rnnu up];
   end
   R1nn=Zm*[Rnny Rnnu];
   if size(R1nn)==[0 0], R1nn=1; end
	%theta computation
   theta_min=inv(R1nn'*R1nn)*R1nn'*Y0;%inv(R1nn)*Y0;
   criter=sum((Y0-R1nn*theta_min).^2)/N;
   VLs=[VLs;criter];
   penal=0.4*nn*log10(N)/N;
   S=[S;penal];
end
VLs=VLs/VLs(1);
VS=VLs+S;
VS=VS/VS(1);

function [VLs,S,VS]=estnBiv(y,u,nmax,d0)
if nargin~=4
    return;
end

sz=size(y);
if(sz(2)~=1), y=y';end;%column
sz=size(u);
if(sz(2)~=1), u=u';end;

VLs=[];
S=[];
VS=[];
L=2*nmax;
N=length(y)+1-L;   
%construction of R(nn) and Z
Z=[];
uz=u;
for k=1:2*nmax,
    Z=[Z uz(2*nmax-k+1:2*nmax-k+N-1)];
end;
y0=y(L+1:length(y));

%QR-triangularization
[Q,R]=qr(Z);
Z1=R(1:L,:);
Zm=Z1'\Z';
Y0=Zm*y0;

for nn=0:nmax
   Rnny=[];
   Rnnu=[];
   for k=1:nmax
      yp=y(L-k+1:L-k+N-1);
      Rnny=[Rnny yp];
   end
   for k=1:nn
      up=u(L-k+1:L-k+N-1);
      Rnnu=[Rnnu up];
   end
   R1nn=Zm*[Rnny Rnnu];
   if size(R1nn)==[0 0], R1nn=1; end
	%theta computation
   theta_min=inv(R1nn'*R1nn)*R1nn'*Y0;%inv(R1nn)*Y0;
   criter=sum((Y0-R1nn*theta_min).^2)/N;
   VLs=[VLs;criter];
   penal=0.4*nn*log10(N)/N;
   S=[S;penal];
end
VLs=VLs/VLs(1);
VS=VLs+S;
VS=VS/VS(1);

function [VLs,S,VS]=estnAiv(y,u,nmax,d0,nB0)
if nargin~=5
    return;
end

sz=size(y);
if(sz(2)~=1), y=y';end;%column
sz=size(u);
if(sz(2)~=1), u=u';end;

VLs=[];
S=[];
VS=[];
L=2*nmax;
N=length(y)+1-L;   
%construction of R(nn) and Z
Z=[];
for k=1:2*nmax,
    Z=[Z u(2*nmax-k+1:2*nmax-k+N-1)];
end;
y0=y(L+1:length(y));

%QR-triangularization
[Q,R]=qr(Z);
Z1=R(1:L,:);
Zm=Z1'\Z';
Y0=Zm*y0;

Rnnu=[];
for k=d0+1:d0+nB0
  up=u(L-k+1:L-k+N-1);
  Rnnu=[Rnnu up];
end

for nn=0:nmax
   Rnny=[];
   for k=1:nn
      yp=y(L-k+1:L-k+N-1);
      Rnny=[Rnny yp];
   end
   R1nn=Zm*[Rnny Rnnu];
   if size(R1nn)==[0 0], R1nn=1; end
   %theta computation
   theta_min=inv(R1nn'*R1nn)*R1nn'*Y0;%inv(R1nn)*Y0;
   criter=sum((Y0-R1nn*theta_min).^2)/N;
   VLs=[VLs;criter];
   penal=0.4*nn*log10(N)/N;
   S=[S;penal];
end
VLs=VLs/VLs(1);
VS=VLs+S;
VS=VS/VS(1);