Best baseline correction method for quasi-static response
24 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi, im trying to create a code for removing drift in a peice of acceleration signal i have. The signal measure a quasi-static response and im integrating it twice to get displacment but of course there are issues with low frequency noise and drift. If there's any advice on what baseline correction methods would suit best id appreciate it. filtering the signal removes some of the displacement so isnt ideal for this type of signal.
1 Kommentar
Image Analyst
am 7 Apr. 2023
If you have any more questions, then attach your data and code to read it in with the paperclip icon after you read this:
Antworten (1)
Mathieu NOE
am 6 Apr. 2023
Bearbeitet: Mathieu NOE
am 6 Apr. 2023
hello
this seems to be a reccurent topic on this forum... so here we go again...
the best advice I can give are :
1/ make sure your sensor is correctly calibrated and correct for offset and sensivity error before you do any integration
2/ look at your acceleration signal spectrum (make a FFT of it) and look at the spectrum shape.
as you will probably do some high pass filtering to remove offset and drift , you have to choose the best corner frequency of your filter :
- not too low or you will integrate some portion of the drift / low freq noise
- not too high or your double integration will give an underestimated displacement
One possible code is to use cumtrapz for the integration followed by a dc wash out filter.
repeat that combo twice for double integration
I use filtfilt for forward and backward filtering so no phase distorsion (if that matters to you)
data file attached and code example below :
data = readmatrix('GEE Earthquake.txt'); % time , accel (g)
t = data(:,1);
dt = mean(diff(t)); % Average dt
fs = 1/dt; % Frequency [Hz] or sampling rate
% Acceleration data
acc = data(:,2)*9.81; % gravity factor ( g to m/s²)
figure(1)
plot(t,acc)
ylabel('Accel [g]')
xlabel('Time [s]')
title('Earthquake Acceleration')
% fft of the Acceleration signal
L = length(acc); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
accSpectrum = fft(acc, NFFT)/L; % detrend the signal before fft to remove the large DC (and near DC) fft output
f = fs/2*linspace(0,1,NFFT/2+1);
figure(2);
loglog(f, abs(accSpectrum(1:NFFT/2+1)))
title('FFT of acceleration')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
% main code
fc = 0.25; % Hz
displ = MN_double_int2(acc,fc,fs); % double integration (cumtrapz) with high pass filtering
figure(3)
plot(t,displ)
ylabel('Dis. [m]')
xlabel('Time [s]')
title('Estimated Displacement')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ydintd] = MN_double_int2(x,fc,Fs)
dt = 1/Fs;
% DC blocking filter (discrete)
fn = fc/(Fs/2); % normalized cutt off frequency
p = (sqrt(3) - 2*sin(pi*fn))/(sin(pi*fn) + sqrt(3)*cos(pi*fn));
b = [1 -1]; % (numerator)
a = [1 -p]; % (denominator)
y = filtfilt(b,a,x); % filter the test data with filtfilt
% accel to velocity integration
yint = dt*cumtrapz(y); % integration
yintd = filtfilt(b,a,yint); % detrend data with filtfilt
% velocity to displ integration
ydint = dt*cumtrapz(yintd); % integration
ydintd = filtfilt(b,a,ydint); % detrend data with filtfilt
end
2 Kommentare
Mathieu NOE
am 7 Apr. 2023
hello
what you call baseline correction does exist but IMHO isnot suited to accelerometer signals
this will do a baseline correction but on data that are only positive values not an AC signal
Siehe auch
Kategorien
Mehr zu Filter Analysis finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!