Site Response Analysis - FFT

Hi,

When I do site response analysis, I make a fourier transform of an earthquake record and then multiply it with an magnification function and get an inverse fourier transform, I get a symmetrical and more cramped time record. What is the problem, how can I fix it. The code is as follows. Thanks.

% Earthquake Data Process
clear all;
dt = 0.005; % Zaman adımı
g = 9.81; % Yerçekimi ivmesi (m/s^2)
data5 = dlmread('Gilroy_No1.txt', '', 4, 0);
data5 = data5';
data = data5(:);
acc = data; % 'g' cinsinden normalize etme işlemi iptal edildi
t = 0.0:dt:dt*(length(data)-1);
% Hız ve yer değiştirme hesaplama
vel = cumsum(acc) * dt * g;
dis = cumsum(vel) * dt;
% Fourier Transform
Ts = dt; 
Fs = 1 / Ts;
ACC = fft(acc, 2^nextpow2(numel(acc)));
N = numel(ACC); 
f = (0:(N-1)) / N * Fs;
m_ACC = abs(ACC) / N;
m_ACC = m_ACC(1:N/2+1);
m_ACC(2:end-1) = 2 * m_ACC(2:end-1); % Tek taraflı spektrum
f = f(1:N/2+1); % Sadece pozitif frekanslar için eksen
% 25 Hz altındaki frekans bileşenlerini almak
indices = f <= 25;
f_limited = f(indices);
m_ACC_limited = m_ACC(indices);
% Zemin Parametreleri (Sert kaya üzerinde uniform, sönümsüz zemin)
h = 3.05; % Kalınlık (m)
Vs = 320; % Vs30 (m/s)
k = h / Vs;
% Amplifikasyon Fonksiyonu
y = 1 ./ abs(cos(k * 2 * pi * f_limited)); % Amplifikasyon Fonksiyonu
% Amplifikasyon Fonksiyonu ile çarpma
y = y';
N_FAS = y .* m_ACC_limited;
% Eksik frekans bileşenlerini sıfır ile doldurma
full_ACC = zeros(size(ACC));
full_ACC(1:length(N_FAS)) = N_FAS;
% Ters Fourier dönüşümü (tek taraflı spektrumdan dönüşüm)
reconstructed_acc = ifft(full_ACC, 'symmetric'); % 'symmetric' seçeneği simetri sağlar
reconstructed_acc= reconstructed_acc*length(reconstructed_acc);
% Hız ve yer değiştirme hesaplama
reconstructed_vel = cumsum(reconstructed_acc) * dt * g;
reconstructed_dis = cumsum(reconstructed_vel) * dt;
% Zaman eksenini yeniden oluşturma
t_reconstructed = (0:length(reconstructed_acc)-1) * dt;
% Sonuçları tek bir plot altında toplama (3x3 layout, transposed)
figure;
subplot(3,3,1), plot(t, acc, '-b');
xlabel('Time (sec)'), ylabel('Original Ground Acceleration (g)');
title('Original Ground Acceleration');
grid on;
subplot(3,3,4), plot(t, vel, '-b');
xlabel('Time (sec)'), ylabel('Velocity (m/sec)');
title('Velocity');
grid on;
subplot(3,3,7), plot(t, dis, '-b');
xlabel('Time (sec)'), ylabel('Displacement (m)');
title('Displacement');
grid on;
subplot(3,3,2), plot(f, m_ACC, '-b');
xlabel('Frequency (Hz)'), ylabel('Original Fourier Amplitude (g-s)');
title('Original Fourier Amplitude (g-s)');
xlim([0 25]);
grid on;
subplot(3,3,5), plot(f_limited, y, '-b');
xlabel('Frequency (Hz)'), ylabel('Amplification Function');
title('Amplification Function');
grid on;
subplot(3,3,8), plot(f_limited, N_FAS, '-b');
xlabel('Frequency (Hz)'), ylabel('Modified Fourier Amplitude (g-s)');
title('Modified Fourier Amplitude (g-s)');
grid on;
subplot(3,3,3), plot(t_reconstructed, reconstructed_acc, '-b');
xlabel('Time (sec)'), ylabel('Reconstructed Ground Acceleration (g)');
title('Reconstructed Ground Acceleration');
grid on;
subplot(3,3,6), plot(t_reconstructed, reconstructed_vel, '-b');
xlabel('Time (sec)'), ylabel('Reconstructed Velocity (m/sec)');
title('Reconstructed Velocity');
grid on;
subplot(3,3,9), plot(t_reconstructed, reconstructed_dis, '-b');
xlabel('Time (sec)'), ylabel('Reconstructed Displacement (m)');
title('Reconstructed Displacement');
grid on;
% Plotları düzenleme
sgtitle('Earthquake Data Analysis');

Answer 1

[Star Strider]

Just looking at your code (without having the data), one problem could be that your amplification function is defined only over the positive half of the Fourier trannsfomr result. I do not know if it is supposed to be symmetric (seismology is not an area of my expertise), howeever making it symmetric by using the appropriate version of flip and then concatenating that to the left side of the orignal vector could be the solution.

To illustrate —

f_limited = 0:0.5:25;
h = 3.05; % Kalınlık (m)
Vs = 320; % Vs30 (m/s)
k = h / Vs;
y = 1 ./ abs(cos(k * 2 * pi * f_limited)); % Amplifikasyon Fonksiyonu

figure
plot(f_limited, y)

f_limited2 = [-flip(f_limited) f_limited];
y2 = [flip(y) y];
grid

figure
plot(f_limited2, y2)
grid

.

[ibrahim can beziklioglu]

Actually, what I want to get is a graph like the one below, because I only amplify at high frequencies in the fourier, all the other data stays the same, but my data ends up more cramped and symmetrical in a way that I don't understand why, the inverse fourier. If you have noticed, both the original fourier and the scaled fourier have a similar character.

I also added earthquake data

[Star Strider]

My initial inclinatiion was to create a FIR filter to filter the signal in the time domain, however while the firls function gave a good accounting of it in the passband, I was not certain how to define the filter above 25 Hz, and I am still not certain of that, so I defaulted to unity gain. In the end, I went with filtering the signal in the frequency domaiin, although that is not ideal.

See if this does what you want —

A1 = readmatrix('Gilroy_No1.txt', 'HeaderLines',4)

A1 = 1600×5

1.0e+00 * 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0004 0.0005 0.0004 0.0004 0.0003 0.0002 0.0002 0.0003 0.0005

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

data5 = A1(:,5);
Fs = 1/0.005; % Zaman adımı
t = linspace(0, size(A1,1)-1, size(A1,1))/Fs;

nnz(ismissing(data5))

ans = 1

data5 = fillmissing(data5,'linear');

figure
plot(t, data5)
grid
xlabel('t')
ylabel('Amplitude')
title('Original Signal')

L = numel(data5);
nfft = 2^nextpow2(L)

nfft = 2048

FTs = fft(data5 .* hann(L), nfft);                                  % Using A Power-Of-2 value For ‘nfft’ Makees ‘fft’ Moree Efficient
FTss = fftshift(FTs);
Fv = linspace(-Fs/2, Fs/2, numel(FTs));

figure
plot(Fv, abs(FTss))
grid

Lidx = (Fv >= -25) & (Fv <= 25);
NF = nnz(Lidx)

NF = 512

% f_limited = 0:0.5:25;
f_limited = linspace(0, 25, fix(NF/2));
Fs

Fs = 200

h = 3.05; % Kalınlık (m)
Vs = 320; % Vs30 (m/s)
k = h / Vs;
% y = zeros(1, Fs/2/max(f_limited)*numel(f_limited))
y = ones(size(Fv));

yf = 1 ./ abs(cos(k * 2 * pi * f_limited)); % Amplifikasyon Fonksiyonu
yf = [flip(yf) yf];                                                     % Create Symmetrical Filter
y(Lidx) = yf;

figure
plot(Fv, y)
grid
xlabel('Frequency')
ylabel('Filter Magnitude')
title('Filter Function (Frequency Domain)')

FTss_filtered = FTss .* y(:);                                       % Filter Signal In Frequency Domain

figure
plot(Fv, abs(FTss_filtered))
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Amplitude Of The Filtered Signal')

FvLv = Fv >= 0;

figure
plot(Fv(FvLv), abs(FTss_filtered(FvLv)))
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Amplitude Of The Filtered Signal (One-Sided)')

data5_filtered = ifft(ifftshift(FTss_filtered));                    % Inveerse Fouriee Transform Of Filtered Signal
data5_filtered = data5_filtered(1:numel(t));                        % Remove Zerro-Padding

[minimag,maximag] = bounds(imag(data5_filtered))

minimag = -0.0023

maximag = 0.0023

figure
plot(t, data5_filtered(1:numel(t)))

Warning: Imaginary parts of complex X and/or Y arguments ignored.

grid
xlabel('t')
ylabel('Amplitude')
title('Filtered Signal')

figure
plot(t, data5, DisplayName='Original')
hold on
plot(t, data5_filtered(1:numel(t)), DisplayName='Filtered')

Warning: Imaginary parts of complex X and/or Y arguments ignored.

hold off
grid
xlabel('t')
ylabel('Amplitude')
title('Original & Filtered Signals')
legend('Location','best')

Ideally, the imaginary part of the inverted siignal should be zero. The reason that it is not is that the filter may not be perfectly symmetrical.

.