How can I deduce velocity and displacement trends from an acceleration signals

Hello everyone, I have an acceleration signal from an accelerometer. I wrote a simple script in order to analyse this signal: I plotted the signal in time domain, I performed a frequency domain analysis, I filtered and realigned the signal and at the end I tried to reconstruct velocity and displacement with the cumtrapz function. The code is in the following:

clc;
clear all;
close all;

file_acc='CT4H1X21-2.xlsx';

x_acc=xlsread(file_acc);

%%Time domain analysis of acceleration signal

Fsa=19200;
Ta=1/Fsa;

xa=x_acc(:,3);
Na=length(xa);
ta=linspace(0,Na*Ta,Na);

figure
plot(ta,xa)
grid on
title('Acc. on the sparger')
xlabel('t (s)')
ylabel('xa(t), Acceleration (g)')

%%Frequency domain analysis of accelerometric signal

NFFT2_acc = 2^nextpow2(Na);          % Next power of 2 from length of y
Xa=fft(xa-mean(xa),NFFT2_acc)/Na;
fa=Fsa/2*linspace(0,1,NFFT2_acc/2+1);

figure;
plot(fa,2*abs(Xa(1:NFFT2_acc/2+1))) 
grid on
title('Single-Sided Amplitude Spectrum of xa(t)')
xlabel('f (Hz)')
ylabel('Xa(f)')

%%Plot the Power Spectral Density of the accelerometric signal

figure;
[pxxa,fxa] = pwelch((xa-mean(xa)),[],[],[],Fsa);
plot(fxa,pxxa);
grid on
title('Power Spectral Density of xa(t)')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

%%Butterworth filter for accelerometric signal

Fca=150;
Fsa=19200;

[ba,aa]=butter(10,Fca/(Fsa/2));

figure;
freqz(ba,aa)
dataIna=xa;

% Apply the filter to Smooth out the Signal 

xafilter=filter(ba,aa,dataIna);

% Overlay the filtered pressure signals on the original signals
% Filtered signal is delayed

figure;
plot(ta,xa,'b',ta,xafilter,'r');
grid on;
title('xa (accelerometric signal)')
xlabel('t (s)')
ylabel('acc. (g)')
legend({'Original xa Signal','Filtered xa Signal'});

%set(gcf,'NumberTitle','Off', 'Name','Filtered Signal vs. Actual Signal');

%%Compare the original accelerometric signal and delay compensated filtered accelerometric signal

figure;
%xfiltfilt = filtfilt(d.sosMatrix,d.ScaleValues,x2);
xfiltfilta=filtfilt(ba,aa,xa);
plot(ta,xa,ta,xfiltfilta);
grid on
title('xa accelerometric signal')
xlabel('t (s)')
ylabel('acc (g)')
legend({'Original xa Signal','Actual xa Signal (filtered and realigned signal)'});

%%Frequency domain analysis of the filtered and realigned accelerometric signal

Xaf=fft(xfiltfilta-mean(xfiltfilta),NFFT2_acc)/Na;
figure;
plot(fa,2*abs(Xaf(1:NFFT2_acc/2+1))) 
grid on
xlim([0 150])
title('Single-Sided Amplitude Spectrum of xa(t) filtered and realigned Signal')
xlabel('f (Hz)')
ylabel('Xa(f)')

%%Velocity

acceleration=xfiltfilta*9.80665;
velocity=cumtrapz(ta,acceleration);
figure;
plot(ta,velocity);
grid on
title('Velocity of the sparger from the acceleration signal')
xlabel('t (s)')
ylabel('Velocity (m/s)')

%%Displacement

displacement=cumtrapz(ta,velocity);
figure;
plot(ta,displacement);
grid on
title('Displacement of the sparger from the velocity signal')
xlabel('t (s)')
ylabel('Displacement (m)')

I'm not sure of what I did; in particular I found displacement of the order of meters and this is clearly an error. May someone take a look at my code? Thank you very much.

Answer 1

[Star Strider]

I am not sure what you did, either.

Without your data, it’s not possible definitively to determine that.

One problem:

xafilter=filter(ba,aa,dataIna);

Use filtfilt! You used it in other parts of your code. It has a maximally flat phase response, so there will be no delays.

Always check your filter performance with the freqz function so that you know that it does what you want it to.

As for the high displacement, you may need to redesign your filter to be a bandpass filter rather than the default lowpass filter. With a bandpass filter, you can completely remove the entire d-c (constant) offset and low-frequency baseline drift. Subtracting the mean as you did here will not do that, because it only removes the d-c offset, not low-frequency baseline drift. The low-frequency baseline drift will also integrate, and could be the reason your displacements are in parsecs rather than microns.

[Guglielmo Giambartolomei]

Thak you Star Strider. I used filtfilt. Could u suggest me a bandpass filter in order to avoid the low-frequency contents? Thank you.

[Star Strider]

You will have to decide on the appropriate frequencies, but a bandpass filter design using the kaiserord and fir1 functions that works is here:

Fs = 44100;                                                     % Sampling frequency
fcuts = [10 20 20E+3 21E+3];                                    % Frequency Vector
mags =   [0 1 0];                                               % Magnitude (Defines Passbands & Stopbands)
devs = [0.05  0.01  0.05];                                      % Allowable Deviations
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs);
n = n + rem(n,2);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'scale');

figure(1)
freqz(hh, 1, 2^14, Fs)

Note that the values in the ‘fcuts’ vector are in Hz, so you only have to change those and the value for ‘Fs’. Do not normalise them by the Nyquist frequency in your code. The kaiserord function produces the normalised frequency vector that fir1 requires.

[Guglielmo Giambartolomei]

You are so kind Star Strider. Thank You very much!

[Star Strider]

My pleasure!