%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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()