How to normalise segregated surface emg signals into the same number of data points?

Hi everyone,

I have a cell array containing multiple gait cycles of an SEMG signal. Each cycle represents the push phase and recovery phase of a wheelchair. I want to get the average of the SEMG signal of all the cycles however, each cycle contains a different amount of data points. How do I normalise the cycles to contain a certain amount of data points without disrupting the signals? I also have a cell array containg the times of propulsion.

If anyone could point me in the right direction that would be great!

Thanks,

[Jacob Mathew]

Hey Renee,

Could you provide some sample or subset of the SEMG Signal data and the times of propulsion cell array?

[Renee Wurfel]

%% SEMG Signal Analysis - Wheelchair Analysis Practice 1
% *This is my practice code to figure out how to export the data into MATLAB and transfer the data into arrays that I can manipulate and do calculations with.  * 
% https://au.mathworks.com/matlabcentral/playground/new  AI chat to help
%-------------------------------------------------------------------------%
clear all;
clc;
%% Importing raw data from csv. file

raw_data = readmatrix('P02_S09_T01.csv');

%% SEMG Time Data
% *Storing the Sub-frame data in arrays* 

fs = 1926; % sampling frequency in Hz

SF_columnIndex = 2; % Specify the index of the column you want to extract (1 for the Frame column, 2 for the Sub Frame column, and so on)
Sub_Frame = []; % Creating an array to store all the sub-frames related to the EMG data

row = 4; % starting position
while ~isnan(raw_data(row, SF_columnIndex))
    number = raw_data(row, SF_columnIndex);
    Sub_Frame = [Sub_Frame, number];
    row = row + 1;
end
% Converting sub-frame array into time/s by dividing by fs (1926 Hz)
Time = (1:length(Sub_Frame))/fs; % converting the sub frame array into time in seconds array

%% Muscle arrays
% Creating muscle arrays to store all the emg data related to that muscle
% in the specified array
% Sensor #1 = Biceps Brachii
% Sensor #2 = Triceps Brachii
% Sensor #3 = Anterior Deltoid
% Sensor #4 = Posterior Deltoid
% Sensor #5 = Upper Trapezius 

%-------------------------------------------------------------------------%
% Sensor #1 - Biceps Brachii
Biceps_Brachii = [];

BB_columnIndex = 3; % Specify the index of the column you want to extract
row = 4; % This is the starting position of the Frame data

% while there is a number in the row store in array
while ~isnan(raw_data(row, BB_columnIndex))
    number = raw_data(row, BB_columnIndex);
    Biceps_Brachii = [Biceps_Brachii, number];
    row = row + 1;
end
%-------------------------------------------------------------------------%
% Sensor #2 - Triceps Brachii
Triceps_Brachii = [];

TB_columnIndex = 4;
row = 4;
while ~isnan(raw_data(row, TB_columnIndex))
    number = raw_data(row, TB_columnIndex);
    Triceps_Brachii = [Triceps_Brachii, number];
    row = row + 1;
end
%-------------------------------------------------------------------------%
% Sensor #3 - Anterioir Deltoid
Anterior_Deltoid = [];

AD_columnIndex = 5;
row = 4;

while ~isnan(raw_data(row, AD_columnIndex))
    number = raw_data(row, AD_columnIndex);
    Anterior_Deltoid = [Anterior_Deltoid, number];
    row = row + 1;
end
%-------------------------------------------------------------------------%
% Sensor #4 - Posterior Deltoid
Posterior_Deltoid = [];

PD_columnIndex = 6;
row = 4;

while ~isnan(raw_data(row, PD_columnIndex))
    number = raw_data(row, PD_columnIndex);
    Posterior_Deltoid = [Posterior_Deltoid, number];
    row = row + 1;
end
%-------------------------------------------------------------------------%
% Sensor #5 - Upper Trapezius
Upper_Trapezius = [];

UT_columnIndex = 7;
row = 4;

while ~isnan(raw_data(row, UT_columnIndex))
    number = raw_data(row, UT_columnIndex);
    Upper_Trapezius = [Upper_Trapezius, number];
    row = row + 1;
end
%-------------------------------------------------------------------------%
% *Plotting raw muscle signals on one plot* 

% plot(Time,Biceps_Brachii)
% hold on
% plot(Time,Triceps_Brachii)
% hold on
% plot(Time,Anterior_Deltoid)
% hold on
% plot(Time,Posterior_Deltoid)
% hold on
% plot(Time,Upper_Trapezius)
% hold on
% %Add labels to the x and y axes
% title('muscle activation');
% xlabel('Time (s)');
% ylabel('Voltage (V)'); % Confirmed it is in Volts
% legend('biceps','tricpes','AD','PD','UT');


%% Acceleration Sensor 1 (bicep) ACCY1

% Find the maximum row number
maxRow = size(raw_data, 1);
acc_fs = 148; %sa/sec

%find sub frame for acceleration
BB_columnIndex = 2; % Specify the index of the column you want to extract (1 for the Frame column, 2 for the Sub Frame column, and so on)
Acc_Sub_Frame = []; % Creating an array to store all the sub-frames related to the EMG data

acc_sub_frame = 77955; % starting row for the acceleration data
while acc_sub_frame <= maxRow
    number = raw_data(acc_sub_frame, BB_columnIndex);
    Acc_Sub_Frame = [Acc_Sub_Frame, number];
    acc_sub_frame = acc_sub_frame + 1;
end
% converting the acceleration sub frame array into time in seconds array
Acc_Time = (1:length(Acc_Sub_Frame))/acc_fs; 

columnIndex4 = 4; % Corresponds to ACCY1 (acceleration in the y-direction of sensor 1). 
ACCY1 = []; % Creating an array to store all the AccY1

acc_row = 77955; 

while acc_row <= maxRow
    num = raw_data(acc_row, columnIndex4);
    ACCY1 = [ACCY1, num];
    acc_row = acc_row + 1;
end

% plot(Acc_Time,ACCY1)
% title('Biceps Brachii Acceleration (Y-direction)');
% 
% % Add labels to the x and y axes
% xlabel('Time (s)');
% ylabel('Acceleration (mm/s²)');
% Applying rms to the raw acceleration signal (easier to find peaks)
acc_rms_signal = [];
window = 50;
acc_rms_signal = sqrt(movmean((ACCY1.^2),window));
% plot(Acc_Time,acc_rms_signal);


%% Identifying Propulsion Cycles
%  
threshold_value = 8500;
PeakDist = 0.5;
[pks, locs] = findpeaks(acc_rms_signal, Acc_Time, 'MinPeakHeight', threshold_value, 'MinPeakDist', PeakDist);

% disp(pks)
% disp(locs)

% plot(Acc_Time,acc_rms_signal,locs,pks,"o")
% xlabel("Time")
% ylabel("Acceleration")
% axis tight
% legend(["Data" "peaks"],Location="NorthWest")

% Finding Cycle locations in array Time

BB_Cycle_location=[];

for i = 1:length(locs)
    ValueToFind=locs(1,i);
    [~, index] = min(abs(Time-ValueToFind));
    closestValue = Time(index);
    BB_Cycle_location(1,i)=find(Time==closestValue);
end
disp(BB_Cycle_location)
disp(Time(BB_Cycle_location))

% Segregating EMG signal into cycles. 

% propulsion started 5s after recording 5*fs = starting pos
%Updated start time at 5.6 after looking at the signal and also found that
%the end time is around 21.7 sec
% It wont matter because I'm going to ignor the first and last cycle... 
Start_Index = round(5.6*fs);

BB_Cycles = {};
Time_Cycles = {};

for i = 1:length(locs)
   
    End_Index=BB_Cycle_location(1,i);
    BB_Cycles{i}=Biceps_Brachii(Start_Index:End_Index);
    Time_Cycles{i}=Time(Start_Index:End_Index);
    Start_Index=BB_Cycle_location(1,i);

end

[Renee Wurfel]

Hi Jacob,

I have attached the file im reading from as well as the code that I have written that segregates the time and semg signals into cycles.

Thanks!

[Umar]

Hi @Renee Wurfel,

Please update if you made any progress.

Answer 1

[Jacob Mathew]

Hey Renee,

I understand that you want to normalise the cycles to contain a certain number of datapoints without disrupting the signal. One way to achieve this is to use interpolation.

We will use the “interp1” function to achieve interpolation. The link to the function’s documentation is below:

https://www.mathworks.com/help/releases/R2024a/matlab/ref/interp1.html

A possible approach is to take the average length of the cycles and interpolate every cycle to have that length. This normalization can be achieved using the following code:

% Calculate the length of each cycle
cycle_lengths = cellfun(@length, BB_Cycles);

% Compute the average length of the cycles
% This will be the length we will normalise all cycles to
average_length = round(mean(cycle_lengths));

% Initialize cell arrays to store the normalized cycles
BB_Cycles_Normalized = cell(1,length(BB_Cycles));
Time_Cycles_Normalized = cell(1,length(BB_Cycles));

% Loop through each cycle to normalize
for i = 1:length(BB_Cycles)
    % Get the current cycle and its time
    current_cycle = BB_Cycles{i};
    current_time = Time_Cycles{i};
    
    % Create a new time vector for the normalized cycle
    normalized_time = linspace(current_time(1), current_time(end), average_length);
    
    % Interpolate the data to the new time vector
    normalized_cycle = interp1(current_time, current_cycle, normalized_time, 'linear');
    
    % Store the normalized cycle and its time
    BB_Cycles_Normalized{i} = normalized_cycle;
    Time_Cycles_Normalized{i} = normalized_time;
end

This will ensure that each cycle will have consistent number of points allowing comparison across cycles.

[Renee Wurfel]

Thanks Jacob much appreciated!

Answer 2

[Umar]

Hi @Renee Wurfel ,

To achieve the normalization of SEMG cycles while preserving the integrity of the signals, you can utilize interpolation techniques. This method will allow you to resample each cycle to a predetermined number of data points without distorting the signal characteristics. Here’s how you can modify your existing code to incorporate these steps:

% Assuming BB_Cycles contains your segmented SEMG cycles
num_cycles = length(BB_Cycles);
target_length = 1000; % Set your desired target length for normalization
normalized_cycles = zeros(num_cycles, target_length); % Preallocate normalized   cycles

% Interpolation for each cycle
for i = 1:num_cycles
  original_cycle = BB_Cycles{i};
  original_length = length(original_cycle);

    % Create a time vector for the original cycle
    original_time = linspace(0, 1, original_length); % Normalize time between 0   
    and 1

    % Create a new time vector for the target length
    new_time = linspace(0, 1, target_length);

    % Interpolate using linear method (you can use 'spline' for smoother results)
    normalized_cycles(i, :) = interp1(original_time, original_cycle, new_time, 
   'linear');
  end

% Calculate average SEMG signal across all normalized cycles
average_SEMG = mean(normalized_cycles, 1);

% Plotting the average SEMG signal
figure;
plot(linspace(0, 1, target_length), average_SEMG);
title('Average SEMG Signal Across Cycles');
xlabel('Normalized Time');
ylabel('Amplitude');
grid on;

% You may also want to visualize individual normalized cycles
figure;
hold on;
for i = 1:num_cycles
  plot(linspace(0, 1, target_length), normalized_cycles(i, :));
end
title('Normalized SEMG Cycles');
xlabel('Normalized Time');
ylabel('Amplitude');
legend(arrayfun(@(x) sprintf('Cycle %d', x), 1:num_cycles, 'UniformOutput',   
false));
hold off;

In addition to this code, the choice of interpolation method (linear vs. spline) can affect the smoothness of your normalized signals. Spline interpolation tends to provide smoother results but may introduce artifacts if not handled carefully. So, make sure that the normalization process does not introduce significant distortion in your SEMG signals. I will advise to visualize both individual cycles and their averages post-normalization. After obtaining the average signal, consider performing statistical analyses to evaluate variability among different cycles or conditions.

[Renee Wurfel]

Thank you very much Umar! I may need assitance with evaluating the variabilty among different cycles, but that is a great idea!

[Umar]

Hi @Renee Wurfel,

I appreciate your positive feedback regarding the idea. I would be more than happy to assist you with evaluating the variability among different cycles.