%共享ERA/DC算法和SSI算法程序
function [fr,damp,v,d0]=eradc(n1,n2,h,d);
%  fr:frequecy matrix
%  damp:damp ratio matrix
%  n2:the order of hankel matrix

% 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));
h1=hankel(r(2:n2+1),r(n2+1:2*n2));

% ERA start
[p0 d0 q0]=svd(h0);

%judge value's validlity
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);
qr=q0(:,1:i);
dr=dd^(-0.5);

a=dr*pr'*h1*qr*dr;
[v,z]=eig(a);

z1=diag(z);
z2=log(z1);
fr=abs(z2)*d;
damp=real(z2)./fr*d;



%n1是自相關函數的平均次數
%21是hankel矩陣的階數
%h是脈沖的時程響應
%d是采樣頻率