function [pp1,pp2,pp3,pp4] = struct_selectARX( Yvet , N , OrdMax )
% struct_selectARX( Yvet , N , OrdMax ) contiene i criteri per la determinazione
% dell'ordine PPCRE , FPE , AIC e MDL di un modello ARX , di cui siano
% specificati N campioni di ingresso e uscite nella matrice Yvet ; i valori
% dei quattro criteri sono graficati per un modello di ordine massimo OrdMax.
% Il modello puo' essere MISO.

L = length(Yvet);

if( N+OrdMax > L), N = L - OrdMax; end;

  [ nr , nc ] = size( Yvet );

  pp1=zeros(1,OrdMax); pp2=zeros(1,OrdMax); pp3=zeros(1,OrdMax); pp4=zeros(1,OrdMax);


           for k= 1:OrdMax,

               Hyi = myhank( Yvet(:,nc) , N + k , k , 0 );

               Hui =  myhank( Yvet(:,1) , N + k , k , 0 );

               if( nc == 1 ),  Hui = [];
               end;

             if( nc > 2 )               
             
               for input = 2:nc-1,

                   Hui = [ Hui , myhank( Yvet(:,input) , N + k , k , 0 ) ];

               end;

             end;
 
        den = det( [Hyi,Hui]' * [ Hyi,Hui] );
 
        if(den == 0.0), den = 1.0; end;

 sk  = (det([Hyi,Hui ,Yvet(k+1:N+k,nc)]' * [Hyi,Hui,Yvet(k+1:N+k,nc) ] ) / den ) ;

               clear Hui; clear Hyi;

           pp1(k) = 100 * sqrt( sk/( [Yvet(k+1:N+k,nc)]'*Yvet(k+1:N+k,nc) ) );

           pp2(k) = sk*( (N+nc*k)/(N*(N-nc*k)) ) ;

           pp3(k) = 2*nc*k + N*log( sk/N ) ;

           pp4(k) = nc*k*log(N) + N*log( sk/N )  ;


           end; %for OrdMax 


   k=OrdMax;   
   figure, subplot(221); % grafico PPCRE.
   plot( [1:k], pp1(1:k),'-',[1:k] , pp1(1:k) , 'o' );
   title('PPCRE %'), xlabel('Model Order')

   subplot(222); % grafico FPE.
   plot( [1:k] , pp2, '-', [1:k] , pp2 , 'o');
   title('FPE'), xlabel('Model Order')

   subplot(223); % grafico AIC.
   plot( [1:k] , pp3 , '-', [1:k] , pp3,'o');
   title('AIC'), xlabel('Model Order')

   subplot(224); %  grafico MDL.
   plot( [1:k] , pp4,'-',[1:k] , pp4,'o');
   title('MDL'), xlabel('Model Order')

return

%%%
%%% Local functions
%%%

function Hn=myhank(Yvet,Nseq,Ord,Off)

   for i=1:Ord,

      Hn(:,i) = Yvet( Off+i: Nseq-Ord+Off+i-1 );

  end;

return