Require the displacement and velocity dataset from the acceleration dataset
8 views (last 30 days)
Show older comments
Md Armanul Hoda
on 6 Nov 2023
Commented: Mathieu NOE
on 20 Dec 2023
I have acceleration dataset and I want to calculate displacement and velocity dataset. I have used several algorithm which is already shared on this community , but I don’t get the satisfactory result.can anyone please help me out.This is the link of dataset: https://drive.google.com/file/d/19n07MOFw2bz-yn3dISXOTf-zf9QpktNZ/view?usp=sharing
0 Comments
Accepted Answer
Mathieu NOE
on 6 Nov 2023
Edited: Mathieu NOE
on 6 Nov 2023
hello
try this
I assumed your data are in g's maybe this must be changed according to your sensor sensivity / acquisition steup
the main issue in accelaration signal integration is to avoid drift caused by DC compoent integration
so detrending / high pass filtering is most of the time needed, but the HP filyer cut off frequency must be carefully chosen (below the main peak of the accel signal FFT spectrum)
as your first peak appear close to 1.5 Hz , I opted for a HPF cut off frequency = 0.25 Hz.
filename = 'data.mat';
data = load(filename);
t = data.data(:,1);
acc = data.data(:,2); % select here which channel to process
dt = mean(diff(t));
fs = 1/dt; % sampling rate
% Acceleration data
acc = acc*9.81; % gravity factor ( g to m/s²)
figure(1)
plot(t,acc)
ylabel('Accel [g]')
xlabel('Time [s]')
title(' 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;
f = fs/2*linspace(0,1,NFFT/2+1);
accSpectrum = abs(accSpectrum(1:NFFT/2+1));
figure(2);
loglog(f, accSpectrum);
title('FFT of acceleration')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
hold on
% add findpeaks to display dominant frequencies
NP = 6;
[PKS,LOCS] = findpeaks(accSpectrum,f,'SortStr','descend','NPeaks',NP);
plot(LOCS,PKS,'dr');
for ck = 1:NP
text(LOCS(ck)*1.1,PKS(ck)*1.2,[num2str(LOCS(ck),3) 'Hz']);
end
hold off
% main code
fc = 0.25; % Hz , this cut off frequency must be below first significant peak frequency in your signal
[vel,displ] = double_int2(acc,fc,fs); % double integration (cumtrapz) with high pass filtering
figure(3)
plot(t,vel)
ylabel('Vel. [m/s]')
xlabel('Time [s]')
title('Estimated Velocity')
figure(4)
plot(t,displ)
ylabel('Dis. [m]')
xlabel('Time [s]')
title('Estimated Displacement')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [vel,displ] = 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)
accel = filtfilt(b,a,x); % filter the test data with filtfilt
% accel to velocity integration
vel = dt*cumtrapz(accel); % integration
vel = filtfilt(b,a,vel); % detrend data with filtfilt
% velocity to displ integration
displ = dt*cumtrapz(vel); % integration
displ = filtfilt(b,a,displ); % detrend data with filtfilt
end
26 Comments
Mathieu NOE
on 20 Dec 2023
hello again
so I wanted to show how the different integration filters can be compared
the simple code below shows the different option
- the black curve is the Bode plot of an ideal integrator (H(s) = 1/s)
- the light blue (cyan) curve is the Bode plot of the cumtrapz equivalent IIR filter
- the red curve is the Bode plot of the cumtrapz equivalent IIR filter followed by the DC blocking filter (1st order)
- the blue curve is the Bode plot of an pseudo integrator ( H(s) = 1/w0 * s / (s^2 + s + 1)) that we want to compare to the cumtrapz / DC blocking filter combination
you see that they are pretty similar above fc frequency , only below fc there is a difference :
- the cumtrapz / DC blocking filter combination will still have a non zero dc gain
- while the pseudo integrator (PIF in short) will have a 1st order high pass behaviour (so no gain at dc frequency)
code
%%%%%%%%%%%%%%%%%%%% inputs %%%%%%%%%%%%%%%
Fs = 1e2;
fc = 0.1;
%%%%%%%%%%%%%%% simulation %%%%%%%%%%%%%%%
freq = logspace(-2,(log10(Fs/2.56)),100);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PIF: H(s) = 1/w0 * s / (s^2 + s + 1) (pseudo integrator filter with built in high pass effect)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w0 = 2*pi*fc;
w0r = w0/Fs;
alpha = sin(w0r)/2;
b0 = alpha;
b1 = 0;
b2 = -alpha;
a0 = 1 + alpha;
a1 = -2*cos(w0r);
a2 = 1 - alpha;
num1z=[b0 b1 b2]/w0; % the corrective gain factor /w0 is to have same gain as true integrator
den1z=[a0 a1 a2];
[g1,p1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
g1db = 20*log10(g1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cumptrapz method equivalent IIR filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num2z= 0.5*[1 1]/Fs; % numerator
den2z= [1 -1]; % numerator
[g2,p2] = dbode(num2z,den2z,1/Fs,2*pi*freq);
g2db = 20*log10(g2);
% 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)
[g3,p3] = dbode(b,a,1/Fs,2*pi*freq);
g3db = 20*log10(g3);
%% plot
figure(1);
subplot(2,1,1),semilogx(freq,-20*log10(2*pi*freq),'-*k',freq,g1db,'b',freq,g2db,'c',freq,g2db+g3db,'r');
title(' PIF: H(s) = 1/w0 * s / (s^2 + s + 1)');
ylabel('dB ');
legend('ideal integrator','PIF pseudo integrator filter','cumptrapz filter','cumptrapz + dc blocking filter')
subplot(2,1,2),semilogx(freq,-90*ones(size(p1)),'-*k',freq,p1,'b',freq,p2,'c',freq,p2+p3,'r');
xlabel('Frequency (Hz)');
ylabel(' phase (°)')
Mathieu NOE
on 20 Dec 2023
So I tested the new filter (PIF) vs the cumtrapz + dc blocking , but at the end the results are very similar.
The PIF filter is not enough (unfortunately) to reduce the error between experimental and numerical data.
could it be that the numerical model must be corrected (have you made some king of numerical model validation and tuning ?)
attached are the two codes FYI
all the best
More Answers (0)
See Also
Categories
Find more on Vibration Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!