HOW TO FFT A NOISY EARTHQUAKE DATA

26 views (last 30 days)
If someone could help me, I need to FFT an actual earthquake data retrieved from the sensor installed in our school. I saw some articles that I first need to remove mean value, apply baseline correction, and filter before doing so. Thank you very much.

Accepted Answer

Star Strider
Star Strider on 6 Jun 2023
@Mary Claire — The problem with your data is that they are really noisy (that appears to be thermal or noise) that probably originates in the instrumentation. There does not appear to be any specific trend in the time-domain data, so detrending or highpass or bandpass filtering (to remove baseline drift, and optionally high-frequency noise) will not improve the result. The noise is broadband, and cannot be removed with a frequency-selective filter. A Savitzky-Golay filter occasionally works to eliminate broadband noise, however it does not with these data, at least when I have tried it.
T1 = readtable('RF68CENE4.6M05305mins.csv');
T1.Var2 = str2double(T1.Var2)
T1 = 30001×2 table
Var1 Var2 _______________________ ___________ 2023-05-30T15:20:03.996 -5.6046e-05 2023-05-30T15:20:04.006 6.8816e-05 2023-05-30T15:20:04.016 -6.9805e-05 2023-05-30T15:20:04.026 8.3372e-05 2023-05-30T15:20:04.036 -8.4646e-05 2023-05-30T15:20:04.046 4.4586e-05 2023-05-30T15:20:04.056 -1.0253e-05 2023-05-30T15:20:04.066 2.7837e-05 2023-05-30T15:20:04.076 -0.000125 2023-05-30T15:20:04.086 6.6923e-05 2023-05-30T15:20:04.096 -7.208e-05 2023-05-30T15:20:04.106 7.5452e-05 2023-05-30T15:20:04.116 -3.7701e-05 2023-05-30T15:20:04.126 -3.6386e-06 2023-05-30T15:20:04.136 2.1781e-05 2023-05-30T15:20:04.146 3.577e-05
figure
plot(T1.Var1, T1.Var2)
grid
figure
plot(T1.Var1, T1.Var2)
grid
xlim([T1.Var1(1) T1.Var1(500)])
% dt = [0; diff(T1.Var1)];
% dt.Format = 'hh:mm:ss.SSSSSS'
% dt_stats = [mean(dt);std(dt);mode(dt)]
Ts = 0.01;
Fs = 1/Ts;
Fn = Fs/2;
L = size(T1,1);
NFFT = 2^nextpow2(L)
NFFT = 32768
FTVar2 = fft((T1.Var2-mean(T1.Var2)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
envu = envelope(abs(FTVar2(Iv))*2, 100, 'peak');
figure
plot(Fv, abs(FTVar2(Iv))*2, 'DisplayName','Fourier Transform')
hold on
plot(Fv, envu, '-r', 'DisplayName','Upper Envelope')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (Units)')
legend('Location','best')
You need to examine the instrtumentation to determine where the noise is originating, and then do what you can to minimise it.
.
  5 Comments
Mary Claire
Mary Claire on 7 Jun 2023
@Star Strider I see, thank you very much for your help! My apologies for the trouble.
Star Strider
Star Strider on 7 Jun 2023
As always, my pleasure!
No apologies necessary!
.

Sign in to comment.

More Answers (2)

William Rose
William Rose on 6 Jun 2023
Edited: William Rose on 6 Jun 2023
[edit: add detrend command example; fix typos]
You are right that removing the mean value first is good.
x=x-mean(x);
will remove the mean value.
If the data shows a trend, then it would be reasonable to remove it. Instead of removing the mean only, you could do
x=detrend(x);
which removes the mean and the linear trend, if any.
As for pre-filtering, I would first look at results without pre-filtering, then decide if some pre-filtering is useful or justified.
There are various ways to compute the power spectrum, and various options when doing so. cpsd() is a pretty good place to start, in my opinion.
pxx=cpsd(x,x,...);
where the "..." represents the various extra arguments.
More later when I have more time.

Matt
Matt on 6 Jun 2023
Edited: Matt on 6 Jun 2023
You can directly compute your signal FFT as star strider did and you will see that something is happening around 5Hz.
An issue with a FFT spectrum is that you are making an harmonical analysis on the entire length of the signal and you lose some insight on when things are happing. For exemple if you record a guitar where you would play 3 different chords at 3 different moment, you would see 3 different tone on your FFT spectrum and you would not know which one was played first or second. The information is not lost but is in the phase of the FFT that you don't see doing a 2D plots.
A way to go aroud this issue is to slice the signal and compute a FFT for each piece, and plot those. The function spectrogram does this :
% uiopen('home/iscat/Downloads/RF68CENE4.6M05305mins.csv',1)
y = table2array(RF68CENE4(:,2));
% lengths of each slice ( so 100 pts=1s used to compute each fft line)
M = 1000;
% Nb of pts of overlap between slices (ie we should see a continious spectrum evolution)
L = 900;
% A weight on the slices to smooth the spectrum ? does not change much the
% result
g = bartlett(M);
Ndft = 2048;
Delta_T = table2array(RF68CENE4(30001,1))-table2array(RF68CENE4(1,1));
Delta_T = seconds(Delta_T/30000);
fs = 1/Delta_T; % instrument is running at 100Hz
[s,f,t] = spectrogram(y,g,L,Ndft,fs);
% plot results
subplot 211
imagesc(t,f,abs(s).^2)
xlabel('time [s]')
ylabel('FFT amplitude [Hz]')
title('Power')
subplot 212
imagesc(t,f,log10(abs(s).^2))
xlabel('time [s]')
ylabel('FFT amplitude [Hz]')
title('log10 Power')
sgtitle('Nice measure !')
You see that you have a vibrations around 5Hz from the beginnign until 200 s and a main movement around 180s.
Last remark, when chosing the length of the slice you are impacting the minimum frequency you will be able to detect (Shannon sampling theorem), so you need to be carefull with this.
  2 Comments
Mary Claire
Mary Claire on 7 Jun 2023
Thank you very much for your advise. However I'm having trouble excecuting the code. It seems to be having an error in line 51. I'm a novice in using Matlab, i'm sorry for the trouble. The file used was also attached.
Matt
Matt on 8 Jun 2023
Line 50 (in green) is what you get when you drag and drop your csv file in the console : matlab generates this command uiopen for you and that load the csv data. When you do this you get a variable ot type table in your workspace, called RF68CENE4 in your case.
Line 51 fails for you because you don't have the variable RF68CENE4 in your workspace. So you can either uncomment line 50 and change the path so that it leads to your csv data, or simply drag and drop the file in your console.
% type this :
clear
uiopen('path_to_data.csv',1); % or drag and drop
% then the code in my previous message

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!