Sampling and resampling data of different lengths

Hi,

I have two sets of data, one from a motion analysis software (Vicon) which is captured at 100Hz and another set of data from a IMU called motion sense. The trouble I am having is when I plot the data against each other for a comparison study they do not match up. The motionsense software is sampled at about 50Hz, however, there are slight variations with in the data of 0.02 second. Making the sample rate 49.8 at times. I want both the vicon and motionsense to start and finish at the same time, but showing all the data.

So if they are at different frequencies I obvsiouly can't just cut them after 100 data points as I lose alot of data.

Currently if I say I want 0.5 seconds of data to display for each then I will have 2 lines on a graph of different lengths.

Is there a way to get the sample rates the same and constant within the file without losing any data as well as ensuring the data sets start and finish at the same time.

I am trying to get the 50Hz sampling constant for the motionsense data set so that I know I can match both graphs up.

The 67CT_log is the data from the motion sense- The time stamps are in the spreadsheet fore reference.

Vicon has a time stamp for each cell starting from 1,2,3,4,ect,ect

The graph I am getting out is as seen in the attachment.

Both graphs should link up as they are recording the same movements at the same time just on different devices.

Any help would be great.

Thank you

[Mathieu NOE]

hello Alexandra

do we really have a problem here ? I thought first that the major problem was plotting two sets of data with one being not truly sampled at fixed rate (motionsense ) , but this is not a plotting issue as if you plot two sets of data together, one can be sampled at one fixed rate while the second can have irregular time stamps and both lines can perfectly overlay ; it's just that the time stamps (either regular or not must) truly match acquisition sampling time stamps;

now looking at the two curves in your figure it seems to see that they are simply time shifted about 0.5 s (looking at the time delta of positive and negative peaks , this is everywhere the same delta) but this may be due to the fact that the two records starts / stops are not synchronized - time duration are also different

beside that, I don't understand how to interpret the time stamps in the csv file as the delta between two consecutive values is either 658 or 659 , in which units ? relation to 50 Hz rate is not obvious to me

all the best

[alexandra ligeti]

Hi,

Thank you for your reply.

The timestamps delta shifts between 658-659 as you pointed out. It is related to 50Hz as 658/2^15 =0.02 sec, which is slighly under 50Hz. I was wondering how I account for these fluctuations? I would like each sample to be at a constant of 50Hz and for the data to completely line up. Is there a way to Timeshift everything and make the sampling rate consistent within the file?

Thank you once again.

[Mathieu NOE]

hello Alexandra

do you have a code (+data) on which I could work to help you ?

tx

[alexandra ligeti]

Hi there,

So those two files link up.The Dynamic_3M is from Vicon which operates at 100Hz and the 67CT imu sensor data operates at about 50Hz however this varies within the file.

My code currently is as shown below. I want to have it time constant so that I can draw up direct comparison between the data to test whether this sensor is accurate to the gold standard vicon data, but I can because the points in time are different so I cannot compare 0.01 seconds from 1 file to the other because of the differeneces in sampling frequency and the variation in the time stamps of the 67CT motion sense IMU file.

Please let me know if I need to provide anything else, or clarify anything.

Thank you

% Comparing Vicon to MotionSense
% Only 1 dynamic trial carried out


clear all
close all
%% Open Desired Files
%Open Static 3 Marker files 
Vicon3Static = readmatrix('C:\Users\lexil\Documents\PhD\Code\ViconData\20210614\Static_3M_Sensor.csv');

%Open Dynamic 3 Marker files 
Vicon3Dynamic = readmatrix('C:\Users\lexil\Documents\PhD\Code\ViconData\20210614\Dynamic_3M_Sensor01.csv');

%Open motion sense file
Sensor_Static = readmatrix('C:\Users\lexil\Documents\PhD\Code\MotionSenseData\14_06_2021\67CT_log_14062021_12-19-01.csv');
Sensor_Dynamic = readmatrix('C:\Users\lexil\Documents\PhD\Code\MotionSenseData\14_06_2021\67CT_log_14062021_12-26-20.csv');

%% Marker Points
% Marker Positions for static marker 
marker_1_S = Vicon3Static(:,3:5)/1000;
marker_2_S = Vicon3Static(:,6:8)/1000;
marker_3_S = Vicon3Static(:,9:11)/1000;

% Marker Positions for Dynamic marker 
marker_1_D = Vicon3Dynamic(:,3:5)/1000;
marker_2_D = Vicon3Dynamic(:,6:8)/1000;
marker_3_D = Vicon3Dynamic(:,9:11)/1000;


%% Create Vectors 
%Static Vectors
Static_vector_1 = marker_2_S - marker_1_S;
Static_vector_2 = marker_3_S - marker_1_S;
Static_vector_3 = cross(Static_vector_1, Static_vector_2);
Static_vector_4 = cross(Static_vector_3, Static_vector_1); %Vector in pendulums swining plane


%Dynamic Vectors
Dynamic_vec_1 = marker_2_D - marker_1_D;
Dynamic_vec_2 = marker_3_D - marker_1_D;
Dynamic_vec_3 = cross(Dynamic_vec_1, Dynamic_vec_2);
Dynamic_vec_4 = cross(Dynamic_vec_3, Dynamic_vec_1); %Vector in pendulums swining plane


%% Correcting sizes of Vectors
%Determine length of vectors
sv4 = size(Static_vector_4);
sv4D = size(Dynamic_vec_4);

%Determine longest vector
pmax = max([size(Static_vector_4,1) size(Dynamic_vec_4,1)]);

%Padding with zero to make same length
Zpad = zeros(pmax,3);
NewStatic_vector_4 = Zpad;
NewDynamic_vec_4 = Zpad;

NewStatic_vector_4(1:size(Static_vector_4),:)= Static_vector_4;
NewDynamic_vec_4(1:size(Dynamic_vec_4),:)= Dynamic_vec_4;


%% Vector Magnitudes for pendulum plane
%Static Vectors
Mag_Static_4 = sqrt((NewStatic_vector_4(:,1)).^2 + (NewStatic_vector_4(:,2)).^2 + (NewStatic_vector_4(:,3)).^2);

%Dynamic Vectors
Mag_Dynamic_vec_4 = sqrt((NewDynamic_vec_4(:,1)).^2 + (NewDynamic_vec_4(:,2)).^2 + (NewDynamic_vec_4(:,3)).^2);
 
%% Unit vectors for the static and dynamic vector in plane of pendulum swing
%Static Vectors
unit_Static_4 = NewStatic_vector_4./(Mag_Static_4);


%Dynamic Vectors
%Small Swing
unit_Dynamic_vec_4 = NewDynamic_vec_4./(Mag_Dynamic_vec_4);


%% Find angle between static reference vector and changing dynamic vector
% Angle between static and slow dynamic
%Sine Rule
SinTheta_1 = cross(unit_Static_4,unit_Dynamic_vec_4);
Theta_1 = 180.*(asin(SinTheta_1))./pi;

%% MotionSense angle
%MotionsenseAngle = MotionSenseDynamic(1:248,6);
MotionsenseAngle = Sensor_Dynamic(:,6);


%% Synchronisation
fs_vicon = 100;
fs_sensor = 50; %approximately 50HZ

dt_v = 1/fs_vicon;
dt_s = 1/fs_sensor;

%Added 18/06/2021- might not be correct way
dt_v_Resample = resample(dt_v,1,1);
%dt_s_Resample = resample(dt_s,1,2); 

%Phil code- does work
% time_v = 0:dt_v:549*dt_v;
% time_v = time_v';
% 
% time_s = 0:dt_s:273*dt_s;
% time_s = time_s';

%Added 18/06/2021- slightly addjusted to account for addition of resampling
time_v = (0:dt_v_Resample:549*dt_v_Resample);%+0.54;
time_v = time_v';

time_s = 0:dt_s:273*dt_s;
time_s = time_s';


%% Offset
x_offset = time_s - time_v(1:274,1);
y_offset = Theta_1(1:274,1) - MotionsenseAngle(:,1);

%% Plot vicon angle against time
figure
x = time_v+0.54;   % Recorded at 47HZ
y = Theta_1(:,1);
plot(x,y) %blue curve 
title('Angle comparison between Vicon and MotionSense')
xlabel('Time in sec')
ylabel('Angle in degrees')
hold on

% x = time_v;
% y = Theta_1(:,2);
% plot(x,y)
% hold on 
% 
% x = time_v;
% y = Theta_1(:,3);
% plot(x,y)
% hold on 

x = time_s;
y = MotionsenseAngle;
plot(x,y);

legend('Vicon','MotionSense');

%% Calculate Error between the two systems
%root mean square error (RMSE)and Bland-Altman plots
Error = Theta_1(1:274,1) - MotionsenseAngle; %Error between Vicon and sensor
SqError = Error.^2;                          %Square error
AvSqError = mean(SqError);                   %Mean Squared Error  
RMSE = sqrt(AvSqError);                      %Root Mean Squared Error

%Pearson’s R correlation 
R = corrcoef(Theta_1(1:274,1),MotionsenseAngle);

%Bland_Altman plot

[Mathieu NOE]

hello Alexandra

seems to me that I need also the "static" measurements files - according to what I undertstand so far in your code :

%% Find angle between static reference vector and changing dynamic vector
% Angle between static and slow dynamic
%Sine Rule
SinTheta_1 = cross(unit_Static_4,unit_Dynamic_vec_4);
Theta_1 = 180.*(asin(SinTheta_1))./pi;

[alexandra ligeti]

Hi Mathieu,

So sorry thought I attached all the files... Here they are. Again the IMU sensor created the 67CT file and vicon has the Static_3M file.

Thank you so much for your troubles.

Answer 1

[Mathieu NOE]

hello Alexandra

this was a good gym exercise for me today. had to spend first time to understand which variables belong to which sensor (i made some rewording afterwards to make the difference clear to me).

then a bunch of resampling and time offset latter , we're good.

the error can be computed as both datas are now related to the same common time axis (322 samples)

hope this helps

the code will display figures that illustrate how we progress in the code

also it uses a subfunction (see attachement)

% Comparing Vicon (reference sensor / data) to MotionSense (other sensor / data)
% Only 1 dynamic trial carried out
clearvars
% close all

%% Open Desired Files
% %Open Static 3 Marker files 
% Vicon3Static = readmatrix('C:\Users\lexil\Documents\PhD\Code\ViconData\20210614\Static_3M_Sensor.csv');
Vicon3Static = readmatrix('Static_3M_Sensor.csv');  % (reference sensor / data)
[mViconS,nViconS] = size(Vicon3Static);
% %Open Dynamic 3 Marker files 
% Vicon3Dynamic = readmatrix('C:\Users\lexil\Documents\PhD\Code\ViconData\20210614\Dynamic_3M_Sensor01.csv');
Vicon3Dynamic = readmatrix('Dynamic_3M_Sensor01.csv');  % (reference sensor / data)
[mViconD,nViconD] = size(Vicon3Dynamic);

    % Vicon time vector 
    fs_vicon = 100;
    dt_v = 1/fs_vicon;
    time_v =  Vicon3Dynamic(:,1)*dt_v;
    time_v =  time_v - time_v(1); % so it starts at zero

% %Open motion sensor file
% Sensor_Static = readmatrix('C:\Users\lexil\Documents\PhD\Code\MotionSenseData\14_06_2021\67CT_log_14062021_12-19-01.csv');
% Sensor_Dynamic = readmatrix('C:\Users\lexil\Documents\PhD\Code\MotionSenseData\14_06_2021\67CT_log_14062021_12-26-20.csv');
MS_Sensor_Static = readmatrix('67CT_log_14062021_12-19-01.csv');   % (other sensor / data)
MS_Sensor_Dynamic = readmatrix('67CT_log_14062021_12-26-20.csv');   % (other sensor / data)


    % now resample the Dynamic data according to "reference" time vector = time_v
    % motion sensor time vector  (approx 50 Hz)
    dt_s = 1/2^15;
    time_s =  MS_Sensor_Dynamic(:,1)*dt_s;
    time_s =  time_s - time_s(1); % so it starts at zero
    MS_Sensor_Dynamic_new = interp1(time_s,MS_Sensor_Dynamic,time_v,'linear');
    
%     new_time_s =  NewDynamic_vec_4(:,1)*dt_v; % (check / debug only)
%     new_time_s =  new_time_s - new_time_s(1); % new_time_s = time_v now , so this proves interpolation works fine (check / debug only)

%% Marker Points (Vicon data)
% % Marker Positions for static marker 
marker_1_S = Vicon3Static(:,3:5)/1000;
marker_2_S = Vicon3Static(:,6:8)/1000;
marker_3_S = Vicon3Static(:,9:11)/1000;
% Marker Positions for Dynamic marker 
marker_1_D = Vicon3Dynamic(:,3:5)/1000;
marker_2_D = Vicon3Dynamic(:,6:8)/1000;
marker_3_D = Vicon3Dynamic(:,9:11)/1000;
%% Create Vectors (Vicon data)
% %Static Vectors
Static_vector_1 = marker_2_S - marker_1_S;
Static_vector_2 = marker_3_S - marker_1_S;
Static_vector_3 = cross(Static_vector_1, Static_vector_2);
Static_vector_4 = cross(Static_vector_3, Static_vector_1); %Vector in pendulums swining plane
%Dynamic Vectors
Dynamic_vec_1 = marker_2_D - marker_1_D;
Dynamic_vec_2 = marker_3_D - marker_1_D;
Dynamic_vec_3 = cross(Dynamic_vec_1, Dynamic_vec_2);
Dynamic_vec_4 = cross(Dynamic_vec_3, Dynamic_vec_1); %Vector in pendulums swining plane


%% Correcting sizes of Vectors
%Determine length of vectors
sv4 = size(Static_vector_4);
sv4D = size(Dynamic_vec_4);
%Determine longest vector
pmax = max([size(Static_vector_4,1) size(Dynamic_vec_4,1)]);
% %Padding with zero to make same length
Zpad = zeros(pmax,3);
NewStatic_vector_4 = Zpad;
NewDynamic_vec_4 = Zpad;
NewStatic_vector_4(1:size(Static_vector_4),:)= Static_vector_4;
NewDynamic_vec_4(1:size(Dynamic_vec_4),:)= Dynamic_vec_4;

%% Vector Magnitudes for pendulum plane
%Static Vectors
Mag_Static_4 = sqrt((Static_vector_4(:,1)).^2 + (Static_vector_4(:,2)).^2 + (Static_vector_4(:,3)).^2);
Mag_Static_4 = mean(Mag_Static_4);

%Dynamic Vectors
Mag_Dynamic_vec_4 = sqrt((Dynamic_vec_4(:,1)).^2 + (Dynamic_vec_4(:,2)).^2 + (Dynamic_vec_4(:,3)).^2);
Mag_Dynamic_vec_4 = mean(Mag_Dynamic_vec_4);


%% Unit vectors for the static and dynamic vector in plane of pendulum swing
%Static Vectors
unit_Static_4 = NewStatic_vector_4./(Mag_Static_4);
%Dynamic Vectors
%Small Swing
unit_Dynamic_vec_4 = NewDynamic_vec_4./(Mag_Dynamic_vec_4);
%% Find angle between static reference vector and changing dynamic vector
% Angle between static and slow dynamic
%Sine Rule
SinTheta_1 = cross(unit_Static_4,unit_Dynamic_vec_4);
Theta_1 = 180.*(asin(SinTheta_1))./pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of (Vicon data) %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% MotionSense Data - angle
MotionsenseAngle = MS_Sensor_Dynamic_new(:,6);
% correct MotionSense y offset
MotionsenseAngle_y_offset = max(MotionsenseAngle) + min(MotionsenseAngle); % MotionsenseAngle_y_offset =  1.7139 so must be taken into account
MotionsenseAngle = MotionsenseAngle - MotionsenseAngle_y_offset/2;


%% Plot vicon angle against time

y = Theta_1(:,1);
Theta_1_y_offset = max(y) + min(y); % Nb : Theta_1_y_offset = 0.0068 so almost neglectable
y = y - Theta_1_y_offset/2; % correction 

figure(1);
plot(time_v,y,'b',time_v,MotionsenseAngle,'r') %blue curve = Vicon , red = MotionSense)
title('Angle comparison between Vicon and MotionSense')
xlabel('Time in sec')
ylabel('Angle in degrees')
legend('Vicon (reference) ','MotionSense (alternative)');
%% Synchronisation
% compute time difference based on positive slope zero crossing time
% indexes
threshold = 0; % zero crossing
[t0_posVicon,s0_posVicon,t0_negVicon,s0_negVicon]= crossing_V7(y,time_v,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
[t0_posMS,s0_posMS,t0_negMS,s0_negMS]= crossing_V7(MotionsenseAngle,time_v,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 

delta_t = t0_posMS(1) - t0_posVicon(1);
% Motionsense : new time vector (shifted) 
new_samples = fix((max(time_v) - delta_t)/dt_v);
new_timeMS = delta_t + (0:new_samples-1)'*dt_v;
% interp the MotionsenseAngle data to the new defined time vector (this has better time
% accuracy than simply shifting time to a round number of samples)
MotionsenseAngle_new = interp1(time_v,MotionsenseAngle,new_timeMS,'linear');
new_timeMS = new_timeMS - new_timeMS(1); % 
ViconAngle_new = interp1(time_v,y,new_timeMS,'linear');

figure(2);
plot(new_timeMS,ViconAngle_new,'b',new_timeMS,MotionsenseAngle_new,'r') %blue curve = Vicon , red = MotionSense)
title('Angle comparison between Vicon and MotionSense')
xlabel('Time in sec')
ylabel('Angle in degrees')
legend('Vicon (reference) ','MotionSense (alternative)');

% now let's take only the first 322 samples ( = min(mViconS,mViconD)) to not have the trailing zeros 
% of Vicon's data in the error computation
ind = 1:min(mViconS,mViconD);
new_timeMS = new_timeMS(ind);
ViconAngle_new = ViconAngle_new(ind);
MotionsenseAngle_new = MotionsenseAngle_new(ind);

figure(3);
plot(new_timeMS,ViconAngle_new,'b',new_timeMS,MotionsenseAngle_new,'r') %blue curve = Vicon , red = MotionSense)
title('Angle comparison between Vicon and MotionSense')
xlabel('Time in sec')
ylabel('Angle in degrees')
legend('Vicon (reference) ','MotionSense (alternative)');


%% Calculate Error between the two systems
%root mean square error (RMSE)and Bland-Altman plots
Error = ViconAngle_new - MotionsenseAngle_new; %Error between Vicon and sensor
SqError = Error.^2;                          %Square error
AvSqError = mean(SqError);                   %Mean Squared Error  
RMSE = sqrt(AvSqError);                      %Root Mean Squared Error
%Pearson’s R correlation 
R = corrcoef(ViconAngle_new,MotionsenseAngle_new);
%Bland_Altman plot

[Mathieu NOE]

and this is another variant, without the need for zero padding the data

of course, this means at the end the common data is shorter .

% Comparing Vicon (reference sensor / data) to MotionSense (other sensor / data)
% Only 1 dynamic trial carried out
clearvars
% close all

%% Open Desired Files
% %Open Static 3 Marker files 
% Vicon3Static = readmatrix('C:\Users\lexil\Documents\PhD\Code\ViconData\20210614\Static_3M_Sensor.csv');
Vicon3Static = readmatrix('Static_3M_Sensor.csv');  % (reference sensor / data)
[mViconS,nViconS] = size(Vicon3Static);
% %Open Dynamic 3 Marker files 
% Vicon3Dynamic = readmatrix('C:\Users\lexil\Documents\PhD\Code\ViconData\20210614\Dynamic_3M_Sensor01.csv');
Vicon3Dynamic = readmatrix('Dynamic_3M_Sensor01.csv');  % (reference sensor / data)
[mViconD,nViconD] = size(Vicon3Dynamic);

% keep only the min common section (zero padding to make the shorter data longer is the other option)

if mViconS <= mViconD
    Vicon3Dynamic = Vicon3Dynamic(1:mViconS,:); % keep only mViconS samples of longer Vicon3Dynamic
else
    Vicon3Static = Vicon3Static(1:mViconD,:); % keep only mViconD samples of longer Vicon3Static (not the case here )
end

    % Vicon time vector 
    fs_vicon = 100;
    dt_v = 1/fs_vicon;
    time_v =  Vicon3Dynamic(:,1)*dt_v;
    time_v =  time_v - time_v(1); % so it starts at zero

% %Open motion sensor file
% Sensor_Static = readmatrix('C:\Users\lexil\Documents\PhD\Code\MotionSenseData\14_06_2021\67CT_log_14062021_12-19-01.csv');
% Sensor_Dynamic = readmatrix('C:\Users\lexil\Documents\PhD\Code\MotionSenseData\14_06_2021\67CT_log_14062021_12-26-20.csv');
MS_Sensor_Static = readmatrix('67CT_log_14062021_12-19-01.csv');   % (other sensor / data)
MS_Sensor_Dynamic = readmatrix('67CT_log_14062021_12-26-20.csv');   % (other sensor / data)

    % now resample the Dynamic data according to "reference" time vector = time_v
    % motion sensor time vector  (approx 50 Hz)
    dt_s = 1/2^15;
    time_s =  MS_Sensor_Dynamic(:,1)*dt_s;
    time_s =  time_s - time_s(1); % so it starts at zero
    MS_Sensor_Dynamic_new = interp1(time_s,MS_Sensor_Dynamic,time_v,'linear');
    
%     new_time_s =  NewDynamic_vec_4(:,1)*dt_v; % (check / debug only)
%     new_time_s =  new_time_s - new_time_s(1); % ok this proves interpolation works fine (check / debug only)

%% Marker Points (Vicon data)
% % Marker Positions for static marker 
marker_1_S = Vicon3Static(:,3:5)/1000;
marker_2_S = Vicon3Static(:,6:8)/1000;
marker_3_S = Vicon3Static(:,9:11)/1000;
% Marker Positions for Dynamic marker 
marker_1_D = Vicon3Dynamic(:,3:5)/1000;
marker_2_D = Vicon3Dynamic(:,6:8)/1000;
marker_3_D = Vicon3Dynamic(:,9:11)/1000;
%% Create Vectors (Vicon data)
% %Static Vectors
Static_vector_1 = marker_2_S - marker_1_S;
Static_vector_2 = marker_3_S - marker_1_S;
Static_vector_3 = cross(Static_vector_1, Static_vector_2);
Static_vector_4 = cross(Static_vector_3, Static_vector_1); %Vector in pendulums swining plane
%Dynamic Vectors
Dynamic_vec_1 = marker_2_D - marker_1_D;
Dynamic_vec_2 = marker_3_D - marker_1_D;
Dynamic_vec_3 = cross(Dynamic_vec_1, Dynamic_vec_2);
Dynamic_vec_4 = cross(Dynamic_vec_3, Dynamic_vec_1); %Vector in pendulums swining plane


%% Vector Magnitudes for pendulum plane
%Static Vectors
Mag_Static_4 = sqrt((Static_vector_4(:,1)).^2 + (Static_vector_4(:,2)).^2 + (Static_vector_4(:,3)).^2);
Mag_Static_4 = mean(Mag_Static_4);

%Dynamic Vectors
Mag_Dynamic_vec_4 = sqrt((Dynamic_vec_4(:,1)).^2 + (Dynamic_vec_4(:,2)).^2 + (Dynamic_vec_4(:,3)).^2);
Mag_Dynamic_vec_4 = mean(Mag_Dynamic_vec_4);


%% Unit vectors for the static and dynamic vector in plane of pendulum swing
%Static Vectors
unit_Static_4 = Static_vector_4./(Mag_Static_4);
%Dynamic Vectors
%Small Swing
unit_Dynamic_vec_4 = Dynamic_vec_4./(Mag_Dynamic_vec_4);
%% Find angle between static reference vector and changing dynamic vector
% Angle between static and slow dynamic
%Sine Rule
SinTheta_1 = cross(unit_Static_4,unit_Dynamic_vec_4);
Theta_1 = 180.*(asin(SinTheta_1))./pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of (Vicon data) %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% MotionSense Data - angle
MotionsenseAngle = MS_Sensor_Dynamic_new(:,6);
% correct slight y offset (MN)
MotionsenseAngle_y_offset = max(MotionsenseAngle) + min(MotionsenseAngle); % MotionsenseAngle_y_offset =  1.7139 so must be taken into account
MotionsenseAngle = MotionsenseAngle - MotionsenseAngle_y_offset/2;


%% Plot vicon angle against time

y = Theta_1(:,1);
Theta_1_y_offset = max(y) + min(y); % Nb : Theta_1_y_offset = 0.0068 so almost neglectable
y = y - Theta_1_y_offset/2; % correction 

figure(1);
plot(time_v,y,'b',time_v,MotionsenseAngle,'r') %blue curve = Vicon , red = MotionSense)
title('Angle comparison between Vicon and MotionSense')
xlabel('Time in sec')
ylabel('Angle in degrees')
legend('Vicon (reference) ','MotionSense (alternative)');
%% Synchronisation
% compute time difference based on positive slope zero crossing time
% indexes
threshold = 0; % zero crossing
[t0_posVicon,s0_posVicon,t0_negVicon,s0_negVicon]= crossing_V7(y,time_v,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
[t0_posMS,s0_posMS,t0_negMS,s0_negMS]= crossing_V7(MotionsenseAngle,time_v,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 

delta_t = t0_posMS(1) - t0_posVicon(1);
% Motionsense : new time vector (shifted) 
new_samples = fix((max(time_v) - delta_t)/dt_v);
new_timeMS = delta_t + (0:new_samples-1)'*dt_v;
% interp the MotionsenseAngle data to the new defined time vector (this has better time
% accuracy than simply shifting time to a round number of samples)
MotionsenseAngle_new = interp1(time_v,MotionsenseAngle,new_timeMS,'linear');
new_timeMS = new_timeMS - new_timeMS(1);
ViconAngle_new = interp1(time_v,y,new_timeMS,'linear');

figure(2);
plot(new_timeMS,ViconAngle_new,'b',new_timeMS,MotionsenseAngle_new,'r') %blue curve = Vicon , red = MotionSense)
title('Angle comparison between Vicon and MotionSense')
xlabel('Time in sec')
ylabel('Angle in degrees')
legend('Vicon (reference) ','MotionSense (alternative)');


%% Calculate Error between the two systems
%root mean square error (RMSE)and Bland-Altman plots
Error = ViconAngle_new - MotionsenseAngle_new; %Error between Vicon and sensor
SqError = Error.^2;                          %Square error
AvSqError = mean(SqError);                   %Mean Squared Error  
RMSE = sqrt(AvSqError);                      %Root Mean Squared Error
%Pearson’s R correlation 
R = corrcoef(ViconAngle_new,MotionsenseAngle_new);
%Bland_Altman plot

[alexandra ligeti]

Thank you so much for all your work.

I was still stuck working on it today. Will spend some time going through your solution to make sure I understand it all. Sorry, for providing all this confusion with the horribly labeled data sets. Really appreciate the help.

Thank you once again!

[Mathieu NOE]

My pleasure !!