load cttfdata % contains m0{k}, k=1,..,10 (earlier randomly generated)
               % contains also the standard input seguence  ud =[u1 u2], 
            % u1 RBS, ur1 RBS (generated from 
            % [u1 ur1] = idinput([1000 2},'RBS',[0.1 0.5]);) 
            % u2 = u1+ur1; makes u1 and u2 50% correlated and u2 having
            % levels -2,0,2

%smpl % one of 'normal, 'fast', 'slow
mdef = idpoly(1,1,1,1,1,'ts',0); % just a default model if some estimation
                                 % fails
snr = 10; % (amplitude) signal to noise ratio
mnumbers = [1:10]; %Model numbers to be tested
MCn = 10; %Number of Monte Carlo runs
 
for smpl ={'normal','fast','slow'}

for k=mnumbers % over the ten models
    bw1=bandwidth(tf(m0{k}(1,1),'m'));
    bw2=bandwidth(tf(m0{k}(1,2),'m'));
    ts0 = 2*pi/(max(bw1,bw2)*10);
    
    switch smpl{1}
        case 'normal'
        ts = ts0;
        case 'fast'
            ts = ts0/10;
        case 'slow'
            ts = ts0*10;
    end
     u =iddata([],ud,ts);
    z0{k} = [sim(m0{k},u),u];
    zz = z0{k};
    y = z0{k}.y;
    ylev=y'*y/length(y);
    for monte = 1:MCn % over the realizations
        [k,monte]
        e = randn(1000,1);
        yn = y+1/snr*sqrt(ylev)*e;
        zz.y = yn;
        try
        m1 = srivc(zz,[m0{k}.nb m0{k}.nf  0 0]); %Contsid toolbox
        catch
            m1=mdef; % just not to error out if any difficulties.
        end
        try
        m2 = coe(zz,[m0{k}.nb m0{k}.nf  0 0],1/zz.ts/2); %Contsid toolbox
        catch
            m2 = mdef;
        end
        try
        m3 = tfest(zz,[m0{k}.nb m0{k}.nf]); % Special code provided
        catch
            m3 = mdef;
        end
        try
        m4 = tfest(zz,m1);
        catch
            m4 = mdef;
        end
 
        [dum,tfc]= compare(z0{k},m1,m2,m3,m4); %SITB
        tfc(find(tfc<-100))   = -100;
        tfc(find(isnan(tfc))) = -100;
        fit(k,monte,:)=tfc(:);
        fit
    end
end
table = squeeze(median(fit,2));
switch smpl{1}
    case 'normal'
        table_normal = table;
    case 'fast'
        table_fast = table;
    case 'slow'
        table_slow = table;
end
end
    

