How can i get frequency domain of an earthquake?

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))

Answer 1

[Star Strider]

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.

[ADNAN KIRAL]

Thanks for the reply, but just wondering why I got a different plot when I use SeismoSignal software?

please see the image attached here.

thanks in advance for the help.

[Star Strider]

You got a different result because SeismoSignal plots the frequency axis in logarithmic units. My plots are in linear units. The SeismoSignal softtware may also have increased the frequency resolution of the fft result by zero-padding it. That is straightforward to do in MATLAB, however there appeared to be no need to do it with your signals, so I did not.

Change the last 2 plots in my code to:

figure
semilogx(Fv, abs(FTacc(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Acceleration Amplitude (Units)')
xlim([0  15])

figure
semilogx(Fv, (abs(FTacc(Iv))*2).^2)
grid
xlabel('Frequency (Hz)')
ylabel('Acceleration Power (Units^2)')
xlim([0  15])

and they will be the same as the SeismoSignal plots.

I only plotted the amplitude and power spectra, and did not zero-pad the fft. Plotting the zero-padded fft result and the phase as well is straightforward, if you want to do that.

[ADNAN KIRAL]

thanks again, but the amplitude is NOT the same as SeismoSignal. why is that ? sorry about disturbing.

matlab= 0.08 ;

SeismoSignal= 0.75;

thanks for reply.

[Star Strider]

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.)

[ADNAN KIRAL]

i see, many thanks for kind reply.

[Star Strider]

As always, my pleasure!

[ADNAN KIRAL]

Dear @Star Strider

My apologies for disturbing you again. I have received an email from SeismoSignal for the calculation of FFT in their program.

may I ask what is the difference between this and MATLAB calculation? can I modify the above code to get the same result? many thanks for the feedback. any suggestion?

the email:

Seismosoft Sales <[e-mail]>

to me:

The FFT Amplitude plotted in Seismosignal is the SQRT(Re^2 + Im^2) where Re and Im are the Real and Imaginary parts of the Fast Fourier Transform of the input signal. As a result it should not be expected equal to the oscillation amplitude of your signal.

In general (especially for noisy signals like earthquake motions) Fourier Amplitude gives more a sense of the energy included in each frequency component of the signal rather than the real oscillation amplitude of this frequency component.

thanks in advance...

[Star Strider]

I agree with the SeismoSignal e-mail.

My code apparently duplicates the SeismoSignal result, at least in doing the same essential presentation. There could be differences in the code SeismoSignal uses to calculate its version of the fast Fourier transform and the one that MATLAB uses. I have no idea how to determine that.

[ADNAN KIRAL]

I see, thanks again.

[Star Strider]

As always, my pleasure!

[Star Strider]

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]

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)