function [pp1,pp2,pp3,pp4] = struct_selectARX_val2( Yiu1 , N , ordmax )

[Lt,nc] = size(Yiu1); % Lunghezza effettiva della sequenza.

L = floor(Lt/2); % Set di dati per l'identificazione

N = L - ordmax-8;

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

for i = 1:ordmax,

 [Theta,Ji,sigma_e,cov_e,err_e,yprev] = minquad3(Yiu1(1:L,:),N,i,0);
   
 [zeta_tmp,epserr] = bianres3( Yiu1(L+1:2*L,:) , N , 8 , i , Theta' , 0 );
 
 sk = (epserr' * epserr)/length(epserr);
 
 clear Theta;

 pp1(i) = 100 * sqrt( sk/( [Yiu1(L+1+i:2*L+i- ordmax-8,nc)]'*Yiu1(L+1+i:2*L+i- ordmax-8,nc) ) );

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

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

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

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  
  
%%%%%%%

function [Theta,J,sigma_e,cov,err,Yprev] = minquad3( Yvet , Nseq , Ord , graph )

  [ nr , nc ] = size( Yvet );    % ci sono ingressi ed uscite;

  if( nc >= 2 ) , %  Prime colonne ingressi e ultima col. uscita.

  H=[ myhank(Yvet(:,nc),Nseq+Ord,Ord,0), myhank(Yvet(:,1),Nseq+Ord , Ord ,0 )];

  else % Nessun ingresso, solo una uscita. 
  
  H= myhank( Yvet(:,1) , Nseq+Ord , Ord , 0 ); % Yvet e' un vettore.

  end;

  if( nc > 2 ),

    for input = 2:nc-1,
     
    H = [ H , myhank( Yvet(:,input) , Nseq+Ord , Ord , 0 )];

    end;

  end;
   
     Theta = pinv( H ) * Yvet( Ord+1:Nseq+Ord , nc) ;
     
   %  Theta = ( inv( H'* H ) )*( H')*Yvet( Ord+1:Nseq+Ord , nc) ;

if (graph == 1),


figure
plot([1:Nseq],H*Theta,'-',[1:Nseq],Yvet( Ord+1:Nseq+Ord ,nc ),'--');
title('Simulazione linea cont. e uscita misurata linea tratt.')
ylabel('Y (-) simulata e (--) misurata'),xlabel('Campioni'), pause;

%print -dps simPe.ps

end;

  
  Yprev = H*Theta;
  errps = Yvet( Ord+1:Nseq+Ord ,nc )-Yprev;

if (graph == 1),

figure
plot([1:Nseq],errps,'-');
title('Errore di previsione')
ylabel('Yv-Yp'),xlabel('Campioni'), pause;

%print -dps errPe.ps

end;


   J = ( errps'* errps )/Nseq;

   sigma_e = ( J * Nseq )/( Nseq - nc*Ord );

   cov = sigma_e * pinv( H' * H );

   err = errps;


return

%%%%

function [zeta,epserr] = bianres3( Y , Ncov , Mgdl , Ord , Theta , graph )

  [ L , nc ] = size( Y ); %L numero di campioni;

  % Ncov lunghezza della sequenza considerata;

  Yo = Y(Ord+1:L,nc);
  
        if( nc >= 2 )

    H = [ myhank( Y(:,nc) , L , Ord , 0 ) , myhank( Y(:,1),L,Ord,0 ) ];

        else H = myhank(  Y(:,1) , L , Ord , 0 );

        end;


        if( nc > 2)

          for input = 2:nc-1,

              H = [ H , myhank(  Y(:,input) , L , Ord , 0 ) ];

          end;

        end;

  Yp = H*Theta'; %vettore delle uscite previste;

if( graph == 1 )
  figure, subplot(211)
  plot([0:Ncov-1],Yo(1:Ncov),'-',[0:Ncov-1],Yp(1:Ncov),':')
% plot([0:L-Ord-1],Yo,'-',[0:L-Ord-1],Yp(1:L-Ord),':')
  title('Prediction (:) and observations (-)')
  xlabel('Samples'), ylabel('y(t)')
end;

  epserr=[]; %vettore dei residui;

  res=[]; %vettore delle covarianze campionarie;

  epserr=Yo-Yp(1:L-Ord);

if( graph == 1 )
  figure, subplot(212)
  plot([0:Ncov-1],epserr(1:Ncov)),title('Identified model residuals')
  xlabel('Samples'), ylabel('\epsilon(t)'), pause;


end;

           for i=1:Mgdl+1,

             res(i) = ( (epserr(1:Ncov))' * epserr(i:Ncov+i-1) )/Ncov;

           end;

if( graph == 1 )
     figure, plot([1:Mgdl+1],res,'-',[1:Mgdl+1],res,'o');
     title('Residual whiteness test'), xlabel('Correlation (\tau)'), pause;
end;

  zeta = Ncov *( res(2:Mgdl+1)*( res(2:Mgdl+1)' ) )/(res(1)*res(1));

return