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

25 visualizzazioni (ultimi 30 giorni)
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

Risposta accettata

William Rose
William Rose il 12 Ago 2022
Modificato: William Rose il 12 Ago 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 Commenti
Changhyun Kim
Changhyun Kim il 12 Ago 2022
Modificato: Changhyun Kim il 12 Ago 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);

Accedi per commentare.

Più risposte (1)

William Rose
William Rose il 12 Ago 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 Commenti
Changhyun Kim
Changhyun Kim il 12 Ago 2022
Modificato: Changhyun Kim il 12 Ago 2022
Thank you very much. As you said, when I converted the seismic wave raw data into python, there seems to be a bit of noise when sampling because of the dt problem or the duration problem. Thank you so much for helping me. I was really struggling with this issue. I am not an intellectual like Sir Rose in fft, so there are a lot of things that I am lacking. I'm new to the power spectral part. And since the MATLAB script file you sent me is very different now, it seems difficult for me to follow by myself. I'm really sorry, but could you please attach the matlabscript file? I hope to closer look and understand.
William Rose
William Rose il 12 Ago 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.

Accedi per commentare.

Community Treasure Hunt

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

Start Hunting!

Translated by