HOW TO FFT A NOISY EARTHQUAKE DATA
26 views (last 30 days)
Show older comments
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.
0 Comments
Accepted Answer
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)
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)
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
More Answers (2)
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.
0 Comments
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
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
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!