using acceleration record, Get velocity and displacement signal & carry out FFT and plot the spectra

Given information

• Time domain signal file : GEE Earthquke.xlsx (attached)
• Sampling rate (frequency): 64Hz
• Acceleration values in g (9.81m/s^2)

1) Velocity and displacement signal

• Consider baseline correction
• Get the relative velocity and displacement signals
• Plot the signals
• Get the amplitude parameters in time domain (peak acceleration, peak velocity and peak displacement)

2) Carry out the FFT and plot the spectra

P.S No matter how much I think and try, I can't code.. plz help me

Answer 1

[Mathieu NOE]

hello

see my suggestion below

clc
clearvars
data = readmatrix('GEE Earthquake.xlsx');
t = data(:,1);
dt = mean(diff(t));       % Average dt
fs = 1/dt;                  % Frequency [Hz] or sampling rate

% Acceleration data
acc = data(:,2)*9.81;         % Add gravity factor ( g to m/s²)
acc = detrend(acc);     % baseline correction

% some additionnal high pass filtering     % baseline correction
N = 2;
fc = 0.25; % Hz
[B,A] = butter(N,2*fc/fs,'high');
acc = filter(B,A,acc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
velocity = cumtrapz(dt,acc);
velocity = detrend(velocity);    % baseline correction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp     = cumtrapz(dt,velocity);
disp = detrend(disp);    % baseline correction

%
figure(1)
subplot(311),plot(t,acc);
title('accel'); ylabel('m/s²');xlabel('time (s)');
subplot(312),plot(t,velocity);
title('velocity'); ylabel('m/s');xlabel('time (s)');
subplot(313),plot(t,disp);
title('disp'); ylabel('m');xlabel('time (s)');


%%fft
nfft = length(acc);
fft_spectrum = (abs(fft(disp,nfft))/nfft);

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*fs/nfft;

figure(2)
plot(freq_vector,fft_spectrum);
title('disp'); ylabel('m');xlabel('time (s)');

xlim([0 10]);

[Changhyun Kim]

thnx for ur help,

Meanwhile, given sample frequency value fs is not 1/dt but 65Hz

After fixing fs to 65, I don't know what criteria to set fc.

fc means arbitrary frequency?

would you plz re-write or revise comment considering this

[Mathieu NOE]

hello again

1/ yes I know fs = 64 Hz (and so it is in my code), simply I computed it from your time data (in your file)

so I don't have to write it explicitely as I prefer to compute it from the data file;

2/ to the second point (fc) :

if I only use detrend, I was not 100% satisfied with the result (disp) , so I added some high pass filter on the acceleration. The cut off frequency fc is choosen by trial and error so that the curve of displacement is almost flat / zero where the acceleration is almost zero. If you reduce fc you can see some more oscillations at very low frequency , if fc is too high thene you will reduce the peak amplitude of the signal. BTW the fft spectrum is a good indication of the frequency content of the disp signal and we can see the energy below 0.2 Hz drops so we can put the cut off frequency around this value

minor bug : the fft window x axis is of course a frequency and not time

figure(2)
plot(freq_vector,fft_spectrum);
title('disp'); ylabel('m');xlabel('Frequency (Hz)');

[Changhyun Kim]

I really wanna thank u for ur help:)

[Mathieu NOE]

my pleasure !

[Changhyun Kim]

in fact,, it doesn't work. The built-in function called butter keeps getting an error and the script cannot be executed.

[Mathieu NOE]

hello @Changhyun Kim

what is the error message ?