fprintf('\n----------------------------------------------------------------------')

for st=40 
    save st st
clear;clf;format short g
load st
Method = 'RGELS';
FF=1;    %  Forgetting factor
sigma=0.1;  %  sigma=0.1, 1.0  % variance of noise
figure(1);
length1=3100;
na = 1; nb =1; nd=1; nc=1;dd=1; %  dd: delay =1
a=[-0.95]; b=[1.52]; d=[-1.0];c=[0.72];
%a=[1, 0.412, 0.309];  b=[0, 0.6804,0.6303];
par0=[a,b,c,d]';  roots([1,a]);
k0=16; %k0=16 
n=length(par0); PP = eye(n)*1e6;   RR = 1;  h = ones(n,1)*1e-6;  par1 = ones(n,1)*1e-6;
a1=[1,a]; b1=[0,b];d1=[1,d];c1=[1,c];
sy=f_integral(a1,b1);      sv=f_integral(c1,d1);
delta_ns = sqrt(sv/sy)*100*sigma;
    ns = [sv,sy,delta_ns]
    fprintf('$\\sigma_v^2=%6.2f^2$, $\\delta_{\\ns}=%6.2f%s \n',sigma,delta_ns,'\%$');  
rand('state',st);   
randn('state',1);
eta=f_rn(length1);  
u=eta(:,1); v=eta(:,2);
Gz=tf(b1,a1,1); Gn=tf(d1,c1,1);e=eta(:,2);
y=lsim(Gz,u) -lsim(Gn,v * sigma);
jj = 0; j1 = 0; j2 = 0;
for k = 20:length1   
    jj=jj+1;
    h = [-y(k-1:-1:k-na); u(k-1:-1:k-nb);-e(k-1:-1:k-nc);v(k-1:-1:k-nd)];   %delay = 1
    yy=y(k);
    [par1,delta,PP,RR] = f_ls(par0,par1,h,yy,FF,'LS',PP,RR);
    e(k)=y(k)+par1(1:na)'*y(k-1:-1:k-na) - par1(na+1:na+nb)'*u(k-1:-1:k-nb);
    v(k)=e(k)+par1(na+nb+1:na+nb+nc)'*e(k-1:-1:k-nc)-par1(na+nb+nc+1:n)'*v(k-1:-1:k-nd);
    ls(jj,:)=[jj, par1',delta];
    if (jj==100)|(jj==200)|(jj==300)|mod(jj,500)==0
        j1 = j1+1;
        ls_100(j1,:)=[jj, par1', delta*100];
    end
    ls_100(j1+1,:)=[ 0, par0', 0];
end
fprintf('\nst = %d,  Method = %s, ($\\sigma_v^2=%5.2f^2$, $\\delta_{\\ns}=%6.2f\\%s$)', st, Method, sigma,delta_ns,'%');  
fprintf('\n %s  \n','$k$  &  $a$  &  $b$  &  $c$  &  $d$  & $\delta\ (\%)$ \\');  
fprintf('%5d &%10.5f &%10.5f &%10.5f &%10.5f &%10.5f & \\\\\n',ls_100');
jk=(17:10:2990)';
plot(ls(jk,1), ls(jk,n+2));
grid;
end
fprintf('=====================================================================\n')
if sigma==0.1
    z1=[ls(:,1), ls(:,n+2)];
    save z1 z1
else 
    load z1
    z0=[z1, ls(:,n+2)];
    figure(2)
    plot(z0(jk,1),z0(jk,2),'k',z0(jk,1),z0(jk,3),'b')
    axis([0, 3000, 0, 0.33])
    text(600,0.058,'{\it\sigma_v^2} = 1.00^2')
    text(600,0.13,'{\it\sigma_v^2} = 0.10^2')
    line([247,620],[0.024,0.119])    
    end
xlabel('\it        t'); ylabel('{\it    \delta}');