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

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.

Answer 1

[Star Strider]

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!

.

[YQ]

[FTs1,Fv] = FFT1(T1.('Peak Acceleration'),T1.Time);

Here, why you've used FFT for time function ?

[PeakVal,idx] = max(abs(FTs1)*2);

Why the absoulte value of FTs1 is multiplied by 2?

[Star Strider]

Because the instructions are:

‘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.’

and to me this implies using the fft function.

The fft calculation splits the energy of the signal into half in the ‘negative’ frequency section and half in the ‘positive’ frequency section. Calculating a one-sided Fourier transform then requires that the magnitude be multipliied by 2 in order to approximate the original signal amplitude.

.

[YQ]

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]

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')

.