How can i get frequency domain of an earthquake?
7 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
Star Strider
on 27 Jan 2021
Edited: Star Strider
on 27 Jan 2021
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 Comments
Star Strider
on 3 Jun 2023
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.
.
More Answers (1)
Paul
on 3 Jun 2023
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 Comments
See Also
Categories
Find more on Vibration Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!