from distribute modal identification algorithm based on NExT and ERA techniques. by Lei
distribute modal identification algorithm based on NExT and ERA techniques.

era.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Eigensystem Realization Algorithm using Matlab 
%
%    input:
%      y(i,(k-1)*m+j)   IRF   Yij(k) timeK input j output i
%        ncols     no. of columns in the hankel matrix
%                  (this should be twice the no. of modes to estimate
%                  don't forget to include some computational modes!)
%        nrows     no. of rows in the hankel matrix
%     fs     sampling frequency in Hz
%     inputs,m  no of inputs;

function [fd1,zeta1,shapes,Lambda1,partfac,EMAC1,sv,A,B,C]=era(Y,fs,ncols,nrows,inputs,varargin)
% Y is in form of Y(ouputs,inputs,data sequence).
% ncols = bcols * inputs.
% nrows = brows * outputs.  ----> brows == number of rows of C?

% Lambda1 and partfac are in continuous-time.

error(nargchk(5,7,nargin));
if nargin==5
    Na = []; % system order
    Msg = 0; % want to display msg?
   % MSG = NARGCHK(LOW,HIGH,N) returns an appropriate error message string if
   % N is not between LOW and HIGH. If it is, NARGCHK returns an empty matrix. 
 
elseif nargin==6
    Na = varargin{1};
    Msg = 0;
elseif nargin ==7
    Na = varargin{1};
    Msg = varargin{2};
end

[outputs,npts]=size(Y);
if outputs > npts % The outputs should be smaller than the data sequence.
    Y=Y';
    [outputs,npts]=size(Y);
end
npts=npts/inputs; % Now 'npts' represents the total data points in Y.
%
%     find block sizes
%
brows=fix(nrows/outputs); % brows = how many output blocks.
nrows=outputs*brows; % Redefine the row numbers.
bcols=fix(ncols/inputs); % bcols = how many time steps.
ncols=inputs*bcols; % Redefine the column numbers.
m=inputs; % inputs number.
q=outputs; % outputs number.
%
%     check no. of data points
%
nused=bcols+brows-1; % Total data points used. Eq.5.44.
if npts < bcols+brows-1 % Is nps greater than the total data points needed.
    fprintf(' \n Not enough Data Points for NROWS NCOLS requested\n')
    return
end
if Msg
    fprintf(' \n %4i time samples used: %7.3f secs of data\n',nused,nused/fs)
end
%
%     Form the Hankel matrix H(0).
%
H0=zeros(nrows,ncols);
for ii=1:brows
    H0((ii-1)*q+1:ii*q,1:bcols*m)=Y(:,(ii-1)*m+1:(bcols+ii-1)*m); % Eq.5.44.
end
%
%   Decompose the data matrix
%
[R1,Sigma1,S1]=svd(H0,0);
sv=diag(Sigma1);

%      Plot the Singular Values and prompt the user for cutoff
%

if isempty(Na)
    subplot(211)
    semilogy(sv)                  % plot singular values
    xlabel('Singular Value')
    ylabel('log(sv)')
    grid
    subplot(212)
    semilogy(sv(1:1/5*length(sv)),'o-')                  % plot singular values
    xlabel('Singular Value')
    ylabel('log(sv)')
    grid
    cut=input('Estimate the Singular Value Cutoff: ');
else
    cut = Na;
end
%cut=8;

%     Truncate the matrices using the cutoff

D=diag(sqrt(sv(1:cut)));      % build square root of the singular values.
Dinv=inv(D);						% (sigma)^(-1/2) ---- Eq.5.34
Rn=R1(:,1:cut);                 % use only the principal eigenvectors Eq.(5.31)
Sn=S1(:,1:cut);                 %   "    "    "    "     "   " Eq.(5.31).
%
%     Build the second Hankel matrix H(1).
%
H1=zeros(nrows,ncols);
for ii=1:brows
    H1((ii-1)*q+1:ii*q,1:bcols*m)=Y(:,ii*m+1:(bcols+ii)*m); % Eq.5.27 when k=1.
end
%
%    Calculate the realization of A and find eigenvalues and eigenvectors
%
A=Dinv*Rn'*H1*Sn*Dinv ;        % build A as per ERA. Eq.(5.34)
B=D*Sn';B=B(:,1:m) ;           % build B as per ERA. Eq.(5.34)
C=Rn*D; C=C(1:q,:);
% build C as per ERA. Eq.(5.34)
clear H0 H1;
[Vectors,Values]=eig(A);
%
%    Extract the freq and damping factor information from eigenvalues
%
Lambda=diag(Values);           % roots in the Z-plane
s=log(Lambda).*fs;           % Laplace roots Page.132 line 9th.
zeta=-real(s)./abs(s)*100;        % damping factors: ksi*100
fd=imag(s)./2./pi;  % damped natural freqs:
%
%  Calculate Mode Shape, Modal Participation factors
%  Modal observ. and Modal controllability for EMAC calculation
%
shapes=C*Vectors; % Mode shapes.

InvV=inv(Vectors);
partfac=InvV*D*Sn(1:m,:)';

%%%by Han
if isempty(find(imag(diag(Values))==0,1))
    Ac = fs*logm(Values);  % continuous-time A
    partfac = inv(Values-eye(size(Values)))*Ac*partfac;
else
    partfac = zeros(size(partfac));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try to calculate the parameters EMAC using what the book does.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ahat=Values; % or Ahat=inv(vectors)*A*vectors; Eq(5.37).
Bhat=InvV*B; % Eq(5.37)
Chat=C*Vectors; % Eq(5.37)

% Eq(5.42) in the book.
Lambda=diag(Values);
qhati=zeros(1,ncols);qhat=zeros(cut,ncols); % Eq(5.42).

for ii=1:cut
    for jj=1:bcols % bols is the total input block in the row of Hankel matrix.
        qhati(1,(jj-1)*m+1:jj*m)=Bhat(ii,:)*Lambda(ii)^(jj-1); % Eq(5.42).
    end
    qhat(ii,:)=qhati;
end

Qbar=InvV*D*Sn'; % Eq.(5.45): Qbar and Eq.(5.47).

for ii=1:cut
    qbar(ii,:)=Qbar(ii,:);
end

EMAC=[];
for ii=1:cut
    qhati=qhat(ii,:); % Eq(5.47).
    qbari=qbar(ii,:);
    EMAC(ii,1)=abs(qbari*qhati')/sqrt(abs(qbari*qbari')*abs(qhati*qhati'))*100; % Eq. 5.49.
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mode Singular Value (MSV): not successful.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    Sort into ascending order
%
[fd,I]=sort(fd);
zeta=zeta(I);
shapes=shapes(:,I);
s=s(I); % Laplace roots Page.132 line 9th
partfac=partfac(I,:);
EMAC=EMAC(I);

%      Remove the negative freqs
%
lower=1;
upper=cut;
for ii=1:cut
    if fd(ii) <= 0
        lower=ii+1;
    end
    if fd(cut-ii+1) >= 0.499*fs
        upper=cut-ii;
    end
end

%fd=fd(lower:upper);
%zeta=zeta(lower:upper);

fd1=fd(lower:upper); % It is the damped natural frequency.
zeta1=zeta(lower:upper);
fd1=fd1./sqrt(1-(zeta1/100).^2)% Calculate the undamped natural frequency: wn=wd/sqrt(1-ksi^2);
Lambda1=s(lower:upper);
shapes=shapes(:,lower:upper);
partfac=partfac(lower:upper,:);
EMAC1=EMAC(lower:upper);

[fd1,ii]=sort(fd1); 
zeta1=zeta1(ii)
EMAC1=EMAC1(ii);
shapes=shapes(:,ii); 
partfac=partfac(ii,:);
Lambda1=Lambda1(ii);
%%%%%%%%%%%%%%%%Ͻģ̬ͱΪʵBenchmark
[m,n]=size(shapes);
for i=1:m
    for j=1:n
shapes(i,j)=sign(angle(shapes(i,j)))*abs(shapes(i,j));
    end
    AA=max(abs(shapes(:,j)));
    shapes(:,j)=shapes(:,j)/AA;
end
shapes=shapes;
% BB=zeros(1,n);
% shapes=[BB;shapes];
% save shapes.mat shapes
% for i=1:9
%     figure
%     plot(shapes(:,i))
% %      figure(gcf);
% end
% Fl=linspace(0,m,m+1);
%    for i=1:m
%        figure(i)
%    plot(shapes(:,i),Fl,'b-','linewidth',2.0)
%    axis([-1,1,0,4])
%    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Msg
    fprintf('\n       fd1     zeta1     EMAC1 \n');
    disp([fd1 zeta1 EMAC1])
end
% for i=1:9
%     figure
%     plot()

Contact us