Site Response Analysis - FFT

38 次查看（过去 30 天）
ibrahim can beziklioglu
ibrahim can beziklioglu 2024-10-29，22:06
评论： Star Strider 2024-11-5，3:11
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');

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2024-10-29，23:46
Ran in:
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
.
  2 个评论
ibrahim can beziklioglu
ibrahim can beziklioglu 2024-10-30，6:59
移动：Star Strider 2024-11-3，22:32
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
Star Strider 2024-11-5，3:11
Ran in:
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.
.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Geology 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by