Require the displacement and velocity dataset from the acceleration dataset

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

Answer 1

[Mathieu NOE]

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

[Md Armanul Hoda]

Thank you for your help, can you provide me some study material so that I can understand high pass. And I don’t clearly get that dc blocking filter part. If I well get literature, that will be good for my understanding.

Thank you in advance.

[Mathieu NOE]

hello again

dc blocking filter is indeed a high pass filter. it's the same thing, its purpose is to remove unwanted dc / low frequency content from a signal

there are tons of info about signal processing (including low, high and bandpass filters on the internet :

FOR DSP PRESENTATION (wustl.edu)

The Scientist and Engineer's Guide to Digital Signal Processing Introduction to Digital Filters (analog.com)

[Md Armanul Hoda]

Thank you for the documents it helps me a lot.

Is that the standard formula for P, an and b in dc blocking filter( because I haven’t find the DC filtering formulation in the book and I am curious to know how it works).

And our output as displacement and velocity should starts from the zero at t=0,

And displacement can be nearer about blue curve in the attached figure.

Please help me out in this.

Thanks in advance.

[Md Armanul Hoda]

[Mathieu NOE]

hello again

first , see function attached for more information about the DC blocking filter

it comes from this Fex submission : The DC Blocking Filter - File Exchange - MATLAB Central (mathworks.com)

[Mathieu NOE]

The blue curve is your "reference " displacement sensor ? can you share this data ?

it is good thing to have a calibrated reference sensor to test and fine tune the double integration process

also, it is recommended to verify the sensors sensivity and acquisition system gains to not introduce errors

[Mathieu NOE]

See below improved code

now you can force the initial conditions for velocity and displacement to be zero

function [vel,displ] = double_int2(x,fc,Fs,IC_flag)

    % IC_flag : if IC_flag > 0  initial condition are forced to 0 for velocity and displacement 

full code

  
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,1); % 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,IC_flag)

    % IC_flag : if IC_flag > 0  initial condition are forced to 0 for velocity and displacement 
    
    
    %% main code
    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)
    
    %% window (IC)
    window = ones(size(x));
    tmp = filter(b,a,window);
    window = window - tmp; % now we have a rising exponential (starts at 0 and finishes at 1) (this could also be done with step function)
    
    %%
    accel = filtfilt(b,a,x);                  % filter the acceleration data with filtfilt
        
    % accel to velocity integration
    vel = dt*cumtrapz(accel);           % integration 
    vel = filtfilt(b,a,vel);     % detrend data with filtfilt
     % apply window (IC flag)
    if IC_flag>0
        vel = vel.*window;
    end
   
    % velocity to displ integration
    displ = dt*cumtrapz(vel);           % integration
    displ = filtfilt(b,a,displ);     % detrend data with filtfilt
    % apply window (IC flag)
    if IC_flag>0
        displ = displ.*window;
    end

end

[Md Armanul Hoda]

I have sent you the data in your inbox.If possible,Please help me out in that.

[Mathieu NOE]

ok, I 'll have a look at it tomorrow

[Mathieu NOE]

hello again

with your reference data , I could slightly tune the cut off frequency of the dc blocking filter (not much : from 0.25 to 0.2 Hz)

the results are not bad IMHO, but maybe you expect a perfect match

I can't really say where the diffeence may come but maybe already each sensor has a different dynmic response , so it's maybe one reason for differences

code :

% reference data
load('u.mat')
u = u';
load('du.mat')
du = du';
channel = 2;
u = u(:,channel); % select here which channel to process
du = du(:,channel); % select here which channel to process


% double integration of accel data 
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
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.2;  % Hz , this cut off frequency must be below first significant peak frequency in your signal
[vel,displ] = double_int2(acc,fc,fs,1); % double integration (cumtrapz) with high pass filtering
figure(3)
plot(t,vel,'r',t,du,'b')
ylabel('Vel. [m/s]')
xlabel('Time [s]')
title('Estimated vs Real Velocity')
legend('Estimated','Real');

figure(4)
plot(t,displ,'r',t,u,'b')
ylabel('Dis. [m]')
xlabel('Time [s]')
title('Estimated vs Real Displacement')
legend('Estimated','Real');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [vel,displ] = double_int2(x,fc,Fs,IC_flag)

    % IC_flag : if IC_flag > 0  initial condition are forced to 0 for velocity and displacement 
    
    
    %% main code
    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)
    
    %% window (IC)
    window = ones(size(x));
    tmp = filter(b,a,window);
    window = window - tmp; % now we have a rising exponential (starts at 0 and finishes at 1) (this could also be done with step function)
    
    %%
    accel = filtfilt(b,a,x);                  % filter the acceleration data with filtfilt
        
    % accel to velocity integration
    vel = dt*cumtrapz(accel);           % integration 
    vel = filtfilt(b,a,vel);     % detrend data with filtfilt
     % apply window (IC flag)
    if IC_flag>0
        vel = vel.*window;
    end
   
    % velocity to displ integration
    displ = dt*cumtrapz(vel);           % integration
    displ = filtfilt(b,a,displ);     % detrend data with filtfilt
    % apply window (IC flag)
    if IC_flag>0
        displ = displ.*window;
    end

end

[Md Armanul Hoda]

Hi,

Sorry for late response thank you for your help!

In this you have taken channel as 2, in selection of acceleration the second channel acceleration is 3rd column of data.

I have provided you a raw displacement dataset. You are comparing it with the numerical finite element methods results displacement which is in u dataset.

Thank you, if possible please look into that.

[Md Armanul Hoda]

Or velocity mach can be okay with numerical one.

[Mathieu NOE]

hello

can you please clarify which datasets (files) must be compared and which channel do you want

I have those files :

data.mat

ddu.mat

du.mat

u.mat

NI_OBS.mat

[Md Armanul Hoda]

 I have extracted velocity from acceleration in this code with a cutoff frequency of 0.35. It somehow matches for the first 50 seconds, but afterward, a significant difference starts appearing. How can we overcome this? Do you have any ideas?
%%%%%
Data:
Experimental_acceleration in m/s^2
du: numerical velocity in m/s

%%%%%%%%code%%%%% clc;clear ;close all; % reference data % displacement load('u.mat') u = u'; load('du.mat') du = du'; channel = 1; u = u(1:18000,channel); % select here which channel to process du = du(:,channel); % select here which channel to process %% accel data load('Experimental_acceleration.mat') t = (0:17999)*0.005; acc = Experimental_acceleration(channel,1:18000)'; % select here which channel to process dt = mean(diff(t)); fs = 1/dt; % sampling rate

%% Acceleration plot figure(1) plot(t,acc) ylabel('Accel [m/s^2]') xlabel('Time [s]') title(' Acceleration') % xlim([0 80])

%% 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.3540; % Hz , this cut off frequency must be below first significant peak frequency in your signal [vel,displ] = double_int2(acc,fc,fs,1); % double integration (cumtrapz) with high pass filtering %% Velecity figure plot(t,vel,'b:','LineWidth',1.5) hold on plot(t,du,'k-','LineWidth',.5) xlabel('t(s)') ylabel('velocity (m/s)') ylim([-0.12,0.12]) % Create a second set of axes on the right side yyaxis right; hold on; plot(t,du-vel,'k:','LineWidth',1.5); ylabel('Error') ylim([-.01,.15]) legend('Highpass','Numerical','Error', 'Location', 'NorthEast','Orientation', 'horizontal') box on grid on; pbaspect([2.5 1 1]) % xlim([0 40]) set(gca, 'FontName', 'Times New Roman','FontSize',12,'YColor', 'k') title(strcat('Velocity @ storey no: ',num2str(channel))) % %% % save("du6.mat","vel") % %% Displacement % figure % plot(t,displ,'b:','LineWidth',1.5) % hold on % plot(t,u,'k-','LineWidth',.5) % xlabel('t(s)') % ylabel('displacement (m)') % ylim([-0.01,0.01]) % % Create a second set of axes on the right side % yyaxis right; % hold on; % plot(t,u-displ,'k:','LineWidth',1.5); % ylabel('Error') % ylim([-0.002,0.015]) % legend('estimated','Numerical','Error', 'Location', 'NorthEast','Orientation', 'horizontal') % box on % grid on; % pbaspect([2.5 1 1]) % xlim([0 80]) % set(gca, 'FontName', 'Times New Roman','FontSize',12,'YColor', 'k') % title(strcat('Displacement @ storey no: ',num2str(channel))) % save("u1_0.45.mat","displ") %%%%%%%%%%%%%%%%%%%%%%%% IC_flag : if IC_flag > 0 initial condition are forced to 0 for velocity and displacement %%%%%%%%%%%%%%% function [vel,displ] = double_int2(x,fc,Fs,IC_flag)

    %% main code
    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)

    %% window (IC)
    window = ones(size(x));
    tmp = filter(b,a,window);
    window = window - tmp; % now we have a rising exponential (starts at 0 and finishes at 1) 

    %%
    accel = filtfilt(b,a,x);                  % filter the acceleration data with filtfilt

    %% accel to velocity integration
    vel = dt*cumtrapz(accel);           % integration 
    vel = filtfilt(b,a,vel);            % detrend data with filtfilt
     % apply window (IC flag)
    if IC_flag>0
        vel = vel.*window;
    end

    %% velocity to displ integration
    displ = dt*cumtrapz(vel);           % integration
    displ = filtfilt(b,a,displ);        % detrend data with filtfilt
    % apply window (IC flag)
    if IC_flag>0
        displ = displ.*window;
    end

end

[Mathieu NOE]

hello again

I searched first to see if there was no (hidden) problem with the numerical integration process but could not find anything

then I noticed something strange : if I still understand right du is the result of a model and we compare that to integration of experimental acceleration

now , if we focus for example on t between 50.45 and t = 50.52 s we can see that the experimental acceleration has a sudden negative peak (and integrating a negative peak gives a decrease on the velocity and that is exactly what we can see on the bottom blue trace).

On the contrary, the black trace (computed velocity) has a positive velocity peak (increase) , and so for the first half of that velocity peak should correspond a positive acceleration, and the second half should correspond to a negative acceleration, which is exactly the contrary of what the experimental data shows

to better explain it , the next plot is the experimental acceleration plotted with the computed derivative of du we get this result and then the question is what explain the difference between experimental and computed acceleration ? I cannot answer that because I don't know how you do the experimental vs computation model in terms of inputs (are they exactly the same ?)

attached both pieces of code to generate the plots above

[Md Armanul Hoda]

Sorry to bother you again.

I got that problem at 50sec and removed that. That problem was with experimental data, which was stored in absolute coordinate. I subtracted the ground acceleration from the dataset to get relative acceleration and now result is coming better that earlier. I have a doubt:Is there any chance to get the error a minimum values as much possible.

[Mathieu NOE]

no problem - but why we have still some error remains a mystery to me

have you tried to change the high pass filter frequency ?

I am not saying that there are not other (better ?) methods to do acceleration integration , but I wonder also if the physical behaviour of your model could be slightly non linear . Can you do some distorsion (THD) analysis or make some FRF sweeps at different excitation levels ?

[Md Armanul Hoda]

I tried to compare model in frequency domain using FFT and first three frequency is matching well. I have compared separately fft of acceleration, And velocity fft by taking the processed data using high pass and numerical model. I have also compared data in time domain which was shared in previous text with high error.

I don’t know where I am lagging, I am attaching the data. If possible please help me out in this. Thank you very much in advance.

Data: u: numerical disp(m) du: numerical velocity(m/s) ddu: numerical acceleration(m/s^2) Relat_Exp_Acc: Experimental acc(m/s^2) Note: need dataset upto 18000 time step, not beyond that.

[Mathieu NOE]

can you share the last code you are using ?

[Md Armanul Hoda]

The code is same which you have shared with IC_flag. That one which I have attached in my last text. For different velocity I have used different cutoff frequencies, in between 0.1 to 0.2.

[Md Armanul Hoda]

should I send that code?

[Mathieu NOE]

yes please

[Md Armanul Hoda]

This is that code which I am using to get velocity data.

clc;clear ;close all;
% reference data
% displacement
load('u.mat')
u = u';
load('du.mat')
du = du';
channel = 2;
L1=18000;

u = u(1:L1,channel); % select here which channel to process
du = du(1:L1,channel); % select here which channel to process

%% accel data 
load('Relat_Exp_Acc.mat')
t = (0:L1-1)*0.005;
acc = Relat_Exp_Acc(channel,1:L1)'; % select here which channel to process
dt = mean(diff(t));
fs = 1/dt; % sampling rate 

%% Acceleration plot
figure(1)
plot(t,acc)
ylabel('Accel [m/s^2]')
xlabel('Time [s]')
title(' Acceleration')
% xlim([0 80])

%% 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.13;  % Hz , this cut off frequency must be below first significant peak frequency in your signal
[vel,displ] = double_int2(acc,fc,fs,1); % double integration (cumtrapz) with high pass filtering
%% Velecity
figure
plot(t,vel,'b:','LineWidth',1.5)
hold on
plot(t,du,'k-','LineWidth',.5)
xlabel('t(s)')
ylabel('velocity (m/s)')
ylim([-0.12,0.12])
% Create a second set of axes on the right side
yyaxis right;
hold on;
plot(t,du-vel,'k:','LineWidth',1.5);
ylabel('Error')
ylim([-.01,.15])
legend('Highpass','Numerical','Error', 'Location', 'NorthEast','Orientation', 'horizontal')
box on
grid on;
pbaspect([2.5 1 1])
xlim([0 90])
set(gca, 'FontName', 'Times New Roman','FontSize',12,'YColor', 'k')
title(strcat('Velocity @ storey no: ',num2str(channel)))
%%
%  save("du6_0.21.mat","vel")
%% Displacement
% figure
% plot(t,displ,'b:','LineWidth',1.5)
% hold on
% plot(t,u,'k-','LineWidth',.5)
% xlabel('t(s)')
% ylabel('displacement (m)')
% ylim([-0.01,0.01])
% % Create a second set of axes on the right side
% yyaxis right;
% hold on;
% plot(t,u-displ,'k:','LineWidth',1.5);
% ylabel('Error')
% ylim([-0.002,0.015])
% legend('estimated','Numerical','Error', 'Location', 'NorthEast','Orientation', 'horizontal')
% box on
% grid on;
% pbaspect([2.5 1 1])
% % xlim([0 80])
% set(gca, 'FontName', 'Times New Roman','FontSize',12,'YColor', 'k')
% title(strcat('Displacement @ storey no: ',num2str(channel)))
% save("u1_0.45.mat","displ")
%%%%%%%%%%%%%%%%%%%%%%%% IC_flag : if IC_flag > 0  initial condition are forced to 0 for velocity and displacement 
    %%%%%%%%%%%%%%%
function [vel,displ] = double_int2(x,fc,Fs,IC_flag)

    
    
    %% main code
    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)
    
    %% window (IC)
    window = ones(size(x));
    tmp = filter(b,a,window);
    window = window - tmp; % now we have a rising exponential (starts at 0 and finishes at 1) 
    
    %%
    accel = filtfilt(b,a,x);                  % filter the acceleration data with filtfilt
        
    %% accel to velocity integration
    vel = dt*cumtrapz(accel);           % integration 
    vel = filtfilt(b,a,vel);            % detrend data with filtfilt
     % apply window (IC flag)
    if IC_flag>0
        vel = vel.*window;
    end
   
    %% velocity to displ integration
    displ = dt*cumtrapz(vel);           % integration
%     displ = filtfilt(b,a,displ);        % detrend data with filtfilt
    % apply window (IC flag)
    if IC_flag>0
        displ = displ.*window;
    end

end

[Mathieu NOE]

hello

a tried to see if chnaging thehigh pass functions from filtfilt to filter would change anything : not much indeed

same if I apply or not the"window" to the acceleration before it get's integrated with cumtrapz

there is not much more I could think of right now to make both curves better match , sorry !!

here the code after some modifications, but IMO not really different results :

clc;clear ;close all;
% reference data
% displacement
load('u.mat')
u = u';
load('du.mat')
du = du';
channel = 2;
L1=18000;
u = u(1:L1,channel); % select here which channel to process
du = du(1:L1,channel); % select here which channel to process
%% accel data 
load('Relat_Exp_Acc.mat')
t = (0:L1-1)*0.005;
acc = Relat_Exp_Acc(channel,1:L1)'; % select here which channel to process
dt = mean(diff(t));
fs = 1/dt; % sampling rate 
%% Acceleration plot
figure(1)
plot(t,acc)
ylabel('Accel [m/s^2]')
xlabel('Time [s]')
title(' Acceleration')
% xlim([0 80])
%% 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.13;  % Hz , this cut off frequency must be below first significant peak frequency in your signal
[vel,displ] = double_int2(acc,fc,fs,1); % double integration (cumtrapz) with high pass filtering
%% Velecity
figure
plot(t,vel,'b:','LineWidth',1.5)
hold on
plot(t,du,'k-','LineWidth',.5)
xlabel('t(s)')
ylabel('velocity (m/s)')
ylim([-0.12,0.12])
% Create a second set of axes on the right side
yyaxis right;
hold on;
plot(t,du-vel,'k:','LineWidth',1.5);
ylabel('Error')
ylim([-.01,.15])
legend('Highpass','Numerical','Error', 'Location', 'NorthEast','Orientation', 'horizontal')
box on
grid on;
pbaspect([2.5 1 1])
xlim([0 90])
set(gca, 'FontName', 'Times New Roman','FontSize',12,'YColor', 'k')
title(strcat('Velocity @ storey no: ',num2str(channel)))
%%
%  save("du6_0.21.mat","vel")
%% Displacement
% figure
% plot(t,displ,'b:','LineWidth',1.5)
% hold on
% plot(t,u,'k-','LineWidth',.5)
% xlabel('t(s)')
% ylabel('displacement (m)')
% ylim([-0.01,0.01])
% % Create a second set of axes on the right side
% yyaxis right;
% hold on;
% plot(t,u-displ,'k:','LineWidth',1.5);
% ylabel('Error')
% ylim([-0.002,0.015])
% legend('estimated','Numerical','Error', 'Location', 'NorthEast','Orientation', 'horizontal')
% box on
% grid on;
% pbaspect([2.5 1 1])
% % xlim([0 80])
% set(gca, 'FontName', 'Times New Roman','FontSize',12,'YColor', 'k')
% title(strcat('Displacement @ storey no: ',num2str(channel)))
% save("u1_0.45.mat","displ")
%%%%%%%%%%%%%%%%%%%%%%%% IC_flag : if IC_flag > 0  initial condition are forced to 0 for velocity and displacement 
    %%%%%%%%%%%%%%%
function [vel,displ] = double_int2(accel,fc,Fs,IC_flag)
    
    
    %% main code
    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)
    
    %% window (IC)
    window = ones(size(accel));
    tmp = filter(b,a,window);
    window = window - tmp; % now we have a rising exponential (starts at 0 and finishes at 1) 
    
    %%
     % apply window (IC flag)
    if IC_flag>0
        accel = accel.*window;
    end
%       accel = filtfilt(b,a,accel);                  % filter the acceleration data with filtfilt
    accel = filter(b,a,accel);                  % filter the acceleration data with filter
         
    %% accel to velocity integration
    vel = dt*cumtrapz(accel);           % integration 
     % apply window (IC flag)
    if IC_flag>0
        vel = vel.*window;
    end
%     vel = filtfilt(b,a,vel);            % detrend data with filtfilt
    vel = filter(b,a,vel);            % detrend data with filter
   
    %% velocity to displ integration
    displ = dt*cumtrapz(vel);           % integration
    % apply window (IC flag)
    if IC_flag>0
        displ = displ.*window;
    end
%         displ = filtfilt(b,a,displ);        % detrend data with filtfilt
        displ = filter(b,a,displ);        % detrend data with filter

end

[Mathieu NOE]

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]

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