%------------------------------------------
% Plotting Global Mode Shape
%------------------------------------------
function [frequency,modeshape,damping]=fdd_plot(freq,df,U,SV)
frequency=[];
modeshape=[];
damping=[];
%
h1=figure;h2=figure;
bedw=5;topw=75;sscrn=get(0,'ScreenSize'); % bed width, top width
set(0,'Units','pixels');
set(h1,'position',[bedw,floor(sscrn(4)/2),sscrn(3)/2-2*bedw,floor(sscrn(4)/2)-(topw)]);
set(h2,'position',[bedw,topw-2*bedw,sscrn(3)/2-2*bedw,ceil(sscrn(4)/2)-(2*topw-bedw)]);
%
iter=1;
while(iter>=1),iter=iter+1;
%
% Singular value plot
%
figure(h1);plot(freq,SV);
set(gca,'yscale','log');
xlabel('Frequency (Hz)');ylabel('Singular Value');
title('Result for Singular Value Decomposition');
%
BUTTON=ones(2,1);
while(BUTTON==ones(2,1))
% Pick 2 points, return 1(left), 2(middle), 3(right)
[x_coord(1),y,BUTTON(1)]=ginput(1);
if ((BUTTON(1)~=1)), break;
end
[x_coord(2),y,BUTTON(2)]=ginput(1);
if ((BUTTON(2)~=1)), break;
end
% Zoom Result of SVD
index_mod=ceil(x_coord/df); % Convert Freq. to index No.
low=index_mod(1);high=index_mod(2);
%
% Pick-Picking method
%
[powermax_PP,idx_tmp_PP]=max(SV(low:high));
idx_f_PP=idx_tmp_PP+low-1;
f_PP=freq(idx_f_PP);
% Plot the extracted frequency
plot(freq(low:high),SV(low:high)); set(gca,'yscale','log');
xlabel('Frequency (Hz)');ylabel('Singular Values');
title('Result of Singular Value Decomposition');
%
% Curve Fitting Method
%
fit_ndt=5;
fit_order=4;
fit_i=(idx_tmp_PP-fit_ndt):(idx_tmp_PP+fit_ndt);
fit_idx=fit_i+low-1;
fit_x=freq(fit_idx);
fit_y=sqrt(SV(fit_idx).^(-1));
fit_p=polyfit(fit_x,fit_y,fit_order);
app_y=polyval(fit_p,fit_x);
app_SV=app_y.^(-2);
% Check fitting
figure(h2);
plot(fit_x,SV(fit_idx),'o-');
hold on;
plot(fit_x,app_SV,'g+-');
legend(' Raw ',' Fitted',1);
% Fitted Natural Frequency
[powermax_CF,idx_tmp_CF]=max(app_SV);
idx_f_CF=(idx_tmp_PP+low-1)-(fit_ndt-idx_tmp_CF+1);
f_CF=freq(idx_f_CF);
% Plot
plot(f_PP,SV(idx_f_PP),'rx');
plot(freq(fit_idx),powermax_PP/2,'rx-');
plot(f_CF,app_SV(idx_tmp_CF),'g+');
plot(freq(fit_idx),powermax_CF/2,'g+-');
% Damping estimation from the half-power bandwidth
fit_frq=roots([fit_p(1:fit_order)*(fit_p(fit_order+1)-sqrt(1/(powermax_CF/2)))]);
%
nreal=0;
for ifit_order=1:fit_order
nreal=nreal+isreal(fit_frq(ifit_order));
end
frq_beta_a=0; frq_beta_b=0;
if (nreal==fit_order)
frq_beta_a=fit_frq(floor(fit_order/2)+1);
frq_beta_b=fit_frq(floor(fit_order/2));
elseif (nreal<fit_order)
for ifit_order=1:fit_order
if((isreal(fit_frq(ifit_order))==1)&&(frq_beta_b==0))
frq_beta_b=fit_frq(ifit_order);
elseif ((isreal(fit_frq(ifit_order))==1)&&(frq_beta_a==0))
frq_beta_a=fit_frq(ifit_order);
end
end
end
%
damp=(frq_beta_b-frq_beta_a)/(2*f_CF)*100;
legend(' Pick Picking',' Curve Fitting',1);
% Check estimation
app_y_beta_a=polyval(fit_p,frq_beta_a);
app_SV_beta_a=app_y_beta_a.^(-2);
app_y_beta_b=polyval(fit_p,frq_beta_b);
app_SV_beta_b=app_y_beta_b.^(-2);
plot(frq_beta_a,app_SV_beta_a,'md',frq_beta_b,app_SV_beta_b,'md',f_CF,app_SV(idx_tmp_CF),'md');
xlabel('Frequency (Hz)');ylabel('Singular Values');
text_out=['Freq=',num2str(f_CF),'Hz;','Damping=',num2str(damp),'%','1/2-P-Width=',num2str(frq_beta_b-frq_beta_a),'Hz'];
title(text_out);hold off;
end
%
% Mode Shape plot
%
mode=sign(imag(U(:,idx_f_CF))).*abs(imag(U(:,idx_f_CF)));
figure(h2);
plot([mode(1:end)],[3:-1:0],'-o');
title( text_out );
% Another Selection ?
[x_coord(1),y,BUTTON(1)]=ginput(1);
if (BUTTON(1)~=1)
dummy=input('Save this and Select another? [Y/N]','s');
if (dummy=='y'||dummy=='Y')
frequency=[frequency,f_CF];
modeshape=[modeshape,mode];
damping=[damping,damp];
elseif (dummy=='n'||dummy=='N')
break;
end
end
end
%
% Sort and save Extracted Frequencies and Mode-Shapes
%
[frequency,idx_sort]=sort(frequency);
modeshape=modeshape(:,idx_sort);
damping=damping(idx_sort);