how to convert the seismic wave data in Excel to Fourier spectrum using filtering and sampling

23 views (last 30 days)
Given information
  • Time domain signal file : GYEONGJU(USN).xlsx (attached)
  • Sampling rate (frequency): 100Hz
  • Acceleration values in g (9.81m/s^2)
I want to convert the ground acceleration data in the attached Excel file to the frequency [Hz] on the horizontal axis and the magnitude on the vertical axis by performing FFT in MATLAB.
Also, it must plot only up to the Nyquist frequency. i.e 100sample(0~50Hz)
In addition, I need the Fourier Spectrum as raw data and the Fourier Spectrum filtered by performing FFT on 100sample seismic wave data with a 20Hz low-pass filter. In other words, I need the 100sample(filtered, 0~20Hz) Fourier Spectrum and the 100sample(raw,0~50Hz) Fourier Spectrum both.
The dominant frequency of this earthquake should be 6 to 10 Hz.
I want to get the Fourier Spectrum for 100sample(filtered, 0~20Hz) and 100sample(raw, 0~50Hz) from the attached jpg file. I tried but failed
I'm sorry, but I'd appreciate your help. plz, help me

Accepted Answer

William Rose
William Rose on 12 Aug 2022
Edited: William Rose on 12 Aug 2022
[edit: I adjusted the frequency-domain filtering to make the filtered X symmetric about the Nyquist frequency. This does not affect the plot shown, because the plot shown only goes up to the Nyquist frequency. However, it has the benefit that if you do the invese FFT of Xfilt, you will now get a real result, as you would desire.]
data=xlsread('GYEONGJU(MKL).xlsx');
N=length(data);
t=data(:,1); %vector of times (s)
x=data(:,2); %vector of x-values
dt=(t(end)-t(1))/(N-1); %sampling interval (s)
fs=1/dt;
df=fs/N;
f=(0:N-1)*df; %vector of frequencies (Hz)
X=fft(x); %compute FFT(x)
%% Plot results
plot(f,abs(X),'-b')
xlabel('Frequency(Hz)'); ylabel('|X(f)|');
grid on
xlim([0,fs/2]);
That is the amplitude spectrum of the unfiiltered seismogram. YOu can tell by inspecting th plot that if we filter frequencies above 20 Hz, it will not be the case that the dominant frequency is 6-10 Hz. The energy seems to be spread rather broadly from 1 to 16 Hz (ssuming the portion above 20 Hz is filtered). You can filter in the frequency domain or th time domain. Filtering in the frequency domain is easier, once you have the FFT: simply set the FFT to zero for frequencies above 20 Hz.
Xfilt=X;
Xfilt(f>20 & f<(fs-20))=0;
hold on; plot(f,abs(Xfilt),'-r');
legend('Unfilt.','Filtered');
Try that.
  3 Comments
Changhyun Kim
Changhyun Kim on 12 Aug 2022
Edited: Changhyun Kim on 12 Aug 2022
Even the coding I wrote down now has the wrong amplitude, but the graph shape is the same.
.
.
clc; clear;
data = xlsread('GYEONGJU(USN).xlsx');
t = data(:,1);
dt = mean(diff(t));
fs = 1/dt;
y = data(:,2);
plot(t,y);
Y = fft(y);
F = abs(Y);
subplot(3,1,1)
plot(F);
n = length(F)
subplot(3,1,2)
F1 = F(1:n/2)
plot(F1);
subplot(3,1,3)
F2 = F1/(n/2)
f = (1:n/2)/(dt*n)
plot(f,F2);

Sign in to comment.

More Answers (1)

William Rose
William Rose on 12 Aug 2022
In the script you posted, the first plot you generated was over-written by the second plot you generated. Adding a "figure" command before the second plot will fix this.
You generated a plot with three panels, but the JPEG has 3 traces on one plot. WHich do you want?
I am attaching a script that generates a time domain plot. This is the plot that was overwritten in your script, and I have added labels to it. It also generates a plot with two amplitude spectra on a single plot, as in the JPEG which yu supplied. The two spectra are the unfiltered 100 Hz signal and the 100 Hz signal, filtered at 20 Hz.
I will discuss resampling, and I will add the plot of the FFT of the resamled signal, in a comment.
  12 Comments
William Rose
William Rose on 12 Aug 2022
You are welcome. Here is the script that generates the figures.
I hve added a number of explanatory comments.
The final section of the script plots results. You will notice that the plots in Figure 2 include normalization of the Y-values by 2.5*N. The factor of 2.5 is chosen to give FFT amplitudes that match the JPEG image.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!