lsce.m

% lsce.m % % This is a script that demonstrates the least squares complex % exponential (LSCE) modal parameter estimation algorithm. % This example expects an impulse response function to be in % the memory variable 'irf' when called as well as the time % increment in a variable called 'deltaT'. 30 roots % (15 complex conjugate pairs) using an overdetermination factor % of 4 will be found by this script in the default condition. clear; load frf MaxRoots=input('Maximum number of roots to try: (30)'); if isempty(MaxRoots),MaxRoots=30;end; OverDet=input('Overdetermination factor: (4)'); if isempty(OverDet),OverDet=4;end; % Choose the solution method disp('The following solution methods are available:') disp(' 1) Normal Equations') disp(' 2) Psuedo-Inverse') disp(' 3) Singular Value Decomposition') SolMeth=input('Choose the solution method: (1)') if isempty(SolMeth),SolMeth=1;end; if (SolMeth<1),SolMeth=1;end; if (SolMeth>3),SolMeth=1;end; Nrows=size(irf,1); Ncols=size(irf,2); LAMBDA=zeros(MaxRoots,MaxRoots); % Loop for increasing number of roots for Nroots=2:MaxRoots; clear data, clear A, clear B, clear x, clear z_lambda clear s_lambda, clear est_lambda Nroots Nroots1=Nroots+1; Neqn=OverDet*Nroots; % Calculate the number of equations that can be found from the % irf data Maxeqn=0.80*(Ncols-(Nroots1)); if(Neqn>Maxeqn), break; end; for ii=1:Neqn; data(ii,:)=irf(ii:ii+Nroots); end; % Formulate normal equations solution if (SolMeth==1), cov_data=data'*data; % Set up Ax=B form A(1:Nroots,1:Nroots)=cov_data(1:Nroots,1:Nroots); B(1:Nroots,1)=-cov_data(1:Nroots,Nroots+1); % Solve for coefficients InvCond(Nroots)=1/cond(A); x=inv(A)*B; end; % Formulate psuedo-inverse solution if (SolMeth==2), A(:,1:Nroots)=data(:,1:Nroots); B(:,1)=-data(:,Nroots1); InvCond(Nroots)=1/cond(A); x=pinv(A)*B; end; % Formulate singular value decomposition solution if (SolMeth==3), break end; % Rearrange coefficients for root solving call x(Nroots1)=1; alpha=flipud(x); alpha(1)=1; % Solve for roots in Z domain z_lambda=roots(alpha); s_lambda=log(z_lambda)./deltaT; [Y,I]=sort(imag(s_lambda)); est_lambda=s_lambda(I); % Save estimated modal frequencies for ii=1:Nroots, LAMBDA(Nroots,ii)=est_lambda(ii); end; end; % Plot error chart (inverse of condition number) fig1=figure(1); clg semilogy(InvCond) xlabel('Model Order') ylabel('Error Chart (Inverse of Condition Number)') title('Model Order Determination') % Plot error chart (Normalized singular values) fig2=figure(2); clg [U,S,V]=svd(A,0); SingVals=diag(S); SingVals=SingVals./SingVals(1); semilogy(SingVals) xlabel('Model Order') ylabel('Error Chart (Normalized Singular Values)') title('Model Order Determination') % Plot frequency stability diagram fig3=figure(3); clg for ii=1:MaxRoots, clear yy; yy=ii*ones(1,MaxRoots); plot(imag(LAMBDA(ii,:))/(2*pi),yy,'+') xlabel('Frequency (Hz)') ylabel('Model Order') title('Frequency Stability Diagram') axis([0,50,0,MaxRoots]) hold on end POWER=FRF(1,:).*conj(FRF(1,:)) PowerMax=max(POWER); POWER=POWER/PowerMax; plot(freq,POWER(1:401)*MaxRoots)