v3_lsce.m

% v3_lsce.m % % This is a script that demonstrates the least squares complex % exponential (LSCE) modal parameter estimation algorithm. % A maximum of 30 roots % (15 complex conjugate pairs) using an overdetermination factor % of 4 will be found by this script in the default condition. %********************************************************************** % Author: Randall J. Allemang % Date: 08-Apr-97 % Structural Dynamics Research Lab % University of Cincinnati % Cincinnati, Ohio 45221-0072 % TEL: 513-556-2725 % FAX: 513-556-3390 % E-MAIL: randy.allemang@uc.edu %********************************************************************* % clear; clg; close all; % % Generate an impulse response function from a known MKC model % pi=3.14159265; Ndof=3; BS=511; counter=6; % 3 dof system mass=[8,0,0;0,10,0;0,0,12]; stiff=10000.*[26,-12,0;-12,22,-10;0,-10,10]; damp=[13,-6,0;-6,11,-5;0,-5,5]; null=[0,0,0;0,0,0;0,0,0]; % Form 2N x 2N state space equation. a=[null,mass;mass,damp]; b=[-mass,null;null,stiff]; [x,d]=eig(b,-a); % Sort Modal Frequencies orig_lambda=diag(d); [Y,I]=sort(imag(orig_lambda)); lambda=orig_lambda(I); lambda xxx=input('Hit any key to continue'); clear a,clear b,clear x,clear d; % Set up DSP parameters max_omega=1.5*max(abs(lambda)); max_freq=max_omega/(2*pi); samp_freq=2*max_freq; deltaT=1/samp_freq; omega=linspace(0,max_omega,BS); freq=omega/(2*pi); % Inversion loop for ii=1:BS; om=omega(ii); B=-mass.*om.*om + damp.*j.*om + stiff; H=inv(B); frf(:,ii)=reshape(H,Ndof*Ndof,1); end subplot(211),semilogy(freq,abs(frf(counter,:))) ylabel('Magnitude'),grid scale=360.0/(2.0*pi); subplot(212),plot(freq,scale.*angle(frf(counter,:))) xlabel('Frequency (Hertz)'),ylabel('Phase'),grid pause POWER=zeros(size(frf(1,:))); for ii=1:9 POWER=POWER+frf(ii,:).*conj(frf(ii,:)); end % Formulate two-sided frf in MatLab required format FRF=[frf(counter,:),0,conj(fliplr(frf(counter,2:BS)))]; irf=ifft(FRF); % Remove small imaginary part of irj due to numerical error irf=irf+conj(irf); clg; plot(real(irf)); pause % Add 1 percent random noise to impulse response function max_irf=max(real(irf)); irf=irf+0.01*max_irf*(rand(1,2*BS)-0.5); plot(real(irf)); pause % % Start Least Squares Complex Exponential Algorithm % 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') pause % 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') pause % 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,max_freq,0,MaxRoots]) hold on end PowerMax=max(POWER); POWER=log10(POWER/PowerMax); PowerMin=min(POWER); POWER=POWER-PowerMin; PowerMax=max(POWER); POWER=POWER/PowerMax; plot(freq,POWER*MaxRoots)