Best baseline correction method for quasi-static response
조회 수: 25 (최근 30일)
이전 댓글 표시
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
Image Analyst
2023년 4월 7일
If you have any more questions, then attach your data and code to read it in with the paperclip icon after you read this:
답변 (1개)
Mathieu NOE
2023년 4월 6일
편집: Mathieu NOE
2023년 4월 6일
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
Mathieu NOE
2023년 4월 7일
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
참고 항목
카테고리
Help Center 및 File Exchange에서 Filter Analysis에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!