function [fr,damp,v,d0]=ssi2(n1,n2,h,d);
% generate R-matrix from h
[r,c]=size(h);
r(1:c-n1)=0;
for q=1:c-n1;
r(q)=0;
for i=1:n1;
r(q)=h(q+i)*h(i)+r(q);
end;
end;
% generate Hankel-matrix from R-matrix
h0=hankel(r(1:n2),r(n2:2*n2-1));
%compute xi
yp=h0(1:n2/2,:);
yf=h0(n2/2+1:n2,:);
pref=yf*yp'*pinv(yp*yp')*yp;
[r,c]=size(pref);
pi=pref(1:r-1,:);
[p0 d0 q0]=svd(pi);
%judge qiyi value
i=4;
cr1=3;
cr2=3;
cr3=3;
while (cr1>1.5 | cr2>1.5 | cr3>1.5) & i<20
cr1=d0(i,i)/d0(i+2,i+2);
cr2=d0(i+2,i+2)/d0(i+4,i+4);
cr3=d0(i+4,i+4)/d0(i+6,i+6);
i=i+2;
end;
i=i-4;
pr=p0(:,1:i);
dd=d0(1:i,1:i);
dr=dd^0.5;
qr=q0(:,1:i);
oi=pr*dr;
xi=pinv(oi)*pi;
%compuet xi1
clear yf p0 d0 q0 pr dd dr qr
pi=pref(2:r,:);
[p0 d0 q0]=svd(pi);
xi1=pinv(oi)*pi;
a=xi1*pinv(xi);
[v,z]=eig(a);
z1=diag(z);
z2=log(z1);
fr=abs(z2)*d;
damp=real(z2)./fr*d; |