How can i get frequency domain of an earthquake?
9 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
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
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.
.
更多回答(1 个)
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)
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Vibration Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!