function [Vfpe,Vaic,Vmdl] = select_struct_ARX(Yvet,N,OrdMax)
% struct_selectARX( Yvet , N , OrdMax ) visualizza le funzioni costo per la
% determinazione dell'ordine di un modello ARX MISO secondo i criteri FPE , 
% AIC e MDL. Devono essere specificati il numero N di campioni di ingresso 
% e uscita nella matrice Yvet e l'ordine massimo OrdMax che deve essere
% testato. Il modello puo' essere MISO.

L = length(Yvet);

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

  [ nr , nc ] = size( Yvet );

  Vfpe=zeros(1,OrdMax); Vaic=zeros(1,OrdMax); Vmdl=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;

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

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

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


           end; %for OrdMax 


   k=OrdMax;   
   figure, subplot(311), % grafico FPE.
   plot( [1:k] , Vfpe, '-', [1:k] , Vfpe , 'o'); [Vmin1,iVfpe] = min(Vfpe);
   title('FPE'), xlabel(sprintf('Model Order. Minimul value of FPE = %f for n = %d. (n_a = n_b = n, n_k = 1).',Vmin1,iVfpe))

   subplot(312); % grafico AIC.
   plot( [1:k] , Vaic , '-', [1:k] , Vaic,'o'); [Vmin2,iVaic] = min(Vaic);
   title('AIC'), xlabel(sprintf('Model Order. Minimul value of AIC = %f for n = %d. (n_a = n_b = n, n_k = 1).',Vmin2,iVaic))

   subplot(313); %  grafico MDL.
   plot( [1:k] , Vmdl,'-',[1:k] , Vmdl,'o'); [Vmin3,iVmdl] = min(Vmdl);
   title('MDL'), xlabel(sprintf('Model Order. Minimul value of MDL = %f for n = %d. (n_a = n_b = n, n_k = 1).',Vmin3,iVmdl))

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