HOW TO FFT A NOISY EARTHQUAKE DATA

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.

Answer 1

[Star Strider]

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

.

[William Rose]

That is really nice analysis and code by @Star Strider.

[Jade Ann Chloe]

Good day @Star Strider. I'm trying to remove the broadband noise in my data. Any suggestions how to code the Savitzky-Golay filter?

[Star Strider]

@Jade Ann Chloe — I generally use an order 3 polynomial and then choose the 'framelen' to produce the result I want. I generally begin with it about 25% the length of the signal, and then reduce it to no less than 1% of the length of the signal, stopping in the interim when I get the result I want. This is highly subjective and absolutely empirical, as is much of signal processing.

[Mary Claire]

@Star Strider I see, thank you very much for your help! My apologies for the trouble.

[Star Strider]

As always, my pleasure!

No apologies necessary!

.

Answer 2

[William Rose]

@Mary Claire,

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

Answer 3

[Matt]

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.

[Mary Claire]

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]

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