How can i get frequency domain of an earthquake?

Question

ADNAN KIRAL 2021-1-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/728648-how-can-i-get-frequency-domain-of-an-earthquake

评论： Star Strider 2023-6-3

采纳的回答： Star Strider

EQ.txt

在 MATLAB Online 中打开

Hi,

I need to get the frequency domain of an earthquake. I have used "FFT" in Matlab, but I could not get correctly.

the earthquake is attached here called EQ.txt.

dt=0.01 sec time step;

Can you please correct that?

Thanks for your help.

load EQ.txt 
ti =EQ(:,1);
acc=EQ(:,1);
N=length(ti); 
y=fft(acc,N);
y=(y)/(N/2);
figure; plot(abs(y))

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2021-1-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/728648-how-can-i-get-frequency-domain-of-an-earthquake#answer_607668

编辑：Star Strider 2021-1-27

在 MATLAB Online 中打开

Try this:

acc = readmatrix('EQ.txt');
L = numel(acc);
Ts = 0.01;
Fs = 1/Ts;
Fn = Fs/2;
t = linspace(0, L, L)*Ts;
FTacc = fft(acc)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(t, acc)
grid
xlabel('Time (sec)')
ylabel('Acceleration (Units)')
figure
plot(Fv, abs(FTacc(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Acceleration Amplitude (Units)')
xlim([0  15])
figure
plot(Fv, (abs(FTacc(Iv))*2).^2)
grid
xlabel('Frequency (Hz)')
ylabel('Acceleration Power (Units^2)')
xlim([0  15])

EDIT — (27 Jan 2021 at 18:50)

Corrected typographical errors.

11 个评论
显示 9更早的评论隐藏 9更早的评论

Star Strider 2021-1-27

在 MATLAB Online 中打开

My pleasure.

I have no idea how SeismoSignal calculates its fft results. If it zero-pads the signal or does other processing on it before it calculates the Fourier transform, the amplitudes may differ slightly from those that the MATLAB fft function calculates. My code correspoinds to the way MATLAB suggests that the fft results be scaled and presented, and it works correctly with other signals. (I have used the same essential code that I posted here for several years.)

I trust MATLAB in this respect. It is always possible that SeismoSignal is producing incorrect results.

Testing my code with a simple sine signal:

acc = sin(2*pi*2*(0:0.01:100));

produces the correct result with the frequency at 2 Hz and an amplitude of 0.99929 for the amplitude, and 0.998585 for the power. (You can easily simulate it with my code as well.) It might be interesting to see what SeismoSignal returns for it. (You would have to create this in MATLAB and write it to a text file to export it to SeismoSignal to read.)

Star Strider 2021-3-3

As always, my pleasure!

Star Strider 2023-6-3

在 MATLAB Online 中打开

EQ.txt

My approach to calculating the fft has changed slightly since I wrote this.

I would now calculate it as —

acc = readmatrix('EQ.txt');

L = numel(acc);

Ts = 0.01;

Fs = 1/Ts;

Fn = Fs/2;

t = linspace(0, L, L)*Ts;

NFFT = 2^nextpow2(L);

FTacc = fft(acc.*hann(L), NFFT)/L;

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

[pks,locs] = findpeaks(abs(FTacc(Iv))*2, 'MinPeakProminence',0.04);

figure

plot(Fv, abs(FTacc(Iv))*2)

grid

xlabel('Frequency (Hz)')

ylabel('Magnitude')

xlim([0 15])

text(Fv(locs),pks, sprintf('\\leftarrow Magnitude = %.3f Units\n Frequency = %.3f Hz', pks, Fv(locs)))

The principal differences from my earlier code are the use of zero-padding to the next power of 2 beyond the length of ‘L’ to improve the efficiency of the fft calculation (it incidentally increases the frequency resolution, that in my opinion is always preferable), and using a window (in this instance the hann window) to correct for the fft being finite rather than infinite. Together, they produce a better result than my earlier code. You can of course do other changes as well, such as using the loglog instead of linear scales, if that’s preferable.

I added the findpeaks and text calls out of interest

.

请先登录，再进行评论。

Answer 2

Paul 2023-6-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/728648-how-can-i-get-frequency-domain-of-an-earthquake#answer_1249439

在 MATLAB Online 中打开

EQ.txt

Hi ADNAN,

I had a comment in this thread that was mysteriously deleted, so I'll repost here as a separate answer.

It looks like the SeismoSignal amplitude graph in this comment can be replicated with zero padding to nextpow2 and multiplying the DFT by Ts to approximate the Continuous Time Fourier Transform.

acc = readmatrix('EQ.txt');

Ts = 0.01; Fs = 1/Ts;

ACC = fft(acc,2^nextpow2(numel(acc)));

N = numel(ACC);

f = (0:N-1)/N*Fs;

figure;

semilogx(f,abs(ACC)*Ts),grid

xlim([0.1 10])

yticks(0:.05:.8)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

How can i get frequency domain of an earthquake?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

11 个评论
显示 9更早的评论隐藏 9更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How can i get frequency domain of an earthquake?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

11 个评论 显示 9更早的评论隐藏 9更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论