%Rational Fractional Polynomial Method for extracting Modal Properties
method = 'rfp';
frf = x_dft;
dofm = menu('Approach','MDOF','SDOF');
if dofm == 1
diary (sprintf('%s_%1.0d_Mdof_results.m',method,E1))
disp('MDOF Simulated Data')
else
diary (sprintf('%s_%1.0d_Sdof_results.m',method,E1))
disp('SDOF Simulated Data')
end
disp([sprintf('\n')])
disp('RFP METHOD')
disp([sprintf('\n')])
disp('Two Close Modes')
disp([sprintf('\n')])
disp('Natural Frequencies and Damping Ratios for the data with two close modes')
Natural_frequency_Damping_ratio = [ fn1 E1 ; fn2 E2 ; fn3 E3 ; fn4 E4 ; fn5 E5 ];
W_TR = length(frf)/2;
N = length(frf)/2;
frf = (frf(1:W_TR)); % The HP analyzer just gives the positive freq components
f = f(1:W_TR);
df = f(3) - f(2);
% Before getting to this point I need to select the data
% and need freq. matrix
%----------- Specifying the freq range of curve fit -----------%
wind = menu('Do you want to specify the freq. range for the curve fit?','Yes','No');
if wind == 1
disp('Frequency Range Specified')
specify = menu('How do you want to specify the freq. range?','Point on Graph','Type it');
if specify == 1
figure(fig + 1)
semilogy(f(1:W_TR),abs(frf(1:W_TR)))
title('Select the first point (minimum frequency)')
[x_frm1,y]=ginput(1);
figure(fig + 1)
semilogy(f(1:W_TR),abs(frf(1:W_TR)))
title('Select the first point (maximum frequency)')
[W_TR,y]=ginput(1);
sprintf('The selected frequency range is:\n\tMinimum freq = %8.4g\n\tMaximum freq = %8.4g',x_frm1,W_TR)
else
figure(fig + 1)
semilogy(f(1:W_TR),abs(frf(1:W_TR)))
x_frm1 = input('Minimum Frequency (Hz): ');
W_TR = input('Maximum Frequency (Hz): ');
sprintf('The selected frequency range is:\n\tMinimum freq = %d\n\tMaximum freq = %d',x_frm1,W_TR)
end
% Isolating the Frequency Range
x_frm1 = round(x_frm1/df + 1);
W_TR = round(W_TR/df + 1);
frf_F1 = ones(x_frm1-1,1); % Putting ones before the isolated FRF
components (mult by mean of irf)
frf_F1(x_frm1:W_TR) = frf(x_frm1:W_TR); % Isolated FRF components
frf_F1(W_TR+1:N) = ones(N-(W_TR),1); % Putting ones after the isolated FRF
components (mult by mean of irf)
%---------------------%
[r,c] = size(frf_F1);
if r>c
frf = conj(frf_F1');
else
frf = frf_F1;
end
clear frf_F1
else
x_frm1 = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(fig + 1)
semilogy(f(x_frm1:W_TR),abs(frf(x_frm1:W_TR)))
dof = input('How many DOF?: '); %%%DOF
diary off % Turns diary OFF
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------- Processing Data --------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w = 2*pi*f(x_frm1:W_TR);
% Scaling the frequency from 0 to 1
% Dividing by the maximum frequency
wi = w/max(w);
n = dof*2;
wt = ones(1,length(frf(x_frm1:W_TR)));
iter = 100;
tol=0;
[A,B] = invfreqs(frf(x_frm1:W_TR),wi,n,n,wt,iter,tol);
[R_rfp,P_rfp,K] = residue(A,B); % Residues & Poles, respectively
%--- Calculating the Natural Freq & Damping Ratio ---%
Damp_ratio_rfp = -real(P_rfp)./(abs(P_rfp));
% Here the natural frequency is multiplied by
% the maximum in because the frequencies were
% scaled from 0 to 1 to avoid problems in
% the invfreqs function
Fn_rfp = abs(P_rfp)*max(w)/(2*pi);
%---------- Calcualting the FRF Curve Fit ----------%
frf_rfp = freqs(A,B,wi);
%--- Adding the conjugate components to the FRF and zeros in the truncated ---%
% Experimental FRF
frf(1:x_frm1-1) = 0;
frf(N+1) = 0;
frf(N+2:2*N) = conj(frf(N:-1:2));
% Curve Fit
frf_rfp(x_frm1:W_TR) = frf_rfp;
frf_rfp(1:x_frm1 -1) = zeros(1,x_frm1 -1);
frf_rfp(N+1) = 0;
frf_rfp(N+2:2*N) = conj(frf_rfp(N:-1:2));
%--- Calculating the Impulse Response Function from the FRF Inverse ---%
irf = real(ifft(frf));
irf_rfp = real(ifft(frf_rfp));
%--------- Calculating the Residual ---------%
%-- FRF --%
Residual = frf(x_frm1:W_TR) - frf_rfp(x_frm1:W_TR);
%-- IRF --%
ResidualT = irf - irf_rfp;
%--------- Plotting ---------%
figure(fig + 2)
semilogy(f(x_frm1:W_TR),abs(frf(x_frm1:W_TR)),f(x_frm1:W_TR),abs(frf_rfp(x_frm1:W_TR)))
figure(fig + 3)
plot(f(x_frm1:W_TR),angle(frf(x_frm1:W_TR)),f(x_frm1:W_TR),angle(frf_rfp(x_frm1:W_TR)))
figure(fig + 4)
plot(f(x_frm1:W_TR),imag(frf(x_frm1:W_TR)),f(x_frm1:W_TR),imag(frf_rfp(x_frm1:W_TR)))
figure(fig + 5)
subplot(2,1,1)
semilogy(f(x_frm1:W_TR),abs(frf(x_frm1:W_TR)))
subplot(2,1,2)
plot(f(x_frm1:W_TR),abs(Residual))
figure(fig + 6)
t = linspace(0,1/df,2*N);
plot(t,irf,t,irf_rfp)
figure(fig + 7)
plot(t,(ResidualT))
diary on % Turns diary ON
%--------- Displaying Results ---------%
Residues_and_Poles_rfp = [ R_rfp P_rfp ];
format long g
Natural_freq_Damping_ratio_rfp = [ Fn_rfp Damp_ratio_rfp ];
%--------- Calculating & Displaying the Standard Deviation ---------%
% Frequecy domain (Residual)
Curvefit_Standard_deviation_FRF = sqrt((Residual*Residual')/length(Residual));
% Time domain (ResidualT)
Curvefit_Standard_deviation_IRF = sqrt((ResidualT*ResidualT')/length(ResidualT));
% Around each natural frequency in the FRF
stdn = menu('Standard Deviation around each nat. freq.','Yes','No');
if stdn == 1
std_nat_freq_rfp
end
diary off % Turns diary off
if dofm == 1
sprintf('The Output file (results) is %s_%1.0d_Mdof_results.m',method,E1)
else
sprintf('The Output file (results) is %s_%1.0d_Sdof_results.m',method,E1)
end
% Hilbert Envelope Method
format long g
%%%%%%%%%%%%%%%%%%%%%%