FRF Issue of 4 DOF Frame

Hi everybody, I am triying to analyze a 4 story shear building under and earthquake. I used state space formulation and obtained the displacement of each story. When I try to find FRF to observe either I can find the modal parameters from response histories or not, I can get no meaningful results. I do fft(Output)/fft(Input) but ı obtained meaningless results. Where ı do a mistake in my code ?

State space portion :

global CC EQ M

K = [10000 -5000 0 0;-5000 10000 -5000 0;0 -5000 10000 -5000;0 0 -5000 5000]; %in N/m2

M = [4000 0 0 0;0 4000 0 0;0 0 4000 0;0 0 0 4000]; % in kg

[V,WW] = eig(K,M);

W = sqrt(WW)*[1;1;1;1];

C = inv(transpose(V))*0.1*[W(1) 0 0 0;0 W(2) 0 0;0 0 W(3) 0;0 0 0 W(4)]*inv(V);

CC = [zeros(4) eye(4);-inv(M)*K -inv(M)*C];

EQ = xlsread('el_centro');

step = 0.02;

time = 0:step:max(EQ(:,1));

y0 = [0 0 0 0 0 0 0 0]; % [d1 d2 d3 d4 v1 v2 v3 v4]

[tsol,ysol] = ode23('testode_2D',time,y0);

plot(time,ysol(:,1:4))

title('El centro Floor Response')

legend('Fisrt','Second','Third','Forth')

xlabel('Time')

ylabel('Displacement (m)')

.............................................................................

function dy = testode_2D(t,y)

global CC EQ M

index = interp1(EQ(:,1),1:length(EQ(:,1)),t,'nearest');

if (index ~= 1) & (index ~= length(EQ(:,1)))

if t > EQ(index,1)

oran = (t - EQ(index,1))/0.02;

ag = (EQ(index+1,2) - EQ(index,2))*oran + EQ(index,2);

else

oran = (t - EQ(index-1,1))/0.02;

ag = (EQ(index,2) - EQ(index-1,2))*oran + EQ(index-1,2);

end

else

ag = EQ(index,2);

end

f1 = -M(1,1)*ag;

f2 = -M(2,2)*ag;

f3 = -M(3,3)*ag;

f4 = -M(4,4)*ag;

A00 = zeros(4);

FF = [A00;inv(M)]*[f1;f2;f3;f4];

dy = CC*y + FF;

FRF part :

EQ = xlsread('el_centro');

EQ(:,2) = 9.81*EQ(:,2);

EQ_floor_dat = xlsread('el_centro_dat');

fs = 50; % sampling freq

T = .02; % sampling period

N = length(EQ_floor_dat(:,1)); % number of signal data

time = EQ_floor_dat(:,1); % time vector

Output = EQ_floor_dat(:,5) ;

Input = EQ(:,2);

Y = fft(Output); % output fft

X = fft(Input); %input fft

P2 = abs(Y/N);

P1 = P2(1:N/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = fs*(0:(N/2))/N;

[Mathieu NOE]

hello

we need the input data as well (el_centro file)

[Can AKYOL]

Hello, here it is sir !

Answer 1

[Mathieu NOE]

hello again

so this is the code a bit reworked and expanded

the "interesting" portion in the FRF plot is in the sub Hertz range.. so the quality of the plot is related to simulation time length. It would require to simulate on longer records to further improve the FRF plot below 1 Hz

Code :

clc
clearvars

global CC EQ M
K = [10000 -5000 0 0;-5000 10000 -5000 0;0 -5000 10000 -5000;0 0 -5000 5000]; %in N/m2
M = [4000 0 0 0;0 4000 0 0;0 0 4000 0;0 0 0 4000]; % in kg
[V,WW] = eig(K,M);
W = sqrt(WW)*[1;1;1;1];
C = inv(transpose(V))*0.1*[W(1) 0 0 0;0 W(2) 0 0;0 0 W(3) 0;0 0 0 W(4)]*inv(V); 
CC = [zeros(4) eye(4);-inv(M)*K -inv(M)*C];
EQ = xlsread('el_centro.xlsx');
% step = 0.02;
% time = 0:step:max(EQ(:,1));
time = EQ(:,1);
Input = EQ(:,2);

step = mean(diff(time));

y0 = [0 0 0 0 0 0 0 0];  % [d1 d2 d3 d4 v1 v2 v3 v4]
% [tsol,ysol] = ode23('testode_2D',time,y0);
[tsol,ysol] = ode23(@testode_2D,time,y0);

%% time plot
figure(1),
subplot(211),plot(time,Input)
title('El centro Base Input')
xlabel('Time')
ylabel('Displacement (m)')
subplot(212),plot(time,ysol(:,1:4))
title('El centro Floor Response')
legend('Fisrt','Second','Third','Forth')
xlabel('Time')
ylabel('Displacement (m)')

%% FRFs computations
% NFFT = 512;
% NOVERLAP = 0.5*NFFT;
% NFFT = length(time);
NFFT = length(time); % maximal frequency resolution (but implies zero overlap)
NOVERLAP = 0;
WINDOW = hanning(NFFT);
Fs = 1/step;

for ci =1:4
    [Txy(:,ci),Freq] = tfestimate(Input,ysol(:,ci),WINDOW,NOVERLAP,NFFT,Fs);
end

figure(2),
loglog(Freq,abs(Txy));
title('El centro Floor Response FRF vs. Base Input')
legend('Fisrt','Second','Third','Forth')
xlabel('Frequency (Hz)')
ylabel('Displacement ratio (m / m)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = testode_2D(t,y)
global CC EQ M
index = interp1(EQ(:,1),1:length(EQ(:,1)),t,'nearest');
if (index ~= 1) & (index ~= length(EQ(:,1)))
    if t > EQ(index,1)
      oran = (t - EQ(index,1))/0.02;
     ag = (EQ(index+1,2) - EQ(index,2))*oran + EQ(index,2);
    else
     oran = (t - EQ(index-1,1))/0.02;
     ag = (EQ(index,2) - EQ(index-1,2))*oran + EQ(index-1,2);
    end
else
    ag = EQ(index,2);
end
f1 = -M(1,1)*ag;
f2 = -M(2,2)*ag;
f3 = -M(3,3)*ag;
f4 = -M(4,4)*ag;
A00 = zeros(4);
FF = [A00;inv(M)]*[f1;f2;f3;f4];
dy = CC*y + FF;
end

[Can AKYOL]

Thank you so so much !

[Can AKYOL]

Actually now I wonder why I cannot get my 4 exact frequencies from that plot. do you have any idea for that ?

[Mathieu NOE]

Hello

I changed the excitation signal to a white noise (random) , so now we can control the sampling frequency and time length; also a white noise has constant power vs frequency, which is not the case of a "true" earthquake signal; so we apply here the same techniques as in modal analysis (excitation signals are either broadband or short pulses or chirp sine)

we can get a much prettier FRF plot now :

code :

clc
clearvars

global CC EQ M
K = [10000 -5000 0 0;-5000 10000 -5000 0;0 -5000 10000 -5000;0 0 -5000 5000]; %in N/m2
M = [4000 0 0 0;0 4000 0 0;0 0 4000 0;0 0 0 4000]; % in kg
[V,WW] = eig(K,M);
W = sqrt(WW)*[1;1;1;1];
C = inv(transpose(V))*0.1*[W(1) 0 0 0;0 W(2) 0 0;0 0 W(3) 0;0 0 0 W(4)]*inv(V); 
CC = [zeros(4) eye(4);-inv(M)*K -inv(M)*C];
% EQ = xlsread('el_centro.xlsx');
% step = 0.02;
% time = 0:step:max(EQ(:,1));

samples = 1e4;
step = 1/2;
time = (0:samples-1)'*step;
Input = randn(samples,1);
EQ = [time Input];


y0 = [0 0 0 0 0 0 0 0];  % [d1 d2 d3 d4 v1 v2 v3 v4]
% [tsol,ysol] = ode23('testode_2D',time,y0);
[tsol,ysol] = ode23(@testode_2D,time,y0);

%% time plot
figure(1),
subplot(211),plot(time,Input)
title('Base Input')
xlabel('Time')
ylabel('Displacement (m)')
subplot(212),plot(time,ysol(:,1:4))
title('Floor Response')
legend('Fisrt','Second','Third','Forth')
xlabel('Time')
ylabel('Displacement (m)')

%% FRFs computations
NFFT = 2048*2;
NOVERLAP = 0.5*NFFT;
% NFFT = length(time); % maximal frequency resolution (but implies zero overlap)
% NOVERLAP = 0;
WINDOW = hanning(NFFT);
Fs = 1/step;

for ci =1:4
    [Txy(:,ci),Freq] = tfestimate(Input,ysol(:,ci),WINDOW,NOVERLAP,NFFT,Fs);
end

figure(2),
loglog(Freq,abs(Txy));
title('Floor Response FRF vs. Base Input')
legend('Fisrt','Second','Third','Forth')
xlabel('Frequency (Hz)')
ylabel('Displacement ratio (m / m)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = testode_2D(t,y)
global CC EQ M
index = interp1(EQ(:,1),1:length(EQ(:,1)),t,'nearest');
if (index ~= 1) & (index ~= length(EQ(:,1)))
    if t > EQ(index,1)
      oran = (t - EQ(index,1))/0.02;
     ag = (EQ(index+1,2) - EQ(index,2))*oran + EQ(index,2);
    else
     oran = (t - EQ(index-1,1))/0.02;
     ag = (EQ(index,2) - EQ(index-1,2))*oran + EQ(index-1,2);
    end
else
    ag = EQ(index,2);
end
f1 = -M(1,1)*ag;
f2 = -M(2,2)*ag;
f3 = -M(3,3)*ag;
f4 = -M(4,4)*ag;
A00 = zeros(4);
FF = [A00;inv(M)]*[f1;f2;f3;f4];
dy = CC*y + FF;
end

[Can AKYOL]

Hello, thank you again,

Now everything is much more clear

[Mathieu NOE]

My pleasure

read some modal analysis publications if you need to refine your skills in that matter :)