   clear;
   clc;
   close all
    
   echo on
   % Rao-Garnier test experiment used to generate Figure 1 in the following paper
   % H. Garnier, Direct continuous-time approaches to system identification.
   % Overview and benefits for practical applications.
   % European Journal of Control, July 2015
   
   % You need to download Version 7.0 of the CONTSID toolbox
   % See www.univ-lorraine.fr/contsid/
   echo off
   
    
    % Rao-Garnier benchmark system definition
    w1=20;
    z1=0.1;
    w2=2;
    z2=0.25;  
    K=1;
    [num1,den1]=ord2(w1,z1);
    [num2,den2]=ord2(w2,z2);
    Nc=conv(K*num1*num2*w2^2*w1^2,[-4 1]);
    Dc=conv(den1,den2);

    %  The sampling period is chosen to be 0.01;
    Ts=0.01;	
    
    M0=idpoly(1,Nc,1,1,Dc,'Ts',0);

    u = prbs(10,7);
    datau = iddata([],u,Ts,'InterSample','zoh');
       
    % Add some Gaussian white noise to the output measurement
    snr=10;
    y = simc(M0,datau,snr);
    data = iddata(y,u,Ts,'InterSample','zoh');

    % Estimate a continuous-time transfer function model
    % by using the CONTSID SRIVC routine from the noisy data object
    % The extra-information needed are:
    nn=[2 4 0];
    Msrivc = srivc(data,nn);

    % Estimate a discrete-time transfer function or state-space model
    % by using the SID ARX, OE and N4SID routines from the noisy data
    % object
    % The number of numerator and denominator parameters and number of samples for the 
    % delay for the discrete-time model (in this ZOH input case) 
    % are [nb nf nk]=[4 4 1]
    nnd=[4 4 1];
    Moe=oe(data,nnd);
    Marx=arx(data,nnd);
    Mn4sid=n4sid(data,4);

    % Note that the SID IV4 routine delivers an error message on the same
    % data. Uncomment the line below if you want to test
    %Miv4=iv4(data,nnd);
   

    % Bode plot responses of the different identified models
    wmin=log10(0.1);
    wmax=log10(2*pi/Ts/2);
    w=logspace(wmin,wmax,200);
    
    [mag_0,phase_0]=bode(M0,w);
    [mag_srivc,phase_srivc]=bode(Msrivc,w);
    
    [mag_arx,phase_arx]=bode(Marx,w);    
    [mag_oe,phase_oe]=bode(Moe,w);
    [mag_n4sid,phase_n4sid]=bode(Mn4sid,w);
    %[mag_iv4,phase_iv4]=bode(Miv4,w);

    fontsize=18;      % Font size in the different plots
    linethickness=1.5;  % linethickness des lignes pour les tracés des courbes vrais
    axis_amp=[w(1) w(length(w)) -70 50];
    axis_phase=[w(1) w(length(w)) -700 200];

    figure,
    subplot(2,1,1)
    semilogx(w(1:4:end),20*log10(mag_arx(1,1:4:end)),'b-*');%,'linewidth',linethickness,'MarkerSize',5);
    hold on
    semilogx(w(1:4:end),20*log10(mag_oe(1,1:4:end)),'b-p');%,'linewidth',linethickness,'MarkerSize',5);
    semilogx(w(1:4:end),20*log10(mag_n4sid(1,1:4:end)),'b-o');%,'linewidth',linethickness,'MarkerSize',5);

    semilogx(w(1:4:end),20*log10(mag_0(1,1:4:end)),'r-');%,'linewidth',linethickness)
    %semilogx(w(1:4:end),20*log10(mag_iv4(1,1:4:end)),'b-*');%,'linewidth',linethickness,'MarkerSize',5);

    legend('DT ARX model','DT OE model','DT N4SID model','Actual system')

    axis(axis_amp)
    ylabel('Magnitude (dB)')
    grid
    set(gca,'FontSize',fontsize)
    set(findall(gcf,'type','text'),'FontSize',fontsize)
    hold off

    subplot(2,1,2)
    semilogx(w(1:4:end),20*log10(mag_srivc(1,1:4:end)),'b-s');%,'linewidth',linethickness,'MarkerSize',5);
    hold on
    semilogx(w(1:4:end),20*log10(mag_0(1,1:4:end)),'r-');%,'linewidth',linethickness)
    ylabel('Magnitude (dB)')
    legend('CT SRIVC model','Actual system')
    xlabel('Frequency (rad/s)')
    grid
    axis(axis_amp)
    hold off

    set(gca,'FontSize',fontsize)
    set(findall(gcf,'type','text'),'FontSize',fontsize)

    echo on
    %   Note that OE can deliver good results sometime.
    %   Run this program several times to get an idea of the average results obtained with OE
    %   This example illustrates some of the well-known difficulties that may appear in DT modelling
    %   - sensitivity to the initialization
    %   - numerical issues in the case of fast sampling
    %   - a priori knowledge of the system relative degree not easy to accommodate
    %   - non inherent data prefiltering
    %   These difficulties require extra care in DT modelling while the direct  
    %   SRIVC approach is free from them !     
    echo off