How can i get frequency domain of an earthquake?

조회 수: 20 (최근 30일)
ADNAN KIRAL
ADNAN KIRAL 2021년 1월 27일
댓글: Star Strider 2023년 6월 3일
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))

채택된 답변

Star Strider
Star Strider 2021년 1월 27일
편집: Star Strider 2021년 1월 27일
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
Star Strider
Star Strider 2021년 3월 3일
As always, my pleasure!
Star Strider
Star Strider 2023년 6월 3일
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
.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Paul
Paul 2023년 6월 3일
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)

카테고리

Help CenterFile Exchange에서 Vibration Analysis에 대해 자세히 알아보기

제품


릴리스

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by