HOW TO FFT A NOISY EARTHQUAKE DATA

27 次查看(过去 30 天)
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.

采纳的回答

Star Strider
Star Strider 2023-6-6
@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 个评论
Mary Claire
Mary Claire 2023-6-7
@Star Strider I see, thank you very much for your help! My apologies for the trouble.
Star Strider
Star Strider 2023-6-7
As always, my pleasure!
No apologies necessary!
.

请先登录,再进行评论。

更多回答(2 个)

William Rose
William Rose 2023-6-6
编辑:William Rose 2023-6-6
[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 2023-6-6
编辑:Matt 2023-6-6
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 个评论
Mary Claire
Mary Claire 2023-6-7
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 2023-6-8
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

请先登录,再进行评论。

Community Treasure Hunt

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

Start Hunting!

Translated by