-
Notifications
You must be signed in to change notification settings - Fork 0
/
EMG_processing_tutorial.m
executable file
·161 lines (139 loc) · 5.11 KB
/
EMG_processing_tutorial.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
% Filename: EMG_processing_tutorial.m
% Author: Samuel Acuña
% Date: 11 May 2018
% Description:
% Guides you through a typical EMG processing routine.
%
% EMG is more about timing of activation, than the amplitude.
%
% Make sure you have the accompanying data file:
% EMG_processing_tutorial_DATA.mat or
% EMG_processing_tutorial_DATA2.mat
%
% Directions: run each section at a time to see steps (usually alt+enter, or command+enter)
%% STEP 1: LOAD EMG DATA
clear; close all; clc;
%load('EMG_processing_tutorial_DATA.mat');
load('EMG_processing_tutorial_DATA2.mat');
% DATA: this is data from the medial gastrocnemius when walking on a treadmill
% DATA2: tibialis anterior, same trial
% Loads matlab structure: trial_data
% trial_data.emg : raw EMG data
% trial_data.label : muscle where EMG data was collected
% trial_data.freq : sampling frequency
% trial_data.time : time of each sample of data
% trial_data.heelStrikes : instances of heel strikes between steps
return
%% STEP 1: EXAMINE RAW EMG DATA
figure(1);
h1 = plot(trial_data.time, trial_data.emg,'color',[ 0 0.4470 0.7410]);
title(['EMG: ' trial_data.label]);
xlabel('time (sec)'); xlim([0 60]);
ylabel('EMG (volts)');
legend([h1],'raw EMG');
%% STEP 2: VISUALIZE HEEL STRIKES
figure(1); hold on;
scale = max(trial_data.emg);
hs_x = [trial_data.heelStrikes, trial_data.heelStrikes]';
hs_y = [zeros(length(trial_data.heelStrikes),1)-(scale/2), zeros(length(trial_data.heelStrikes),1)+(scale/2)]';
h2 = plot(hs_x,hs_y,'r');
hold off;
legend([h1 h2(1)],{'raw EMG', 'heel strikes'});
%% STEP 3: BANDPASS FILTER EMG
% the bandpass gets rid of any noise that isnt part of EMG data. Most EMG
% power is between 5-500 Hz
% low cutoff: want to be high enough to get rid of motion artifacts and drift. No higher than 10 hz, per ISEK guidelines
% high cutoff: no lower than 350 Hz, per ISEK guidelines
BP = [10 500]; % bandpass filter parameters in Hz [low cutoff, high cutoff]
[b_BP,a_BP]=butter(4,BP/(trial_data.freq/2)); % bandpass filter coefficients. (4th order butterworth)
emg_BP = filtfilt(b_BP,a_BP,trial_data.emg); % zero phase lag filter
figure(1); hold on;
h3 = plot(trial_data.time,emg_BP,'color',[0.6875 0.7656 0.8672]);
hold off;
legend([h1 h2(1) h3],{'raw EMG', 'heel strikes', 'bandpass EMG'});
%% STEP 4: FULL-WAVE RECTIFY EMG
% just the absolute value
emg_ABS = abs(emg_BP);
figure(1); hold on;
h4 = plot(trial_data.time,emg_ABS,'color',[0.6758 0.8438 0.8984]);
hold off;
legend([h1 h2(1) h3 h4],{'raw EMG', 'heel strikes', 'bandpass EMG', 'rectify EMG'});
%% STEP 5: LINEAR ENVELOPE
% low pass filter
LP = 10; % low pass filter for linear envelope, in Hz
[b_LP,a_LP]=butter(4,LP/(trial_data.freq/2),'low'); % linear envelope filter
emg_ENV = filtfilt(b_LP,a_LP,emg_ABS);
figure(1); hold on;
h5 = plot(trial_data.time,emg_ENV,'color',[0 0 0],'LineWidth',2);
hold off;
legend([h1 h2(1) h3 h4 h5],{'raw EMG', 'heel strikes', 'bandpass EMG', 'rectify EMG', 'envelope EMG'});
%% STEP 6: NORMALIZE AMPLITUDE
% scaled to the max value (could also do RMS, or max voluntary contraction)
emg_NORM = emg_ENV/max(emg_ENV);
figure(2);
h6 = plot(trial_data.time,emg_NORM,'color',[ 0 0.4470 0.7410],'LineWidth',2);
title(['EMG: ' trial_data.label]);
xlabel('time (sec)'); xlim([0 60]);
ylabel('normalized EMG'); ylim([0 1.5]);
legend([h6],'normalized linear envelope of EMG');
%% STEP 7: TIME NORMALIZE BY A STRIDE
% divides the EMG signals insto strides
npts = 101; % points per gait cycle
nStrides = length(trial_data.heelStrikes)-1; % number of complete gait cycles
emg_Strides = zeros(101,nStrides); %preallocate
for j = 1:nStrides
j1 = find(trial_data.time>trial_data.heelStrikes(j),1); % get index of time at first heel strike
j2 = find(trial_data.time>trial_data.heelStrikes(j+1),1); % get index of time at second heel strike
emg_Strides(:,j) = normcycle(emg_NORM(j1:j2),npts); % time normalize
end
% FIND AVERAGE EMG over a stride
emg_AVG = mean(emg_Strides,2); %average EMG of each stride
emg_STD = std(emg_Strides')'; % standard deviation
% PLOT
figure(3); subplot(2,1,1);
shadedErrorBar([0:100]',emg_AVG,emg_STD);
title(['EMG: ' trial_data.label]);
xlabel('Gait Cycle (0-100%)'); xlim([0 100]);
ylabel('normalized EMG'); ylim([0 1.5]);
legend('EMG (AVG ± STD)')
figure(3); subplot(2,1,2); hold on;
for i = 1:nStrides
plot([0:100]',emg_Strides(:,i));
end
hold off;
title(['EMG for every stride']);
xlabel('Gait Cycle (0-100%)'); xlim([0 100]);
ylabel('normalized EMG'); ylim([0 1.5]);
function yf = normcycle(y,n,x)
% yf = normcycle(y,n,x)
% Convert a signal y to n even-spaced data points over a cycle
% Often used for presentation of gait data, default for n is 101 points
% can specify an indpendent variable x (optional)
if ~exist('n','var')
n=101;
end
[nr,nc]=size(y);
if nc==1 && nr>1
ny=1;
nx=nr;
elseif nr==1 && nc>1
y=y';
ny=1;
nx=nc;
elseif nr>1 && nc>1
ny=nc;
nx=nr;
else
disp('normcycle does not work on a scalar value');
yf=[];
return
end
if ~exist('x','var')
x=[0:(nx-1)]/(nx-1);
else
nx=length(x);
x=(x-x(1))/(x(end)-x(1));
end
kk=[0:(n-1)]/(n-1);
yf=interp1(x,y,kk,'*pchip');
end