EMG data comparison problem

22 views (last 30 days)
Jose Colon
Jose Colon on 2 Feb 2020
Answered: Francisco Romero on 11 Feb 2022
I have some sEMG data sets that I want to compare. How can I do this with Matlab?
%% Analysis of EMG data
%An EMG signal has a frequency between 0 to 500 Hz; but the usable
%frequency is between 50 and 150 Hz.
% EMG signal is produced using an instrument called the electromyograph that records electric activity
% produced by muscles.
% The amplitude for surface EMG signals is within the range of 0 to 10 mV.
% sEMG (surace EMG) signals tend to have a DC bias, therefore a high pass
% filter is applied to a cutoff frequency of 50Hz. This high pass filter
% removes any DC bias.
% The root mean-square (RMS): calculated by finding the square root of: squaring each data point,
%summing the squares, dividing the sum by the number of observations.
% RMS value of a signal is a measure of the power of the signal and it
% reflects on the activity in the motor unit during contraction.
% This makes RMS useful in producing waveforms.
% Amplitude of an EMG signal is random with Gaussian Distribution/
% RMS and MA are two parameters used to measure the amplitude.
% MA short for: Mean Absolute value: is a mean of the absolute value of
% EMG signal
% MA is a measure of the area under the rectified EMG.
% In order to calculate MA remove all of the negative phases of the raw EMG
% This is known as: half-wave rectification.
% The MA and RMS provide an estimate of the amplitude of the raw EMG signal.
clear,clc
load('A1_Ball_high_t2.mat')
t = t2(:,1);
y1 = t2(:,2);
N = length(y1);% find the length of the data per second
ls = size(y1); %% size
f = 1/N;% find the sampling rate or frequency
%Note: sampling frequency should be at least twice the incoming signal.
fs = 2000;
T = 1/fs % period between each sample
t1 = (0 : N-1) *T;
t = (0:1:length(y1)-1)/fs; % sampling period
% nyquist frequency is half of the sampling rate of a signal. Also known as:
% folding frequency.
Nyquist = fs/2;
figure;
%
subplot (3,1,1), plot(t,y1,'b');
title ('EMG signal of single muscle 25 years old patient ');
xlabel ('time (sec)');
ylabel ('Amplitute (V)');
grid on;
Y= abs(fft(y1));
Y(1) = [];
power = abs(Y(1:N/2)).^2;
nyquist = 350;
freq = (1:N/2)/(N/2)*nyquist;
subplot(212), plot(freq,power), grid on
xlabel('Sample number (in Frequency)')
ylabel('Power spectrumen');
title({'Single-sided Power spectrum' ...
' (Frequency in shown on a log scale)'});
axis tight
%%% RMS of the signal
rms_y1 = sqrt(mean(y1.^2));
msgbox(strcat('RMS of EMG signal is = ',mat2str(rms_y1), ''));
rms_emg = rms (y1);
%%%%% MA of the signal
ma_y1 = abs(mean(y1));
msgbox(strcat('MA of EMG signal is = ',mat2str(ma_y1), ''));
%remove any DC offset of the signal
% use the detrend function.
% It subtracts the mean or a best-fit line from data.
y2 = detrend(y1); %% y2 is the singal after DC offset removed.
figure;
rec_y = abs(y2);
plot (rec_y);
xlabel('Sample number (in Frequency)')
ylabel('Rectified EMG signal');
title({'Rectified EMG signal' ...
' (Frequency shown on a log scale)'});
figure;
xdft = fft(y1);
xdft = xdft(1:N/2+1);
psdx = (1/(fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:fs/length(y1):fs/2;
plot(freq,10*log10(psdx))
grid on
title(' Power spectrum FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
%amplitude of the entire signal
amplitude = max(y2(:)) - min(y2(:))
%%% DFT to describe the signal in the frequency
NFFT = 2 ^ nextpow2(N);
% The reason I used this function above:
%The execution time for fft depends on the length of the transform.
% It is fastest for powers of two. It is almost as fast for lengths that
%have only small prime factors. It is typically several times slower for
%lengths that are prime or which have large prime factors.
Y = fft(y1, NFFT) / N;
f = (fs / 2 * linspace(0, 1, NFFT / 2+1))'; % Vector containing frequencies in Hz
amp = ( 2 * abs(Y(1: NFFT / 2+1))) % Vector containing corresponding amplitudes
figure;
plot (f, amp);
title ('plot single-sided amplitude spectrum of the EMG signal')
xlabel ('frequency (Hz)')
ylabel ('|y(f)|')
grid on;
References:
Ali H. Al-Timemy, Rami N. Khushaba, Guido Bugmann, and Javier Escudero, " Improving the Performance against Force Variation of EMG Controlled Multifunctional Upper-Limb Prostheses for Transradial Amputees", IEEE Transactions on Neural Systems and Rehabilitation Engineering, Accepted for Publication, 2015.
http://biomedicalsignalandimage.blogspot.com/2016/02/matlab-code-to-study-emg-signal.html

Answers (1)

Francisco Romero
Francisco Romero on 11 Feb 2022
Hi José,
I saw to this question by looking for other things.
You may find this link useful:

Tags

Community Treasure Hunt

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

Start Hunting!