clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [rn, LR, UR, frfm] = lsfd(w, frf, fn, xin)
% Least-Squares Frequency-Domain (LSFD) estimator for the residues
%
% k=n
% -----
% \ r_k LR
% H(s) = . -------------------------------------- + ------ + UR
% / s^2 + 2 s w_k xi_k + w_k^2 s^2
% -----
% k=0
%
% w natural frequency vector in rad/s
% frf complex frequency response function vector
% fn eigen frequency in Hz obtained using lsfc.m
% xin modal damping factor obtained using lsfc.m
%
% r_k complex residues
% LR and UR lower and upper residual terms
% frfm identified FRF in modal basis using s-model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('Input_Run1_-X_26June2023');
load('FRF_DisSmooth_9June2023');
load('FRF_Amp_DR_27July2023');
%%%%%%%%%%%%%%%%%%%%%%
DataX=[FRF_DisSmooth_9June2023(:,1) FRF_DisSmooth_9June2023(:,2)];
DataY=[FRF_DisSmooth_9June2023(:,3) FRF_DisSmooth_9June2023(:,4)];
DataZ=[FRF_DisSmooth_9June2023(:,5) FRF_DisSmooth_9June2023(:,6)];
Impulse = Time_Point1__X.y_values.values;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wn_Hz = FRF_Amp_DR_27July2023{1,1}(1,:);
zita = FRF_Amp_DR_27July2023{1,1}(2,:);
FRF_Amplitude = FRF_Amp_DR_27July2023{1,1}(3,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FFT analysis
FRF_tot = zeros(size(DataX,1), 6);
Frequency_tot = zeros(6, size(DataX,1));
N = size(DataX,2);
i=1;
[m,n]=size(DataX(:,i));
Fs = 1024;
Ts = 1/Fs;
L1 = m;
% t = (0:L1-1)*Ts;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = DataX(1:L1,i);
P2 = abs(X/L1);
P1 = P2(1:L1);
% P1(2:end-1) = 2*P1(2:end-1);
P11=2*P1;
f1 = Fs*(0:(L1)-1)/L1;
fxx(i,:) = f1(1,:)*2*pi; % rad/s;
frfx(:,i) = P11(:,1)*Fs;
FRF_tot(:,i*3-2) = frfx(:,i);
Frequency_tot(i*3-2,:) = fxx(i,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[w1, frf1] = select_frf(fxx(1,:)', frfx(:,1), 1, 700);
figure(1)
hold on
plot(w1/(2*pi), 20*log10(abs(frf1)), 'LineWidth', 1, ...
'Color', [0, 1, 1], 'DisplayName', 'Measured FRF')
xlim([1 700])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wn_Hz = wn_Hz(1,[1:end]);
zita = zita(1,[1:end]);
fns_x = wn_Hz;
xin = zita;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = 1j*w1;
wn = 2*pi*fns_x ;
nwn = length(wn) ;
L = length(w1) ;
% linear square modal constant estimation
Cls = zeros(L, nwn+2) ;
for ic=1:nwn+2
if ic <= nwn
Cls(:, ic) = 1./(s.^2+2*s*wn(ic)*xin(ic)+wn(ic)^2) ;
elseif ic == nwn+1
Cls(:, ic) = 1./(s.^2) ;
elseif ic == nwn+2
Cls(:, ic) = 1 ;
end
end
rntot = linsolve(Cls, frf1) ;
% residue
rn = rntot(1:nwn, 1) ;
% lower residual term
LR = rntot(nwn+1, 1) ;
% upper residual term
UR = rntot(nwn+2, 1) ;
% identified modal FRF in s-domain
frfm = zeros(L, nwn) ;
for ic=1:nwn
frfm(:, ic) = rn(ic)./(s.^2+2*s*wn(ic)*xin(ic)+wn(ic)^2) ;
end
frfm = sum(frfm, 2)+LR./(s.^2)+UR ;
plot(w1/(2*pi), 20*log10(abs(frfm)), ...
'r', 'DisplayName', 'Identified FRF using s-model')
legend('location', 'southwest')
