era.m

来自「本文提出的特征系统实现算法是多输入多输出的时域模态参数识别方法,」· M 代码 · 共 56 行

%共享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是采样频率