肌电信号的识别.doc

燕山大学
课 程 设 计 说 明 书

题目： 肌电信号分析及动作识别

学院（系）：电气工程学院
年级专业： 10级仪表三班
学 号：
学生姓名：
指导教师：
教师职称： 教授 讲 师

电气工程学院《课程设计》任务书
基层教学单位：自动化仪表系 指导教师：谢平 张淑清
学号

学生姓名

（专业）班级
10级仪表三班
设计题目
肌电信号分析及动作识别
设
计
技
术
参
数

1、独立完成设计任务
2、编程，上机调试。
3、连接硬件实验线路，实现所要求的功能。
4、完成设计，提交课程设计报告。

设
计
要
求
1、利用肌电传感器采集人体肌电信号，利用Matlab软件对肌电信号进行时域和频域(包括FFT，功率谱，倒谱)分析，计算均值、标准差、方差、积分肌电值IEMG、均方根RMS等时域指标参数和平均功率频率MPF、中值频率MF等频域指标。比较不同动作下肌电信号特征，根据设定阈值给出动作状态，设计相应的信号分析及显示界面。
2、基于肌电信号分析结果，通过串口发送命令给单片机系统，根据肌电信号动作状态控制相应的数码管或LED发光显示。
3、扩展：也可通过动作模式驱动电机转动或其他控制输出模块。
工
作
量

软件编程与硬件调试相结合，绘制设计流程图，编制相应软件界面，实现单片机控制与信号处理任务的综合应用
参
考
资
料
1）《微型计算机控制系统》 赖寿宏，机械工业出版社（教材）
2）《单片机及应用》李大友，高等教育出版社（教材）
3）《信号处理原理及应用》 谢平 等 机械工业出版社（教材）
4)《Matlab程序设计及其在信号处理中的应用》聂祥飞等 西南交通大学出版社
周次
第一周
第二周

应
完
成
内
容
熟悉伟福单片机编程环境，调试单片机各基本功能模块；熟悉matlab信号处理工具箱，信号处理系统基本功能模块学习和调试
单片机系统与信号处理系统综合进行硬件调试，
撰写课程设计报告

指导教
师签字

基层教学单位主任签字

目 录

第一章 摘要 ……………………………………………………………… 2

第二章 系统总体设计方案………………………………………………… 3

第三章 肌电信号的时域参数处理及其分析 …………………………… 4

第四章 肌电信号的频域处理方法及其分析 …………………………… 7

3.1 FFT分析 ……………………………………………………………7

3.2 功率谱分析…………………………………………………………8

3.3 倒谱分析……………………………………………………………9

3.4 平均功率频率MPF和中值频率………………………………… 10

第五章 Matlab程序及GUI…………………………………………………11

第六章 系统整体调试及结果说明 ………………………………………24

第七章 学习心得………………………………………………………… 24

参考文献 ……………………………………………………………………25

第一章 摘要
肌电信号是产生肌肉力的电信号根源，它是肌肉中很多运动单元动作电位在时间和空间上的叠加，反映了神经，肌肉的功能状态，在基础医学研究、临床诊断和康复工程中有广泛的应用。
其种类重要有两种：一，临床肌电图检查多采用针电极插入肌肉检测肌电图，其优点是干扰小，定位性好，易识别，但由于它是一种有创伤的检测方法，其应用收到了一定的限制。二，表面肌电则是从人体皮肤表面通过电极记录下来的神经肌肉活动时发放的生物电信号，属于无创伤性，操作简单，病人易接受，有着广泛的应用前景。
主要应用领域有：一，仿生学。提出肌肉生理模型来判别肌肉的动作以来, 电子假肢的研究进入了新的发展时期, 过去电子假肢的控制靠使用者人为开关和选择运动模式来完成 现在则可通过检测人体残肢表面肌电信号, 提取出肢体的动作特征, 来自动控制假肢运动, 利用残肢表面肌电信号的肌电假肢研制在国内外都取得较大进展。二，康复工程。如利用表面肌电信号提取出的特征作为功能性电刺激的控制信号, 帮助瘫痪的肢体恢复运动功能。通过检测表面肌电信号, 并将其作为反馈信号提供给病人和医生, 便于进行合理的治疗和训练。三，运动医学。表面肌电信号在运动医学中也可发挥重要作用, 通过检测运动员运动时的表面肌电信号,及时反映出肌肉的疲劳和兴奋状态, 有助于建立科学的训练方法。
本次课程设计的主要任务就是对微弱的肌电信号进行时域和频域的处理及分析，运用数字处理及matlab的知识进行“屈”和“伸”动作识别。然后通过串口将数据发送到单片机下行微机进行显示。

第二章 系统总体设计方案
根据课程设计要求在上微机利用matble分析肌电信号并处理，基于肌电信号分析结果，通过串口发送命令给单片机系统，根据肌电信号动作状态控制相应的数码管显示。并增加了扩展模块，通过动作模式驱动电机转动或其他控制输出模块。
通过分析上位机matlab中对信号处理的结果，我们可以得到一系列的信号特征值，其中我们选取了具有代表意义积分肌电值来进行处理，并给出对于屈伸动作的阈值。通过对阈值的判断，使数码管显示积分肌电值，使点阵模块显示相应的“屈”和“伸”字样，使电机根据动作进行正反转，蜂鸣器在“屈”动作是发出鸣响。
最后，将上微机的处理数据通过串口通信发送到下微机显示，得到动作识别的要求。

第三章 肌电信号的时域参数处理及其分析

均值：对于一组随机变量来说，均值是一个很重要的数值特征，用来描述一组变量的平均水平。其严格的数学定义非常简单，就是一个随机变量关于概率测度的积分。因此，在此处，均值表示肌电信号的平均水平。公式如下：

（2）方差：方差是各个数据与平均数之差的平方的平均数。

在数理统计和概率论中方差用来度量随机变量和其数学期望（即均值）之间的偏离程度。

在此处，研究信号的随机变量和均值之间的偏离程度有着很重要的意义。

其求解公式如下：

（3）标准差：标准差也称均方差，是各数据偏离

平均数的距离的平均数，它是离均差平方和平均后的方根，用u表示。标准差是方差的

算术平方根。标准差能反映一个数据集的离散程度。所以能很好的反映肌电信号的离散程度。

求解公式如下：

积分肌电值IEMG：积分肌电值就是对所有信号取绝对值后尽心均值的求解，由于对肌电信号直接求均值，均值近似为零，无法表征信号间的差异。若对肌电信号取绝对值后再进行均值运算后，均值恒大于零，因而可用于提取肌电信号的特征。公式如下：

均方根RMS：均方根就是一组数据的平方和除以数据的个数再开方，均方根是最理想的平方滤波方式的典型，让滤波更平滑，更大限度的滤掉噪声。因此，对肌电信号求均方根，可以滤除信号中的噪声，使滤波后的信号更平滑、更明显。公式如下：

（6） 原始信号的时域及上述参数值的Matlab程序

clear;
close all;

a=load('qu.txt');//('shen.txt');
N=10000;
b=0:N-1;
axis([0,10000,-1,1]);

plot(b,a);
xlabel('时间 (s)');
ylabel(

'被测变量y');
grid on;
fprintf('\n数据基本信息：\n')

printf(' 均值 = %7.5f \n',mean(a));
fprintf(

' 标准差 = %7.5f \n', sqrt(var(a)));

fprintf(' 方差= %7.5f \n', var(a));
fprintf(

' 积分肌电值IEMG = %7.5f \n', mean(abs(a)));
fprintf(' 均方根有效值RMS= %7.5f \n', sqrt(mean(a.^2)) );

屈信号数据基本信息：
均值 = 0.03502
标准差 = 0.05775
方差= 0.00334
积分肌电值IEMG = 0.05437
均方根RMS= 0.17246

伸信号数据基本信息：
均值 = 0.00337
标准差 = 0.24421
方差= 0.05964
积分肌电值IEMG = 0.12826
均方根RMS= 0.19993

第四章 肌电信号的频域处理方法及其分析

3.1 FFT分析：FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来。采样得到的数字信号FFT变换。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。而每个点的相位就是在该频率下的信号的相位。如果要要提高频率分辨率，就需要增加采样点数，也即采样

时间。频率分辨率和采样时间是倒数关系。

Matlab实现：
clear;
close

all;
a=load('qu.txt');/'shen.txt'
y=fft(a,1024); %做1024点傅立叶变换
fs=1000;
N=length(y);
mag=abs(y);
f=(0:N-1)/N*fs;
figure;
plot(f,mag);%幅频谱
xlabel('频率');
ylabel('幅值');

title('（屈/伸）肌电幅频');
figure;
plot(f,angle(y));% 相频谱
xlabel('频率');
ylabel('相位');
title('（屈/伸）肌电相频');
grid on;

3.2 功率谱分析：功率谱是信号或噪声的自相关函数的傅里叶变换。如果一确定信号平均功率为有限的，则该信号的自相关函数存在，如随机信号或随机噪声是由二阶随机平稳函数表示的，则其自相关函数存在。

为了更好得描述能量信号、功率信号，我们引入能量谱密度和功率谱密度概念。能量谱密度、功率谱密度函数表示信号的能量、功率密度随频率变化的情况。通过研究功率谱密度，可以帮助了解信号的功率分布情况，确定信号的频带等

。
信号的功率谱即上述FFT分析后，幅频值mag的平方再除以2得到。功率密度谱虽然描述了随机信号的功率在各个不同频率上的分布，但因为它仅与幅度频谱有关没有相位信息，所以从已知功率谱还难以完整地恢复原来的功率信号。

Matlab实现：
a=load('qu.txt');//('shen.txt'）
y=fft(a,1024); %做1024点傅里叶变换

fs=1000;
N=length(y);
mag=abs(y);
f=(0:N-1)/N*fs;
power1=(mag.^2)/2;%周期图法求功率谱
plot(f,power1)
xlabel('频谱');
ylabel('功率谱');
title('（屈/伸）肌电信号功率谱');

3.3 倒谱分析：倒谱是

信号的傅里叶变换谱经对数运算后再进行的傅里叶反变换或者功率谱的对数值的逆傅氏变换称为倒谱。倒频谱函数C(q)(power cepstrum)其数学表达式为：C(q) = | IF(log(s(f))) |^2

。其中，s(f)是信号s(t)的傅里叶变换，log（）为取对数，IF为逆傅里叶变换。
Matlab实现：前100点的倒谱变化比较明显，所以取前100点波形图。

a=load('qu.txt');//('shen.txt')
y=fft(a,1024); %做1024点傅里叶变换
fs=1000;
N=length(y);

mag=abs(y);
f=(0:N-1)/N*fs;
w=rceps(a);%求倒谱
n=[1:100];
plot(n,w(1:100));
xlabel('时间');
ylabel('倒谱');
title('（屈/伸）肌电信号倒谱');
grid on;

3.4 平均功率频率MPF和中值频率MF:平均功率频率是总功率除以总时间。 中值频率是各个时间段的功率的平均值。公式如下：

Matlab实现：
a=load('qu.txt');
y=fft(a,1024); %做10000点傅里叶变换
fs=1000
N=length(y);
mag=abs(y);
f=(0:N-1)/N*fs;
power=(mag.^2)/2;
ss=sum(power);
M2=0.5*ss;
df=fs/N;
M1=0.5*df*(sum(power(1:N-1))+sum(power(2:N)));
MPF=M1/M2;
MF=M2/2;

第五章 Matlab程序及GUI

MATLAB是一种面向工程和科学计算的交互式计算软件，它以矩阵运算为基础，把计算、可视化、程序设计融合到了一个简单易用的交互式工作环境中。同时由于MATLAB是一个数据分析和处理功能十分强大的工程实用软件，它的信号处理与分析工具箱为肌电信号分析提供了十分丰富的功能函数,利用这些功能函数可以快捷而又方便地完成肌电信号的处理和分析，再通过GUI界面设计，显示信号处理后的结果，方便简明。
打开GUI设计界面，根据要求添加相应的按钮，并对按钮进行定义，调整界面整体布局，将界面运行，会产生相应的程序。在程序中，寻找相应的按钮定义，在定义下添加该按钮要实现的功能函数。然后对整个程序进行调试验证。得到可视化GUI界面及功能的实现。

程序如下：
function varargout = untitled20(varargin)
% UNTITLED20 M-file for untitled20.fig
% UNTITLED20, by itself, creates a new UNTITLED20 or raises the existing
% singleton*.
%
% H = UNTITLED20 returns the handle to a new UNTITLED20 or the handle to
% the existing singleton*.
%
% UNTITLED20('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in UNTITLED20.M with the given input arguments.
%
% UNTITLED20('Property','Value',...) creates a new UNTITLED20 or raises the
% existing singleton* Starting from the left, property value pairs are
% applied to the GUI before untitled20_OpeningFunction gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to untitled20_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help untitled20

% Last Modified by GUIDE v2.5 01-Jul-2013 16:55:28

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @untitled20_OpeningFcn, ...
'gui_OutputFcn', @untitled20_OutputFcn, ...
'guiLayoutFcn', [] , ...

'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before untitled20 is made visible.
function untitled20_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to untitled20 (see VARARGIN)

% Choose default command line output for untitled20
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

UIWAIT makes untitled20 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = untitled20_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure

varargout{1} = handles.output;

% --- Executes on button press in qu.
function qu_Callback(hObject, eventdata, handles)
% hObject handle to qu (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
a=load('qu.txt');
axes(handles.axes1);
plot(a);
xlabel('时间 (s)');
ylabel('被测变量y');
title('(屈)原始信号');
grid on;
A=mean(a); % 平均值
B=std(a); % 标准差
C=Var(a); % 方差
D=mean(abs(a)); % 积分机电值
E=sqrt(mean(a.^2));%均方根有效值RMS
y=fft(a,10000); %做10000点傅里叶变换
fs=1000;
N=length(y);
mag=abs(y);
f=(0:N-1)/N*fs;
power=(mag.^2)/10000;%功率谱
ss=sum(power);
M2=0.5ss;
df=fs/N;
M1=0.5*df*(sum(power(1:N-1))+sum(power(2:N)));
F=M1/M2;
G=M2/2;
% [A,B,C,D,E,F]=index_calculation(a);
% 设置显示
A1=A*100000;
A2=round(A1);
A3=A2/100000;%round函数是对数据取整，以此限制数据的有效位
set(handles.pingjunzhi,'string',num2str(A3));

B1=B*100000;
B2=round(B1);
B3=B2/100000;
set(handles.biaozhuncha,'string',num2str(B3));
C1=C*100000;
C2=round(C1);
C3=C2/100000;
set(handles.fangcha,'string',num2str(C3));
D1=D*100000;
D2=round(D1);
D3=D2/100000;
set(handles.jifenjidianzhi,'string',num2str(D3));
E1=E*100000;
E2=round(E1);
E3=E2/100000;
set(handles.junfanggen,'string'num2str(E3));
F1=F*100000;
F2=round(F1);
F3=F2/100000;
set(handles.pingjungonglv,'string',num2str(F3));
G1=G*100000;
G2=round(G1);
G3=G2/100000;
set(handles.zhongweipinglv,'string',num2str(G3));

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over pingjunzhi.
function pingjunzhi_ButtonDownFcn(hObject, eventdata, handles)
% hObject handle to pingjunzhi (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over biaozhuncha.
function biaozhuncha_ButtonDownFcn(hObject, eventdata, handles)

% hObject handle to biaozhuncha (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over jifenjidianzhi.
function jifenjidianzhi_ButtonDownFcn(hObject, eventdata, handles)
% hObject handle to jifenjidianzhi (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over pingjungonglv.
function pingjungonglv_ButtonDownFcn(hObject, eventdata, handles)
% hObject handle to pingjungonglv (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over fangcha.
function fangcha_ButtonDownFcn(hObject, eventdata, handles)
% hObject handle to fangcha (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% --- Executes on button press in shen.
function shen_Callback(hObject, eventdata, handles)
% hObject handle to shen (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
a=load('shen.txt');
axes(handles.axes2);
plot(a);

xlabel('时间 (s)');
ylabel('被测变量y');
title('（伸）原始信号');
grid on;
A=mean(a); % 平均值
B=std(a); % 标准差
C=Var(a); % 方差
D=mean(abs(a)); % 积分机电值
E=sqrt(mean(a.^2));%均方根有效值RMS
y=fft(a,10000); %做10000点傅里叶变换
fs=1000;
N=length(y);
mag=abs(y);
f=(0:N-1)/N*fs;
power=(mag.^2)/10000;%周期图法功率谱
ss=sum(power);
M2=0.5*ss;
df=fs/N;
M1=0.5*df*(sum(power(1:N-1))+sum(power(2:N)));
F=M1/M2;
G=M2/2;
% [A,B,C,D,E,F]=index_calculation(a);
% 设置显示
A1=A*100000;
A2=round(A1);
A3=A2/100000;
set(handles.pingjunzhi,'string',num2str(A3));
B1=B*100000;
B2=round(B1);
B3=B2/100000;
set(handles.biaozhuncha,'string',num2str(B3));
C1=C*100000;
C2=round(C1);
C3=C2/100000;
set(handles.fangcha,'string',num2str(C3));
D1=D*100000;
D2=round(D1);
D3=D2/100000;
set(handles.jifenjidianzhi,'string',num2str(D3));

E1=E*100000;
E2=round(E1);
E3=E2/100000;
set(handles.junfanggen,'string',num2str(E3));
F1=F*100000;
F2=round(F1);
F3=F2/100000;
set(handles.pingjungonglv,'string',num2str(F3));
G1=G*100000;
G2=round(G1);
G3=G2/100000