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

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

Answer 1

[William Rose]

@Changhyun Kim.,

[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.

[Changhyun Kim]

I'm really sorry and thank you. By the way, I have the correct seismic wave data. But, that is so different from the picture from the book & thesis. What the Rose teacher did was also different from the picture. especially amplitude. Could you please look at the jpg file I attached and help me a little more? plz

[William Rose]

@Changhyun Kim,

Now I see your comment. Let me think about it before I reply.

[Changhyun Kim]

@William Rose

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

Answer 2

[William Rose]

@Changhyun Kim,

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.

[Changhyun Kim]

@William Rose

I want to get the same with jpeg. And there are still problems with different amplitudes. And thank you so much for your concern. You don't know how grateful I am. Where did the problem with amplitude come from? Let's think together. And thank you so much for finding the error in the coding I wrote. It really is

[Changhyun Kim]

@William Rose

But, I don't need 20sample(raw,0~10Hz) in jpg. Because the seismic wave data is fs=100Hz.

so, I hope to 2 traces[(100sample(filtered, 0~20Hz), 100sample(raw, 0~50Hz)] on one plot.

The most important thing is to make the amplitude the same as the jpg, but it keeps failing now

[William Rose]

[edit: forgot to attach the script.]

The script I posted above generates two plots, as shown here:

The third plot in your JPEG was the amplitude spectrum of a signal sampled at 20 Hz. When sampling the signal at 20 Hz, it is essential to first remove frequencies above the Nyquist frequency (), to avoid aliasing. I tried using Matlab's resample(), which has a built-in defualt algorithm to minimize aliasing, but I was not happy iwth the results - too much aliasing. Therefore I filtered the signal myself before resampling. I choose to use filtfilt() with a 6 pole Butterworth. I choose filtfilt() because it has zero phase at all frequencies, and (therefore) does not introduce a time delay. THat is a subject for a separate discussion or reading. I would prefer to oversample, i.e. to do lowpass filtering, then sample at 4 to 5 times the cutoff frequency. This gives me confidence that there will not be aliasing in my sampled signal. But the JPEG you posted for the 20-Hz sample signal shows significnat power in the spectrum up to about 9 Hz, which is very close to the Nyquist frequency of 10 Hz. Therefore, in the script attached, I choose a cutoff frequency of 8 Hz for the anti-aliasing filter. This is closer to the Nyquist frequency than I would normally choose. The script generates the figures below. NOte that the second image is quite similar to the JPEG in your posting.

[William Rose]

@Changhyun Kim,

My answers are laggin behind your comments. I worked on a script to adress your initial request and by the time I finished it and posted it, I saw your next comment. SO I will adress that now.

The reason the amplituds are different in the plots I have generated and the plots in the JPEG is that there ar different ways of normalizinf the FFT. On way is to divide the FFT by N. When I do this, I get the plot below:

Note that the vertical axis has "x10^-3" on it.

[William Rose]

The scle above is closer to the vertical scale in your JPEG but still not the same.

Another option is to normalize by the frequency bin width. THis is done especially with power spectra, so that if you integrate across the spectrum, the units for the integral end up being ower, and not power times Hertz.

[William Rose]

Esitmate the power spectral denisty with cpsd(). Take the square root of thee PSD to get the amplitude spectral density. Units will be , where "pwer" has units of g in your case. Let me think about it.

[Changhyun Kim]

@William Rose

After fft, the frequency axis is enlarged by (sampling time * number of data), and the amplitude axis is enlarged by (half of the number of data). Therefore, scaling is required. However, The amplitude result value of jpg is different from the coding amplitude result that I wrote and you wrote.

[Changhyun Kim]

@William Rose

You divided the amplitude scaling by the number of data, and when i do that, when i input the other seismic wave data, the amplitude result of the coding is completely different from the amplitude of the other seismic wave data Fouries spectrum in the thesis.

Also, As I mentioned it, I don't need 20sample(raw,0~10Hz) in jpg. Because the seismic wave data which i uploaded is fs=100Hz. so, I hope to 2 traces[(100sample(filtered, 0~20Hz), 100sample(raw, 0~50Hz)] on one plot.

[Changhyun Kim]

After reading all the posts sir.Rose sent me, thinking about it myself, and commenting again, I think it caused a misunderstanding in communication.

Although I did not make it in a sophisticated form like the sir Rose's coding, the result of the coding I wrote had the same shape as the jpeg. so, the Fourier spectrum I wanted is the Fourier spectrum that is not much different from the amplitude values of the Fourier spectrum in the thesis.

[William Rose]

@Changhyun Kim,

You can compute an estimate of the power spectral density with

[Pxx,fpsd]=cpsd(x,x,[],[],[],fs);

where x and fs are defined in the code posted previously.

Then you can plot the ampltude spectrum (where amplitude=sqrt(power)) as follows:

figure;
plot(fpsd,sqrt(Pxx));
xlabel('Frequency (Hz)'); ylabel('Amplitude [g/sqrt(Hz)]'); grid on

This gives the plot below. The vertical axis values are sigificantly different than the values in the JPEG you posted.

If you normalize the amplitude of the FFT (see previous comment) by 2.5*N, istead of normalizing by N, you get a spectrum whose values are very similar to the JPEG. See comparison below of the JPEG and the plot normalized by 25*N. Are you confident that the time domain signal you are using is the same as the signal used to generate the JPEG? I think the signal used to generate the JPEG is not the same as the signal in GYEONGJU(USN).xlsx. I say that becausethe relative heights of peaks is not quite the same for the two. See the circled regions in the plots below. Note that the middle peak in the circled regions have the same amplitude(0.0025), but the peaks on either side are somewhat dfferent. Perhaps the amplitude or duration of the JPG signal differs from GYEONGJU(USN). (The duration matters: the inclusion of a longer section of the "noise" before or after the event will affect the spectrum amplitude.)

[Changhyun Kim]

@William Rose

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]

@Changhyun Kim,

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.