function [fd,zeta,shapes,partfac,EMAC,sv,A,B,C]=s (Y,fs,ncols,nrows,inputs)
% [fd,zeta,shapes,partfac,EMAC,sv]=era(Y,fs,ncols,nrows,inputs)
%
% Eigensystem Realization Algorithm using Matlab for a MIMO case!
% written by JJ Hollkamp June 91
%
% 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;
%
shift=10; % Time shift for EMAC calc
[outputs,npts]=size(Y);
if outputs > npts
Y=Y';
[outputs,npts]=size(Y);
end
npts=npts/inputs;
%
% find block sizes
%
brows=fix(nrows/outputs);
nrows=outputs*brows;
bcols=fix(ncols/inputs);
ncols=inputs*bcols;
m=inputs;
q=outputs;
%
% check no. of data points
%
nused=bcols+brows-2+2*shift;
if npts < bcols+brows-2+shift+shift
fprintf(' \n Not enough Data Points for NROWS NCOLS requested\n')
return
end
fprintf(' \n %4i time samples used: %7.3f secs of data\n',nused,nused/fs)
%
% Form the Hankel matrix (shift last block column)
%
H0=zeros(nrows,ncols);
for ii=1:brows-1
H0((ii-1)*q+1:ii*q,1:(bcols-1)*m)=Y(:,(ii-1)*m+1:(bcols+ii-2)*m);
H0((ii-1)*q+1:ii*q,(bcols-1)*m+1:ncols)=Y(:,(bcols+ii-3+shift)*m+1:(bcols+ii-2+shift)*m);
end
%
% Shift last block row for EMAC calculation
%
H0((brows-1)*q+1:nrows,1:(bcols-1)*m)=Y(:,(brows-2+shift)*m+1:(brows+bcols-3+shift)*m);
H0((brows-1)*q+1:nrows,(bcols-1)*m+1:ncols)=Y(:,(brows+bcols-4+2*shift)*m+1:(brows+bcols-3+2*shift)*m);
%
% Decompose the data matrix
%
[P,S,Q]=svd(H0,0);
sv=diag(S);
%
% Plot the Singular Values and prompt the user for cutoff
%
subplot(211)
semilogy(sv) % plot singular values
xlabel('Singular Value')
ylabel('log(sv)')
grid
subplot(212)
semilogy(sv(1:30)) % plot singular values
xlabel('Singular Value')
ylabel('log(sv)')
grid
cut=input('Estimate the Singular Value Cutoff: ');
%
% Truncate the matrices using the cutoff
%
D=diag(sqrt(sv(1:cut))); % build sqrt(S)
Dinv=inv(D);
P=P(:,1:cut); % use only the principal eigenvectors
Q=Q(:,1:cut); % " " " " " "
%
% Build the second Hankel matrix
%
H1=zeros(nrows,ncols);
for ii=1:brows-1
H1((ii-1)*q+1:ii*q,1:(bcols-1)*m)=Y(:,ii*m+1:(bcols+ii-1)*m);
H1((ii-1)*q+1:ii*q,(bcols-1)*m+1:ncols)=Y(:,(bcols+ii-2+shift)*m+1:(bcols+ii-1+shift)*m);
end
%
% Shift last block row for EMAC calculation
%
H1((brows-1)*q+1:nrows,1:(bcols-1)*m)=Y(:,(brows-1+shift)*m+1:(brows+bcols-2+shift)*m);
H1((brows-1)*q+1:nrows,(bcols-1)*m+1:ncols)=Y(:,(brows+bcols-3+2*shift)*m+1:(brows+bcols-2+2*shift)*m);
%
% Calculate the realization of A and find eigenvalues and eigenvectors
%
A=Dinv*P'*H1*Q*Dinv; % build A as per ERA
B=D*Q';B=B(:,1:m); % build B as per ERA
C=P*D; C=C(1:q,:); % build C as per ERA
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
fd=imag(s)/2./pi; % damped natural freqs
zeta=-real(s) ./abs(s)*100; % damping factors
%
% Calculate Mode Shape, Modal Participation factors
% Modal observ. and Modal controllability for EMAC calculation
%
shapes=P(1:q,:)*D*Vectors;
InvV=inv(Vectors);
partfac=InvV*D*Q(1:m,:)';
%
% Calculate the EXTENDED modal amplitude coherence
%
mom=P((brows-1)*q+1:nrows,:)*D*Vectors;
mcm=InvV*D*Q((bcols-1)*m+1:ncols,:)';
%clear P D Q InvV;
for jh=1:cut
%
% observability EMAC
%
EMAC1=abs(shapes(:,jh)*Lambda(jh)^(brows-2+shift))./abs(mom(:,jh));
for ii=1:q
if EMAC1(ii) > 1.
EMAC1(ii) =1/EMAC1(ii);
end
end
%
% controllability EMAC
%
EMAC2=abs(partfac(jh,:)*Lambda(jh)^(bcols-2+shift))./abs(mcm(jh,:));
for ii=1:m
if EMAC2(ii) > 1.
EMAC2(ii) =1/EMAC2(ii);
end
end
EMAC(jh)=100*sum(EMAC1)/q*sum(EMAC2)/m; % Average EMAC
end
EMAC=EMAC';
%
% Sort into ascending order
%
[fd,I]=sort(fd);
zeta=zeta(I);
shapes=shapes(:,I);
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);
shapes=shapes(:,lower:upper);
partfac=partfac(lower:upper,:);
EMAC=EMAC(lower:upper);
[fd zeta EMAC];
