Sampling and resampling data of different lengths

13 vues (au cours des 30 derniers jours)
alexandra ligeti
alexandra ligeti le 18 Juin 2021
Commenté : Mathieu NOE le 23 Juin 2021
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
  6 commentaires
Mathieu NOE
Mathieu NOE le 23 Juin 2021
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
alexandra ligeti le 23 Juin 2021
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.

Connectez-vous pour commenter.

Réponse acceptée

Mathieu NOE
Mathieu NOE le 23 Juin 2021
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
  3 commentaires
alexandra ligeti
alexandra ligeti le 23 Juin 2021
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
Mathieu NOE le 23 Juin 2021
My pleasure !!

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by