MATLAB and C Programs - Digital Signal Processing: A Practitioner's Approach - page 164

yk=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:1000;
Dz=[1 -a1 -a2];
R=roots(Dz);
rho(1,k)=abs(R(1)); theta(1,k)=angle(R(1));
rho(2,k)=abs(R(2)); theta(2,k)=angle(R(2));
alpha=0.98*alpha1;
yk1=yk;
yk=z(k);
xk2=xk1; xk1=xk;
xk=a1*xk1+a2*xk2;
sk3=sk2; sk2=sk1; sk1=sk;
sk=a1*sk1+a2*sk2+yk;
ek=yk-xk;
Jk=ek^2-v;
Obj(k)=Jk;
gk=[2*ek*sk1 2*ek*sk2 1]';
Hk_new=Hk/alpha;
delH_num=Hk_new*gk*gk'*Hk_new;
delH_den=1+gk'*Hk_new*gk;
delH=delH_num/delH_den;
Hk=Hk_new-delH;
old_delp_k=delp_k;
P(:,k)=pk;
pk=pk+delp_k;
a1=pk(1); a2=pk(2); v=pk(3);
Jk_hat=Jk-gk'*old_delp_k;
delp_k=Hk*gk*Jk_hat;
del_a1=delp_k(1);
del_a2=delp_k(2);
delx_k=sk1*del_a1+sk2*del_a2;
xk_corkt=xk+delx_k;
xk=xk_corkt;
alpha1=alpha1*0.9995+0.0005;
x(k)=xk;
y(k)=yk;
T(k)=k;
%%%
A(k)=a1; B(k)=a2; C(k)=v;
end
subplot(221);plot(T,0.5*x,T,y);grid;zoom;
subplot(222);plot(T,A,T,B,T,C);grid
subplot(212);polar (theta,rho,'o');grid;