How to split my EMG signal into individual cycles?

7 views (last 30 days)
Hi I have recorded an emg signal of 5 different muscles on someone performing manual wheelchair propulsion.
Sensor 1 = biceps
Sensor 2 = triceps
Sensor 3 = Anterioir Deltoid
Sensor 4 = Posterior Deltoid
Sensor 5 = Upper Trapezius
Each sensor also records the acceleration in the x, y and z directions.
I want to split the emg signal into seperate propulsion cycles (when the hand pushes and releases to when the hand hits the push rim again.) so that I can average each cycle, as well as convert time (s) to percentage of cycle.
I know there were 13-14 cycles during the recording as i videoed the subject while gathering data.
I have only looked at the bicep data as shown below and was wondeing if there is a way to do this with the information I have.
Please let me know if you need anymore information.
%% SEMG Signal Analysis - Row Test 1 Participant 1 28/06/2024
% *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. *
clear all;
clc;
%% Importing raw data from csv. file
% *BOLD TEXT*
raw_data = readmatrix('P02_S09_T01.csv');
Error using readmatrix (line 171)
Unable to find or open 'P02_S09_T01.csv'. Check the path and filename or file permissions.
%% SEMG Time Data
% *Storing the Sub-frame data in arrays*
fs = 1926; % sampling frequency in Hz
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;
while ~isnan(raw_data(row, columnIndex))
number = raw_data(row, columnIndex);
Sub_Frame = [Sub_Frame, number];
row = row + 1;
end
Time = (1:length(Sub_Frame))/fs; % converting the sub frame array into time in seconds array
%% 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
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_row = 77955;
while acc_row <= maxRow
number = raw_data(acc_row, columnIndex);
Acc_Sub_Frame = [Acc_Sub_Frame, number];
acc_row = acc_row + 1;
end
Acc_Time = (1:length(Acc_Sub_Frame))/acc_fs; % converting the sub frame array into time in seconds array
columnIndex4 = 4; % Corresponds to ACCY1 (acceleration in the y-direction of sensor 1).
ACCY1 = []; % Creating an array to store all the AccY1
row2 = 77955;
while row2 <= maxRow
num = raw_data(row2, columnIndex4);
ACCY1 = [ACCY1, num];
row2 = row2 + 1;
end
acc_rms_signal = [];
window = 50;
acc_rms_signal = sqrt(movmean((ACCY1.^2),window));
%% Muscle arrays
% Sensor #1 - Biceps Brachii
% Sensor #2 - Triceps Brachii
% Sensor #3 - Anterioir Deltoid
% Sensor #4 - Posterior Deltoid
% Sensor #5 - Upper Trapezius
% *Storing the muscle data in arrays*
columnIndex = 3; % Specify the index of the column you want to extract (1 for the Frame column, 2 for the Sub Frame column, and so on)
Biceps_Brachii = [];% Creating an array to store all the frames related to the EMG data
row = 4; % This is the starting position of the Frame data
% Create a for loop to run through every sensor and store in each array.
while ~isnan(raw_data(row, columnIndex))
number = raw_data(row, columnIndex);
Biceps_Brachii = [Biceps_Brachii, number];
row = row + 1;
end
%% BP filter
% filtering from 20 to Hz
fnyq = fs/2;
fcuthigh = 20;
fcutlow = 450;
%4th order Butterworth BP filter
[b,a]=butter(4,[fcuthigh,fcutlow]/fnyq,'bandpass');
butterworth_filter_signal=filtfilt(b,a,Biceps_Brachii);
%% RMS
%
rms_signal = [];
window = 50;
rms_signal = sqrt(movmean((rec_signal.^2),window));
% Create a figure with two subplots
figure;
% Plot the first graph on the top
subplot(2, 1, 1);
plot(Time, rms_signal);
title('EMG');
% Plot the second graph on the bottom
subplot(2, 1, 2);
plot(Acc_Time, acc_rms_signal);
title('Acceleration');

Accepted Answer

Umar
Umar on 12 Aug 2024
Edited: Umar on 12 Aug 2024

Hi @Renee Wurfel ,

When analyzing your code, I did notice the error you're encountering when trying to read the CSV file in MATLAB, let's break down the key areas to check and troubleshoot. Make sure that the file `P02_S09_T01.csv` is located in the current working directory of MATLAB. You can check your current directory by using the command:

    pwd

If the file is not in this directory, you can either move it there or specify an absolute path in your `readmatrix` function. For example:

raw_data = readmatrix('C:\path\to\your\file\P02_S09_T01.csv');

Now, addressing your query regarding, “To split the EMG signal into propulsion cycles, you need to identify the start and end points of each cycle. Since you have video documentation, you can visually confirm the timing of hand pushes and releases. Typically, this can be done by detecting peaks in the acceleration data or specific patterns in the EMG signal that correspond to the push phase.”

Please see my step by step instructions below.

Identify Propulsion Cycles

I do agree with @Taylor comments, “I would recommend using the findpeaks function. You can set a minimum amplitude value (Threshold) and minimum time between cycles (MinPeakDistance) if you have a rough idea of the frequency of contration cycles.” Use MATLAB’s `findpeaks` function on the acceleration or filtered EMG data to identify when the hand is pushing against the rim. For example:

   [pks, locs] = findpeaks(acc_rms_signal, 'MinPeakHeight', threshold_value);

Here, `threshold_value` should be set based on your data characteristics to ensure only significant peaks are detected.

Define Cycles

Once peaks are identified, define each cycle by pairing successive peaks (push release) and calculate their indices for segmentation.

Segment EMG Data

With the identified cycle indices, segment your EMG data accordingly:

   cycles = []; % Initialize a matrix to store segmented cycles
   for i = 1:length(locs)-1
    start_idx = locs(i);
    end_idx = locs(i+1);
    % Store each cycle in a cell array
     cycles{i} = Biceps_Brachii(start_idx:end_idx);
     end

Average Each Cycle

After segmentation, compute the average for each cycle:

average_cycles = cellfun(@mean, cycles);

This will give you a vector of average EMG values for each propulsion cycle.

Convert Time to Percentage of Cycle

To convert time into percentage of cycle duration:

cycle_durations = cellfun(@length, cycles);
percent_cycles = cellfun(@(x) (1:length(x)) / length(x) * 100, cycles, 
‘UniformOutput', false);

This will yield a cell array where each entry corresponds to the percentage time for that specific cycle.Also, visualizing each segmented cycle alongside their averages can provide insights into muscle activation patterns throughout propulsion. I will use plotting functions to create clear visual representations:

     figure;
     hold on;
     for i = 1:length(cycles)
       plot(percent_cycles{i}, cycles{i});
     end
     plot(average_cycles, 'k--', 'LineWidth', 2); % Plot average as dashed line
     title('EMG Signal During Propulsion Cycles');
     xlabel('Percentage of Cycle (%)');
     ylabel('EMG Signal (mV)');
     hold off;

Repeat similar steps for other muscle sensors (triceps, deltoids, trapezius) using their respective data arrays after confirming consistency in cycle identification across all muscles. Also, depending on your research goals, consider analyzing frequency domain features using techniques like Fast Fourier Transform (FFT) or investigating muscle co-activation patterns during propulsion. If you still require further assistance or specific code implementations for additional analyses, feel free to ask!

  5 Comments
Renee Wurfel
Renee Wurfel on 19 Aug 2024
Thank you for your help Umar it is much appreciated i ended up adding a minpeakdistance and it worked!
threshold_value = 8500;
[pks, locs] = findpeaks(acc_rms_signal, Acc_Time, 'MinPeakHeight', threshold_value, 'MinPeakDist', 0.5);
plot(Acc_Time,acc_rms_signal,locs,pks,"o")
xlabel("Time")
ylabel("Acceleration")
axis tight
legend(["Data" "peaks"],Location="NorthWest")
Umar
Umar on 19 Aug 2024
Hi @ Renee Wurfel,
Thank you for your kind words. I am pleased to hear that adding the minpeakdistance resolved your issue. It’s always rewarding to know that my assistance has been helpful. If you have any further questions or need additional support, please don’t hesitate to reach out.

Sign in to comment.

More Answers (1)

Taylor
Taylor on 12 Aug 2024
I would recommend using the findpeaks function. You can set a minimum amplitude value (Threshold) and minimum time between cycles (MinPeakDistance) if you have a rough idea of the frequency of contration cycles.

Products


Release

R2023b

Community Treasure Hunt

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

Start Hunting!