Effacer les filtres
Effacer les filtres

How can i get frequency domain of an earthquake?

18 vues (au cours des 30 derniers jours)
ADNAN KIRAL
ADNAN KIRAL le 27 Jan 2021
Commenté : Star Strider le 3 Juin 2023
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))

Réponse acceptée

Star Strider
Star Strider le 27 Jan 2021
Modifié(e) : Star Strider le 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 commentaires
Star Strider
Star Strider le 3 Mar 2021
As always, my pleasure!
Star Strider
Star Strider le 3 Juin 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.
I added the findpeaks and text calls out of interest
.

Connectez-vous pour commenter.

Plus de réponses (1)

Paul
Paul le 3 Juin 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)

Catégories

En savoir plus sur Vibration Analysis dans Help Center et File Exchange

Produits


Version

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by