Prepare Fourier Amplitude Spectrum from ground motion record (Peak-Acceleration vs time-period)

5 次查看(过去 30 天)
Prepare Fourier Amplitude Spectrum from ground motion record (Peak-Acceleration vs time-period) with the input .txt files which contains peak-acceleration in 1st column & time-period in 2nd column.

回答(1 个)

Star Strider
Star Strider 2024-4-24
Perhaps this —
T1 = readtable('india.19911019...at_output.txt');
T1.Properties.VariableNames = {'Peak Acceleration','Time'}
T1 = 1808x2 table
Peak Acceleration Time _________________ ____ -0.0146 0.02 0.0116 0.04 -0.0029 0.06 0.0296 0.08 0.0185 0.1 0.0295 0.12 -0.0229 0.14 -0.0848 0.16 -0.0696 0.18 -0.048 0.2 -0.0352 0.22 -0.0146 0.24 -0.00165 0.26 0.0493 0.28 0.0701 0.3 0.0688 0.32
[FTs1,Fv] = FFT1(T1.('Peak Acceleration'),T1.Time);
[PeakVal,idx] = max(abs(FTs1)*2);
fprintf('\nMaximum acceleration of about %.4f mm/s^2 occurs at about %0.4f Hz\n', PeakVal, Fv(idx))
Maximum acceleration of about 0.0585 mm/s^2 occurs at about 1.9775 Hz
figure
plot(Fv, abs(FTs1)*2)
grid
xlabel('Frequency')
ylabel('Magnitude')
function [FTs1,Fv] = FFT1(s,t)
t = t(:);
L = numel(t);
if size(s,2) == L
s = s.';
end
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
Fv = Fs*(0:(NFFT/2))/NFFT;
% Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv,:);
end
Experiment!
.
  4 个评论
YQ
YQ 2024-4-27
编辑:YQ 2024-4-27
Thank you for your response. Although I've a deveoped a similar code which gives almost same result as of yours but the problem is when I'm comparing these results from the one obtained from 'Sesimo-Signal' software, the response is completely different. For comparison, I'm uploading the ordinates for the reference purpose.
Star Strider
Star Strider 2024-4-27
My pleasure!
I would not say that the response is completely different, however your non-MATLAB software gives a different magnitude result than my MATLAB code.
Taking ths of that result gives —
SeismoSignalVal = log10(1.4939)
SeismoSignalVal = 0.1743
MyCode = 10^(0.0585)
MyCode = 1.1442
These are not strictly comparable, however they are reasonably close. I am not certain what the ‘SeismoSignal’ code does, or how it produces its Fourier transform results. It also appears to use a 25 Hz lowpass filter prior to calculating the Fourier transform, since there is no energy higher than that in the ‘SeismoSignal’ record. I did not use any filters in my code.
Also, the ‘SeismoSignal’ Bode magnitude plot uses a truncated logarithmic frequency axis. (I have reproduced that here.) My plot uses a simple linear frequency axis.
For all intents and purposes, the ‘SeismoSignal’ result and my result are the same.
My analysis —
figure
imshow(imread('Screenshot 202...27 092404.png'))
title('Image of the ‘SeismoSignal’ Fourier Magnitude Plot')
T2 = readtable('FAS_Ordinates.xlsx', 'VariableNamingRule','preserve')
T2 = 2048x5 table
Frequency Period Fourier Amplitude Fourier Phase Power Spectrum _________ ______ _________________ _____________ ______________ 0.02441 40.96 0.02911 -0.22077 6e-05 0.04883 20.48 0.03319 -1.2565 8e-05 0.07324 13.653 0.03426 0.7691 8e-05 0.09766 10.24 0.01155 -0.37581 1e-05 0.12207 8.192 0.01922 -0.93403 3e-05 0.14648 6.8267 0.02871 -1.4569 6e-05 0.1709 5.8514 0.08533 0.18696 0.00051 0.19531 5.12 0.10611 -1.1405 0.00079 0.21973 4.5511 0.14562 1.3089 0.00149 0.24414 4.096 0.18425 0.44059 0.00238 0.26855 3.7236 0.1323 -0.18164 0.00123 0.29297 3.4133 0.15615 -0.7626 0.00171 0.31738 3.1508 0.28456 1.5609 0.00568 0.3418 2.9257 0.09585 -1.527 0.00064 0.36621 2.7307 0.18708 0.26333 0.00246 0.39063 2.56 0.11011 0.65927 0.00085
VN2 = T2.Properties.VariableNames;
figure
plot(T2.Frequency, T2.('Fourier Amplitude'))
grid
xlabel(VN2{1})
ylabel(VN2{3})
title('Linear Frequency Axis')
[PeakVal,idx] = max(T2.('Fourier Amplitude'));
fprintf('\nMaximum acceleration of about %.4f mm/s^2 occurs at about %0.4f Hz\n', PeakVal, T2.Frequency(idx))
Maximum acceleration of about 1.4939 mm/s^2 occurs at about 1.9775 Hz
figure
semilogx(T2.Frequency, T2.('Fourier Amplitude'))
grid
xlabel(VN2{1})
ylabel(VN2{3})
title('Logarithmic Frequency Axis')
figure
semilogx(T2.Frequency, T2.('Fourier Amplitude'))
grid
xlabel(VN2{1})
ylabel(VN2{3})
xlim([1 max(T2.Frequency)])
title('Truncated Logarithmic Frequency Axis')
.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Frequency Transformations 的更多信息

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by