手部动作的肌电信号分析与建模
手部动作的肌电信号分析与建模
摘要
随着人机交互技术的发展，机器的智能程度不断提高，利用肌肤表面的电位信号
判断手部动作就是典型的智能技术，其具有广泛的应用前景和实用价值。本文根据实
验一、实验二中所给的数据进行分析，对每种动作所反馈的不同识别信号进行研究并
建出模型。
1、针对问题一，在MATLAB 中利用plot 绘图函数将所得数据绘制肌电信号图（图
5.1—5.6），比较各个手部动作的肌电信号之间的差异，从图像可以看出各个手部动
作的肌电信号的差异性比较明显，总结差异。这说明建立的模型的准确性、可靠性。
2、针对问题二，进行时域、频域、小波变换、预处理分析，利用MATLAB 软件计
算肌电信号的均值、标准差、方差、积分肌电值IEMG 、均方根有效值RMS 等时域指标
参数和平均功率频率MPF 、中值频率MF 等频域指标，得到时域、频域、小波变换图像，
得出各个手部动作的肌电信号的特征数值。
3、针对问题三和问题四，创建AR 识别模型，在MATLAB 中利用plot 函数将剩余
数据绘制肌电信号图（图5.15—5.20），六种动作的“能量-功率”特征差异是非常
明显的，从而准确的识别出六种动作。用最小二乘拟合算法来检验识别模型的好坏，
得到的拟合曲线几乎完全重合，拟合效果不明显。再次对数据进行若干次小波分解，
将数据简单化、特征化，计算出每层分解的标准差σi, 到达八级小波分解时，数据
趋于固定，得到无偏估计值图表（表5.2—5.3），然后对其做lsqcurvefi 散点-曲线
拟合处理，并将拟合图像叠加，得到图像（图5.23—5.24），显然所得图像和参数范
围中各个动作差异明显，不同动作的识别准确率有一定差异。针对模型阶数、编程、
噪音、测试仪器等因素，分析影响识别准确率的原因。
4、针对问题五，模型建立方法与识别方法与问题三、四是相同的，关键在于如
何优化第三天的识别率，我们综合考虑的可能影响识别率的模型、噪声、时间、预处
理等因素，对其做最大化的单项优化，结合后即为系统最佳优化方案。
关键词：肌电信号识别；时-频域分析；AR 模型；最小二乘拟合算法；无偏估计
一、 问题重述
随着人机交互技术的发展，机器的智能程度不断提高，肌电生理研究的深入和肌
电检测技术也达到了一个前所未有的高度，其中利用肌肤表面的电位信号判断手部动
作就是典型的智能技术。表面肌电是从人体皮肤表面通过电极记录下来的神经肌肉活
动时发放的生物电信号，属于无创伤性，操作简单，病人易接受，在基础医学研究、
临床诊断和康复工程中有着广泛的应用前景。
当检测信号与相应的手部动作不符时，会使人产生错误的判断，因而能否根据电
信号准确识别所表达的动作成为此项技术的关键。为了肌电信号识别模型的建立，我
们分析题目给出的实验数据，利用这些数据解决以下问题：
1. 利用MATLAB 中的load 函数，载入数据到内存，利用绘制函数绘制出信号，
观察各种动作之间的差异，总结差异在哪些方面
2. 利用合适的时域或（和）频域特征表达肌电信号，建立特征向量，写出提取
特征向量的具体方法和程序代码。
3. 针对Database 1的female_1人，利用6个动作的部分数据建立识别模型，
然后利用剩余的数据识别信号对应的6种的不同动作类别，分析影响识别准
确率的原因。
4. 针对Database 1中的所有人，利用female_1、female_2和male_1的所有数
据建立识别模型，然后利用female_3和male_2的所有数据检验识别模型的
好坏，即识别信号对应的6种的不同动作类别，分析影响识别准确率的原
因。
5. 为了分析不同时间测量对识别动作的准确性的影响，利用Database 2中的
male_day_1的所有数据建立识别模型，然后分别利用male_day_2和
male_day_3的所有数据检验识别模型的好坏，分析两天识别准确率的差异。
考虑如何提高第3天的识别准确率。
二、 问题分析
在该题目中我们需要将两组实验所测肌电信号进行分析与建模，从而利用模型判
断所测信号对应的手部动作，需要考虑性别、肌肉疲劳等对所测信号的影响。问题的
关键在于如何排除一些关系不大的因素，获取有用信息，找出肌电信号与手部动作的
关系。
2.1不同手部动作肌电信号差异分析
问题要求观察不同动作之间的信号差异，为结果更准确，我们采取单一变量法，
即分析某一个人不同手部动作的肌电信号。为了防止奇异数据影响正常结果，我们将
六种手部动作的30次重复实验所得双通道数据分别取均值，在MATLAB 中利用plot
绘图函数将所得数据绘制肌电信号图，做出相应观察与比对。 2.2肌电信号的特征向量的建立及对特征向量的提取方法分析
利用MATLAB 软件对肌电信号进行时域和频域的特征建立，计算均值、标准差、
方差、积分肌电值IEMG 、均方根有效值RMS 等时域指标参数和平均功率频率MPF 、中
值频率MF 等频域指标，得到时域、频域图像，比较不同动作下的肌电信号特征。从
特征完善性分析，我们重新对数据进行小波变换处理，进行多次分解，得到分解后
的特征图像。
2.3六种手部动作识别模型的建立及识别准确率分析
根据实际肌电信号的随机性特征, 基于AR(Autoregressive)模型, 引入能量、功率
概念，做出了“能量-功率”识别图像，得到AR 模型的各项参数, 分析此参数和对应
肌肉活动所确定的肢体动作之间的关系, 从而得到基于动作模式的表面肌电信号
(EMG)AR模型参数。
在识别准确率方面，针对问题3，利用MATLAB 将female-1六种手部动作原始数
据的部分数据做出模型，做出模型的能量-功率图像进行识别；针对问题4，利用
female_1、female_2和male_1的所有数据建立识别模型，然后利用female_3和
male_2的所有数据与其做拟合曲线，从拟合度分析模型的识别准确率。另外，原始数
据预处理程度、AR 模型性质两个方面的因素会影响识别准确率。
2.4提高模型的识别准确率方法分析
模型识别准确率受很多因素影响，包括噪声，模型局限性，傅立叶变换及均方
误差等问题，因此在实验及数据处理上注意减小或消除误差，可提高模型的识别准确
率。
三、 模型假设
（1） 假设各样本数据能真实客观地反映电信号与手部动作的关系；
（2） 假设每位受试者做重复动作时动作是完全一致的；
（3） 假设各动作所产生电信号不因男女性别不同而产生差异；
（4） 假设所选的时频域能够正确地表达出肌电信号；
（5） 假设在MATLAB 中影响不同手部动作肌电信号的因素已考虑完全。
四、 符号说明
五、 模型的建立与求解
5.1
问题1的模型建立与求解
5.1.1手部动作肌电信号分析处理
在采用单一变量法原则后，对于六种动作、每种动作两个通道的数据成像仍十分
繁杂，考虑到采取肌电信号的两个传感器均在前臂皮肤上，在同种动作下两个位置肌
肉所表达的能量总体特征是一致的，我们即用六种动作单一通道的数据成像来分析不
同动作肌电信号差异
。
5.1.2用plot （）绘图函数绘制肌电信号图
在MATLAB 环境中，用plot （）函数将处理后的数据进行绘图，得到了六种手部
动作相应肌电信号图形。
图5.1 握圆柱物体 图5.2 提手柄
图5.3 夹片状物体 图5.4 手执长物体
图5.5 握球 图5.6 掐小颗粒物体
5.1.3图像差异分析
通过以上6副图可以判断，在不同手部动作下，测得的肌电信号的幅值是不同
的，并且与手形具有一定的关联：比如掐、执、夹时手部靠食指与拇指发力，其图像
幅值的Max 在0.3附近；而握、提时发力的不仅仅是两只手指，所以图像幅值的Max
要大，在1附近波动。
5.2问题2的模型建立与求解
5.2.1时域特征
时域特征的提取方法是将肌电信号看成均值为零、方差随信号强度变化而变化的
随机信号，由于时域特征的提取比较简单，故在肌电信号应用领域得到比较广泛的应
用，最常用的肌电信号的时域特征有以下几种：
均值：
对于一组随机变量来说，均值是一个很重要的数值特征，用来描述一组变量的平
均水平。其严格的数学定义非常简单，就是一个随机变量关于概率测度的积分。
因此，在此处，均值表示肌电信号的平均水平。公式如下
方差：
方差是各个数据与平均数之差的平方的平均数。在概率论和数理统计中，方差（英
文Variance ）用来度量随机变量和其数学期望（即均值）之间的偏离程度。在许
多实际问题中，研究随机变量和均值之间的偏离程度有着很重要的意义。其求解
公式如下：
积分肌电值IEMG ：
积分肌电值就是对所有信号取绝对值后尽心均值的求解，由于对肌电信号直接求
均值，均值近似为零，无法表征信号间的差异。若对肌电信号取绝对值后再进行
均值运算后，均值恒大于零，因而可用于提取肌电信号的特征。公式如下：
均方根RMS ：
均方根就是一组数据的平方和除以数据的个数再开方，均方根是最理想的
平方滤波方式的典型，让滤波更平滑，更大限度的滤掉噪声。因此，对肌电信号
求均方根，可以滤除信号中的噪声，使滤波后的信号更平滑、更明显。公式如下：
5.2.2频域特征
FFT 分析方法：
FFT 是离散傅立叶变换的快速算法，可以将一个信号变换到频域。当某些信号特
征在时域表现不明显时，变换到频域之后有可能就很容易看出特征了，这就是很多信
号分析采用FFT 变换的原因，并且在频谱分析方面FFT 可以将一个信号的频谱提取出
来。为了方便进行FFT 运算，通常N 取2的整数次方。
实验中采样频率为Fs=500Hz，设信号频率F ，采样点数为N 。那么FFT 之后结果
就是一个为N 点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率
值下的幅度特性。假设原始信号的峰值为A ，那么FFT 的结果的每个点（除了第一个
点直流分量之外）的模值就是A 的N/2倍。而第一个点就是直流分量，它的模值就是
直流分量的N 倍。而每个点的相位就是在该频率下的信号的相位。如果要要提高频率
分辨率，就需要增加采样点数，也即采样时间。
功率谱分析方法：
定义信号f(t)的能量（作归一化处理）：由电压f(t)（或者电流f(t)）在1Ω电
阻上消耗的能量
E =⎰f 2(t ) dt -∞∞
其中 E =u ⋅i =u 2
=u 2，若积分值存在，信号的能量为有限值，则称f(t)为能量信
号。
对于能量无限大的信号（如周期信号），我们考虑能量的时间平均值，这显然就是信号的平均功率。这种信号称为（平均）功率信号。
定义信号f(t)的平均功率：电压f(t)在1Ω电阻上消耗的平均功率（简称功率） P =lim T →∞122T f (t ) dt ⎰-T 2T
式中，T 是为求平均的时间区间。
在这两个理论模型的基础上，我们通过MATLAB 做出了六种动作的时频域图像，这里我们仅拿出握圆柱体与提手柄的图像做参考：
图5.7 握圆柱物体 图5.8 提手柄
平均功率频率MPF 和中值频率MF:平均功率频率是总功率除以总时间。 中值频率是各个时间段的功率的平均值。公式如下：
倒谱分析法：功率谱的对数值的逆傅氏变换称为倒谱。信号进行3000点短时傅
ˆ定义如下：里叶变换，经对数运算后再进行的傅里叶反变换及加窗处理。如序列倒谱c ˆ=F -1[log F [x )] c
其中F 和F -1分别表示傅里叶变换和逆傅里叶变换。计算LPC 系数，建立肌电信号的正则方程，对倒谱进行预测误差。在MATLAB 中编程实现6种肌电信号倒谱图的绘制，我们取握圆柱物体与夹片状物体的图像做参考：
图5.9 握圆柱物体 图5.10 夹片状物体
5.2.3小波变换提取数据特征：
功率谱中的频率成份表征出了信号的重要信息，但是却没有反映出这些频率成分对应的时域信息。小波分析不同于傅里叶变换，它对时域和频域信息均有很好的反映，是应用较为广泛的一种特征提取方法。
离散小波变换可以看成是信号的层层分解，首先信号被一个个像滤波器组g(n)和h(n)进行高通、低通滤波，滤波器的截止频率均为Fs/4，继而信号被二倍下采样，由此得到信号的近似分量cA 和细节分量cD 。第二层对第一层的近似分量再进行类似分解, 滤波器组的截止频率变为Fs/8，层层分解直到指定阶段。以三级分解为例，如图5.13，离散信号经尺度j=1,2,3的小波分解，得到逼近系数cA3和小波系数cD3，cD2，cD1。
图5.11
我们通过MATLAB 对6种手部动作肌电信号分别进行5次小波变换，确定阈值和去噪。其中采用启发式（heusure ）和极大极小原理（minimaxi ）得到阈值，去噪则采用一维信号的自动消噪，以此获取特征图像，我们取握球、握圆柱物体的特征图像做参考：
图5.12 握球 图5.13 握圆柱物体
5.2.4预处理分析
因现有表面肌电信号采集结构，肌电信号在测量、记录阶段都存在一定的噪声影响（见图5.14），包括环境噪声、运动产生噪声、电极噪声等，采取必要的措施对噪声加以抑制，提高肌电信号的信噪比，是增强后续动作识别效果的必要措施。在进行肌电信号特征向量的建立之前，我们对所有实验数据进行了小波去噪，一个含噪声的一维信号的模型可以表示为如下形式：
其中f(i)为真实信号，e （i ）为噪声，s （i ）为含噪声的信号。
图5.14表面肌电信号采集结构图
5.2.5基于以上分析方法所得六种不同手部动作各特征向量数值
我们根据题目实验数据，利用MATLAB 软件对肌电信号进行时域和频域(包括FFT ，功率谱，倒谱) 分析，计算均值、标准差、方差、积分肌电值IEMG 、均方根有效值RMS 等时域指标参数和平均功率频率MPF 、中值频率MF 等频域指标，具体见下表：
表5.1 通道1的肌电信号各项特征数值
5.3问题3、4的模型的建立与求解
5.3.1创建基于AR 模型的能量-功率识别图像
AR 模型（自回归模型，Auto Regression Model ）是典型的现代参数功模型，是一个线性的、二阶矩平稳模型，比较适合短数据分析，而且运算便利，应用于表面肌电信号的识别。AR 模型可以表示如下：
其中ai 为AR 系数，n 为模型阶次，k 表示肢体动作，M 为动作数，Wk 代表观测信号与拟合数据的误差，当Wk 为白噪声时模型线性最优。
AR 模型的Yule-Walker 方程和Levinson-Durbin 递推算法：在MATLAB 中，函数levinson 和aryule 都采用Levinson-Durbin 递推算法来求解AR 模型的参数a1,a2,„„,ap及Wk 白噪声序列的方差。
为了更好得描述能量信号、功率信号，我们引入能量谱密度和功率谱密度概念。能量谱密度、功率谱密度函数表示信号的能量、功率密度随频率变化的情况。
对于肌电信号而言我们是将每次采集时间定在一定时间段内，因此我们对过程的样本函数应用截取函数。设过程ξ(t )的截取函数ξT (t )（截取的随机过程）为：
⎧⎪ξ(t ) ξT (t )=⎨⎪⎩0则截取函数的傅里叶变换为：
t ≤
T
2 其他
ξT (ω) =⎰ξT (t ) e
-∞
∞
-j ωt
dt =⎰ξT (t ) ⋅e -j ωt d ω
T 2T -2
平稳随机过程ξT (ω) 的平均功率为
⎧1S ξ=E ⎨lim
⎩T →∞T
⎰
T 2T -2
ξ2(t ) dt ⎬=
⎭
⎫
12π
⎰
∞
-∞
P ξ(ω) ⋅d ω
ξ(t )的功率谱密度为P ξ(ω)
P ξ(ω) =lim
T →∞
E [ξT (ω) ]
T
2
这样的平均功率等于各个频率分量（统计值）单独贡献出的功率之连续和，是在频率域上描述随机过程统计特性的最主要数字特征。
基于AR 模型的原理，我们计算真实的自相关值时，采用逆Levinson-Durbin 递归方法，利用Matla 求得滤波器的阶数、各参数值，再计算他们自相关值。将各个手部动作肌电信号数据进行3000点快速傅立叶变换，其中我们引入了消去直流分量的概念，使功率谱更能体现能量、功率等有效信息，以图像的形式表达出各个手部动作特征，因而从图像模型中我们可以直观的识别其对应哪个手部动作。
5.3.2六种手部动作AR 模型识别
基于AR 模型与“能量-功率”图像成像特征，在MATLAB 环境下我们分别得到了六种手部动作图像，具体如下：
图5.15 握圆柱物体 图5.16 提手柄
图5.17夹片状物体 图5.18 手执长物体
图5.19 握球 图5.20掐小颗粒物体
显然这六种动作的“能量-功率”特征差异是非常明显的，因此对其余数据进行手部动作识别时，通过其成像特征就可以识别出该肌电信号对应的手部动作。
5.3.3模型好坏分析
数据拟合：
在5.3.2中，我们得到了关于六种手部动作的AR 识别模型，但是其模型好坏及识别率并未分析，下面进行的就是针对这部分的分析。
首先，我们将female_1、female_2、male_1与female_3、male_2的所有原始数据进行lsqcurvefi 散点-曲线拟合处理，然后以握圆柱物体（cyl ）为例，建立相应功能函数fun （内容见附录），所采用拟合方程为：
f=1-8/9.8696.*exp(-9.8696.*c(1).*x/(c(2).^2))
得到第一幅图拟合参数 c(1)=-13.4730，c(2)=97.1325； 得到第二幅的拟合参数 c(1)=-13.4786, c(2)=97.4249；
图5.21中第一幅散点图为female_1、female_2、male_1的数据拟合曲线，第二幅为female_3、male_2的数据拟合曲线，第三幅则为将两者做叠加对比的数据拟合曲线。把第三幅拟合曲线放大后即图5.22，我们发现两条曲线几乎完全重合，这显然是不符合常理的。用MATLAB 中toolbox 自带的cftool 拟合函数进行重复实验验证，得到相同结果，证明我们所采用的拟合程序并没有错误。
分析原因：
我们所采用的数据量较大，散点分布状况特征相似，得到的拟合曲线参数极为相似，故得到的拟合曲线几乎完全重合。
图5.21 lsqcurvefi散点-曲线拟合图像
图5.22 放大后的叠加对比散点-曲线拟合图像
解决办法：
我们引入概率统计的有关知识，在有限次测量时，用其算术平均值X mean 来代替总体均值（真值）μ，理论和实践均证明算术平均值X mean 最接近真值，且n 越大越接近真值，X mean 是μ的无偏估计，用X mean 代替μ就可得到σ的常用估计式：
σ是表征测量值分散特性的一个度量指标，σ越大，样本落在真值附近的概率越小，即样本分散；反之越集中。
其次我们再对数据进行若干次小波分解，将数据简单化、特征化，计算出每层分解的标准差σ经实验反复尝试，当分解层数较少时（小于3层），其标准差差异明显，而当到达八级小波分解时，其标准差趋于固定，故我们采用八级小波分解后的数据，计算其标准差，结果如下表：
i ，
表5.2 female_1、female_2、male_1的8次小波分解后的无偏估计
表5.3 female_3、male_2的8次小波分解后的无偏估计
对female_1、female_2、male_1与female_3、male_2的实验数据分别进行8次小波分解得到的无偏估计值，然后对其做lsqcurvefi 散点-曲线拟合处理，并将拟合图像叠加，得到图像，下面以cyl 与hook 为例：
图5.23 cyl的无偏估计拟合曲线
图5.24 hook的无偏估计拟合曲线
模型评价：
显然我们所得图像中两种动作（cyl 与hook ）的特征参数差异明显，相同动作的两条曲线具有相似的特征，说明我们建立的识别模型具有一定的识别能力，可以根据不同的肌电信号找到其对应的手部动作。但cyl 的拟合程度较hook 差，说明不同动作的识别准确率有一定差异，且差异发生在不同的范围内。
5.3.4影响识别准确率的原因分析
①表面肌电信号十分微弱，易受环境噪声、运动产生噪声、电极噪声等影响，而小波去噪处理并不能完全消除噪声。
②AR 模型的阶数n 选取不合理，n 太低，谱估计的分辨率不够；太高则会出现许多实际不存在的虚假细节。
③傅立叶变换本身具有一定的局限性，即仅适用于平稳的非突变信号及变换伴有位置信息的丢失
5.4问题5的模型建立与求解 5.4.1模型建立、求解与评价
利用Database 2中的male_day_1的所有数据建立AR 识别模型，模型与之前的模型一样，只是改变模型取值域与参数，这里我们仍旧分析cyl 的三天的肌电信号图像差异。
通过MATLAB
编程实现的功率图像如下：
图5.25 male的cyl 三天功率谱图
由图像直观得到与第一天相比，随着天数增加，其功率图像的差异越大。
接着对cyl 的3天的原始数据进行8次小波变换，得到对应的无偏估计值，进行lsqcurvefi 散点-曲线拟合处理，得到图5.26如下
图5.26 male的cyl 三天拟合曲线图
5.4.2识别准确率优化问题
①针对AR 模型，为更加有效地实现模式分类，有些研究将AR 模型系数变换为反射系数、倒谱系数和对数面积比等新的参数，将其作为特征适量以提高聚类分离度。许多实验表明，AR 模型的阶数取4时，对信号的分析和识别性能是最好的，阶次更高不会改善分类结果，只会增大运算量。
②在功率谱分析上，为了解决频率分辨率低与频谱能量向旁侧泄露的问题，提出一系列的周期图法的改进方法或最大似然法实现高分辨率功率谱估计。
③在数据预处理上，不仅使用小波去噪，而且在小波去噪的阈值选取采用自适应无偏估计，最后使用加窗处理使数据更具科学性。
④将预处理后的肌电信号，采用SOM 、LVO 、BP 、人工神经网络等网络进行识别，采用对多的是SOM 网络，其他网络的学习、训练建模过程相对复杂，且AR 模型进行识别已有足够的准确度，故本文不予讨论。
六、 模型评价
6.1模型的优点和缺点 6.1.1模型的优点
（1) 我们基于AR 模型构建的“能量-功率”图像可以直观地反应各个手部动作之间的差异，对于表面肌电信号的识别也更加方便。
（2）调用MATLAB 提取的肌电动作信号清晰，对比性强；plot()函数绘制的声音信号图直观、易懂。
（3）数据处理后肌电信号频谱的稳定性。
（4）AR 模型是一个线性的、二阶矩平稳模型, 比较适合短数据分析, 而且运算便利, 特别适合肌电控制假肢的实时处理。
（5）运用小波去噪，加窗来优化参数，极大的提高了识别的性能和正确率。
（6）利用原始数据的部分数据做出模型，然后用其余数据对其做拟合曲线，从拟合
度分析模型的识别准确率，提高了问题的直观性。 6.1.2模型的缺点
（1）应用MATLAB 时，不能把完全影响正常和非正常的表面肌电信号因素考虑进去。 （2）时域方法最早应用于肌电信号分析，但其不稳定、易受干扰等不足限制了它在肌电信号处理领域的应用。
（3）由于肌电信号在本质上是非平稳的, 统计特性随时间变化, 各种非平稳信号模型、自适应处理和时变系统辨识方法也被纷纷用于肌电研究；因此AR 模型对于肌电信号的非平稳特性具有局限性。
（4）预处理过程中不能完全去噪。 6.2模型的意义
以小波变换为代表的时-频分析方法因结合了时域、频域两方法的特性，在肌电信号分析方面颇有潜力，正受到越来越多研究者的青睐；针对肌电信号的非线性，现已有学者采用非线性动力学方法对其进行分析，取得了一些成果，但仍处于起步阶段，有待更深入的研究．肌电信号的复杂性，决定了单独使用某类方法进行肌电信号处理都可能无法充分利用所得信息，而计算机技术、信息处理技术等的飞速发展，为综合应用各种技术进行肌电信号的处理提供了可能性，也为肌电信号研究的进一步发展奠定了基础。
6.3模型的推广
肌电信号的分析与识别对于康复医学和临床诊断都具有深远的意义。随着计算机技术和信号处理技术的不断发展 , 采用多元时间序列分析、人工神经网络及维格纳分布、小波变换等现代分析方法对肌电信号进行研究, 有可能成为可靠有效地提取信号特征, 实现肌电模式的识别, 从而为肌电控制假肢、功能性神经电刺激的发展以及人体神经系统疾病的诊断提供更加新颖、有效的手段。肌电分析被广泛地应用于肌肉疲劳、重症肌无力、肌强直和肌萎缩等各种肌肉疾病的临床诊断；肌电图也是运动员训练中动作和强度分析的依据；另外，还可利用人体表面肌电的某些特征进行模式分类，进而驱动假肢或其它运动机械(如机器人) 的各种动作，实现仿生控制．由此可见，提取有效的肌电信号特征对肌电信号分析意义重大。随着检测技术、信息处理技术及计算机技术的发展，使肌电信号的定量分析及更深入的研究成为可能。目前，肌电信号已经深入应用到临床医学、运动医学、生物医学工程等领域。
参考文献
[1]章绍辉，数学建模，北京：科学出版社，2010年
[2]丁玉美，阔永红，数字信号处理，西安：西安电子科技大学出版社，2002年 [3]王丹力，MATLAB 控制系统，北京：中国电力出版社，2007年
[4]韩中庚，数学建模竞赛——获奖论文精选与点评，北京：科学出版社，2007年 [5]蔡立羽，王志中，肌电信号分析方法的研究及进展，上海交通大学，（200030） [6]罗志增，杨广映，表面肌电信号的AR 参数模型分析方法，
http://wenku.baidu.com/link?url=XWxPDtyEaR3QpBkicZMXakBIqs0oafZ1KRSJp940s0QwAE_1kueoWtAqbCYwAmnLTK9J2m2acJ51JUFkMyt3lohRClLpJbDaRb7geobJPRa，2015年5月31日
[7]丛玉良，王志宏，数字信号处理原理，北京：电子工业出版社，2009年 [8]同济大学应用数学系编，高等数学，北京：高等教育出版社，2006年
附录
注：所有解决不同数据问题的源代码直接改动load 的提取地址
问题1：
clear;clc;
load('E:\A\Database 1\female_1.mat')
arv_cyl_ch2=mean(cyl_ch2,1) %按列平均，即重复30次的平均数值在6秒的平均数值 plot(arv_cyl_ch2);grid on;
set(gca,'xticklabel','0|1|2|3|4|5|6'); xlabel('时间/s'); ylabel('数字信号');
title('female1 cyl ch2');
load('E:\A\Database 1\female_1.mat')
arv_hook_ch2=mean(hook_ch2,1) %按列平均，即重复30次的平均数值在6秒的平均数值 plot(arv_hook_ch2);grid on;
set(gca,'xticklabel','0|1|2|3|4|5|6'); xlabel('时间/s'); ylabel('数字信号');
title('female1 hook ch2');
load('E:\A\Database 1\female_1.mat')
arv_lat_ch2=mean(lat_ch2,1) %按列平均，即重复30次的平均数值在6秒的平均数值 plot(arv_lat_ch2);grid on;
set(gca,'xticklabel','0|1|2|3|4|5|6'); xlabel('时间/s'); ylabel('数字信号');
title('female1 lat ch2');
load('E:\A\Database 1\female_1.mat')
arv_palm_ch2=mean(palm_ch2,1) %按列平均，即重复30次的平均数值在6秒的平均数值 plot(arv_palm_ch2);grid on;
set(gca,'xticklabel','0|1|2|3|4|5|6'); xlabel('时间/s'); ylabel('数字信号');
title('female1 palm ch2');
load('E:\A\Database 1\female_1.mat')
arv_spher_ch2=mean(spher_ch2,1) %按列平均，即重复30次的平均数值在6秒的平均数值 plot(arv_spher_ch2);grid on;
set(gca,'xticklabel','0|1|2|3|4|5|6'); xlabel('时间/s'); ylabel('数字信号');
title('female1 spher ch2');
load('E:\A\Database 1\female_1.mat')
arv_tip_ch2=mean(tip_ch2,1) %按列平均，即重复30次的平均数值在6秒的平均数值 plot(arv_tip_ch2);grid on;
set(gca,'xticklabel','0|1|2|3|4|5|6');
xlabel('时间/s');
ylabel('数字信号');
title('female1 tip ch2');
问题2：
%FFT变换, 获得采样数据基本信息，时域图，频域图，进行数据归一化处理
clear;
close all;
load 'E:\A\Database 1\female_1.mat' %通过仪器测量的原始数据
A=mean(cyl_ch1,1); %求平均值
B=[0.002:0.002:6]; %构建时间序列
C=[B;A]; %构建时间数值序列
x=C(1,:); %将C 中的第一列赋值给x ，形成时间序列 y=C(2,:); %将C 中的第二列赋值给y ，形成被测量序列
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf(' 采样点数 = %7.0f \n',length(x)) %输出采样数据个数
fprintf(' 采样时间 = %7.3f s\n',max(x)-min(x)) %输出采样耗时 fprintf(' 采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x))) %输出采样频率 fprintf(' 最小数值 = %7.3f \n',min(y)) %输出本次采样被测量最小值
fprintf(' 平均数值 = %7.3f \n',mean(y)) %输出本次采样被测量平均值
fprintf(' 最大数值 = %7.3f \n',max(y)) %输出本次采样被测量最大值
fprintf(' 标准方差 = %7.3f \n',std(y)) %输出本次采样数据标准差
fprintf(' 协 方 差 = %7.3f \n',cov(y)) %输出本次采样数据协方差
fprintf(' 积分肌电值IEMG = %7.5f \n', mean(abs(y)));
fprintf(' 均方根有效值RMS= %7.5f \n', sqrt(mean(y.^2)) );
%显示原始数据曲线图（时域）
subplot(4,1,1);
plot(x,y); %显示原始数据曲线图
set(gca,'xticklabel','0|1|2|3|4|5|6');
xlabel('时间 (s)');
ylabel('被测变量y');
title('原始信号(时域图)');
grid on;
%快速傅立叶变换
y=y-mean(y); %消去直流分量，使频谱更能体现有效信息
Fs=500; %得到原始数据female_1.mat时，仪器的采样频率。其实就是length(x)/(max(x)-min(x));
N=3000; %female_1.mat中的被测量个数，即采样个数。其实就是length(y);
z=fft(y,3000); %做10000点傅里叶变换
%频谱分析
f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N; %幅值，单位同被测变量y
Pyy=Mag.^2; %能量；对实数系列X ，有
X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
%显示频谱图(频域)
subplot(4,1,2)
plot(f(1:N/2),Mag(1:N/2),'r'); %显示频谱图
xlabel('频率 (Hz)')
ylabel('幅值')
title('频谱图(频域图)')
grid on;
%显示能量图(频域图)
subplot(4,1,3)
plot(f(1:N/2),Pyy(1:N/2),'y');
xlabel('频率 (Hz)')
ylabel('能量')
title('能量图(频域图)')
grid on;
%求功率谱
Power=Pyy/N;
subplot(4,1,4);
plot(f(1:N/2),Power(1:N/2),'b');
xlabel('频率(Hz)');
ylabel('功率');
title('功率谱图(频域图');
grid;
%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N));
fprintf('\n傅立叶变换结果：\n')
fprintf(' FFT_f = %1.3f Hz\n',f(b)) %输出最大值对应的频率 fprintf(' FFT_T = %1.3f s\n',1/f(b)) %输出最大值对应的周期
%平均功率频率MPF 和中值频率MF
load 'E:\A\Database 1\female_1.mat' %通过仪器测量的原始数据
A=cyl_ch1; %将测量数据赋给A ，A 可以为其他数据 A=mean(cyl_ch1,1);
B=[0.002:0.002:6]; %构建时间序列
C=[B;A]; %构建时间数值序列
x=C(1,:); %将C 中的第一列赋值给x ，形成时间序列 y=C(2,:); %将C 中的第二列赋值给y ，形成被测量序列 z=fft(y,3000); %做10000点傅里叶变换
Fs=500;
N=length(z);
Mag=2*abs(z)/N;
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; %平均功率频率MPF
MF=M2/2; %中值频率MF
fprintf('\n平均功率频率MPF 和中值频率MF 计算结果：\n')
fprintf(' 平均功率频率MPF = %1.3f \n',MPF) %输出平均功率频率MPF fprintf(' 中值频率MF = %1.3f \n',MF) %输出中值频率MF
倒谱分析
clear;
close all;
load 'E:\A\Database 1\female_1.mat'
A=cyl_ch1; %将测量数据赋给A ，A 可以为其他数据 A=mean(cyl_ch1,1);
plot(A);
title('原始数据');
%对指定帧位置进行加窗处理
Q = A'; N = 64; % 窗长
Hamm = hamming(N); % 加窗
frame = 60;%需要处理的帧位置
M = Q(((frame - 1) * (N / 2) + 1):((frame - 1) * (N / 2) + N));
Frame = M .* Hamm;%加窗后的肌电信号谱
[B,F,T] = specgram(A,N,N/2,N);
[m,n] = size(B);
for i = 1:m
FTframe1(i) = B(i,frame);
end
ai = lpc(Frame,10);
% 计算lpc 系数
LP = filter([0 -ai(2:end)],1,Frame); % 建立肌电信号的正则方程
FFTlp = fft(LP);
E = Frame - LP; % 预测误差
subplot(3,1,1);
plot(1:N,Frame,1:N,LP,'-r');
grid on;
title('原始数据和预测数据')
subplot(3,1,2);
plot(E);
grid on;
title('预测误差');
y=fft(A,3000); %做3000点傅里叶变换
Fs=500;
N=3000;
f=(0:N-1)/N*Fs;
w=rceps(A);%求倒谱
n=[1:100];
subplot(3,1,3);
plot(n,w(1:100));
xlabel('时间');
ylabel('倒谱');
title('肌电信号倒谱');
grid on;
小波分析
load 'E:\A\Database 1\female_1.mat' %通过仪器测量的原始数据
A=cyl_ch1; %将测量数据赋给A ，A 可以为其他数据 A=mean(cyl_ch1,1);
B=[0.002:0.002:6]; %构建时间序列
subplot(311);
plot(B,A);
xlabel('采样时间');
ylabel('female cly ch1原始信号幅值');
grid on;
s1=wden(A,'minimaxi','s','one',5,'db3');
subplot(312);
plot(B,s1);
xlabel('采样时间');
ylabel('db3 小波变换后的信号幅');
grid on;
s2=wden(A,'heursure','s','one',5,'sym8');
subplot(313);
plot(B,s2);
xlabel('采样时间');
ylabel('sym8小波变换后的信号幅');
grid on;
问题3：
%FFT变换, 获得采样数据的能量图和功率图
clear;
close all;
load 'E:\A\Database 1\female_1.mat' %通过仪器测量的原始数据
A=cyl_ch1; %将测量数据赋给A ，A 可以为其他数据 A=mean(cyl_ch1,1); %求平均值
B=[0.002:0.002:6]; %构建时间序列
C=[B;A]; %构建时间数值序列
x=C(1,:); %将C 中的第一列赋值给x ，形成时间序列 y=C(2,:); %将C 中的第二列赋值给y ，形成被测量序列
%快速傅立叶变换
y=y-mean(y); %消去直流分量，使频谱更能体现有效信息
Fs=500; %得到原始数据female_1.mat时，仪器的采样频率。其实就是length(x)/(max(x)-min(x));
N=3000; %female_1.mat中的被测量个数，即采样个数。其实就是length(y);
z=fft(y,3000); %做10000点傅里叶变换
%频谱分析
f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N; %幅值，单位同被测变量y
Pyy=Mag.^2; %能量；对实数系列X ，有
X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
%显示能量图(频域图)
subplot(2,1,1)
plot(f(1:N/2),Pyy(1:N/2),'r');
xlabel('频率 (Hz)')
ylabel('能量')
title('能量图(频域图)')
grid on;
%求功率谱
Power=Pyy/N;
subplot(2,1,2);
plot(f(1:N/2),Power(1:N/2),'b');
xlabel('频率(Hz)');
ylabel('功率');
title('功率谱图(频域图');
grid on;
fprintf(' 积分肌电值IEMG = %7.5f \n', mean(abs(y)));
fprintf(' 均方根有效值RMS= %7.5f \n', sqrt(mean(y.^2)) );
问题4：
女12男1综合数据的小波变换后的无偏估计：
clc;clear;
load 'E:\A\Database 1\female_1.mat' %载入仪器测量的female1的原始数据，存储在.mat 文件中
A1=mean(palm_ch1,1); %求cyl 动作在6秒的电位均值
load 'E:\A\Database 1\female_2.mat' %载入仪器测量的female2的原始数据，存储在.mat 文件中
A2=mean(palm_ch1,1); %求cyl 动作在6秒的电位均值
load 'E:\A\Database 1\male_1.mat' %载入仪器测量的female3的原始数据，存储在.mat 文件中
A3=mean(palm_ch1,1); %求cyl 动作在6秒的电位均值
a=[A1;A2;A3]; %将三组均值合在一个数组中
A=mean(a,1); %求数组均值
%进行n 小波变换，确定阈值和去噪
s1=wden(A,'rigrsure','s','one',5,'sym8'); %自适应阈值选择采用无偏估计
%plot(s1);
%求取原始数据N 层小波分析后的无偏估计值
N=3000;
u=sqrt(sum((s1-mean(s1)).^2)/N);
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf(' 平均数值 = %7.3f \n',mean(s1)) %输出本次采样被测量平均值
fprintf(' 标准方差 = %7.3f \n',std(s1)) %输出本次采样数据标准差
fprintf(' n层小波分析的无偏估计值 = %7.3f \n',u) %输出本次采样数无偏估计值
女3男2小波变换后的无偏估计
clc;clear;
load 'E:\A\Database 1\female_3.mat' %载入仪器测量的female1的原始数据，存储在.mat 文件中
A1=mean(lat_ch1,1); %求lat 动作在6秒的电位均值
load 'E:\A\Database 1\male_2.mat' %载入仪器测量的female2的原始数据，存储在.mat 文件中
A2=mean(lat_ch1,1); %求lat 动作在6秒的电位均值
a=[A1;A2]; %将三组均值合在一个数组中
A=mean(a,1); %求数组均值
%进行n 小波变换，确定阈值和去噪
s1=wden(A,'rigrsure','s','one',1,'sym8'); %自适应阈值选择采用无偏估计
%plot(s1);
%求取原始数据N 层小波分析后的无偏估计值
N=3000;
u=sqrt(sum((s1-mean(s1)).^2)/N);
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf(' 平均数值 = %7.3f \n',mean(s1)) %输出本次采样被测量平均值
fprintf(' 标准方差 = %7.3f \n',std(s1)) %输出本次采样数据标准差
fprintf(' n层小波分析的无偏估计值 = %7.3f \n',u) %输出本次采样数无偏估计值
女12男1和女3男2的小波分析后的拟合曲线对比
%运行以下程序:
clc;clear;
A=[0.071 0.043 0.014 0.003 0.002 0.001 0.001 0.001];%小波变换后的无偏估计值
B=[1:8]; %构建时间序列
xdata=B; %将C 中的第一列赋值给x ，形成时间序列 ydata=A; %将C 中的第二列赋值给y ，形成被测量序列 x0=[1 1 1 1]'; %初始点选为全1向量
x=lsqcurvefit('cf',x0,xdata,ydata)
plot(xdata,ydata,'bo');
xi=0.002:0.002:8;
y=cf(x,xi);
grid on;
hold on;
plot(xi,y,'b')
xlabel('x')
ylabel('y')
clear;
A=[0.058 0.034 0.011 0.003 0.002 0.001 0.001 0.001];%小波变换后的无偏估计值
B=[1:8]; %构建时间序列
xdata=B; %将C 中的第一列赋值给x ，形成时间序列 ydata=A; %将C 中的第二列赋值给y ，形成被测量序列 x0=[1 1 1 1]'; %初始点选为全1向量
x=lsqcurvefit('cf',x0,xdata,ydata)
plot(xdata,ydata,'rx');
xi=0.002:0.002:8;
y=cf(x,xi);
grid on;
hold on;
plot(xi,y,'r')
xlabel('x')
ylabel('y')
title('lsqcurvefit函数曲线拟合')
legend('原始数据点',' 原始拟合曲线',' 待识别数据点',' 待识别拟合曲线')
问题5：
day123的cly 的功率谱
%FFT变换, 获得采样数据的能量图和功率图
clear;
close all;
%求功率谱day1
load 'E:\A\Database 2\male_day_1' %通过仪器测量的原始数据
A=cyl_ch1; %将测量数据赋给A ，A 可以为其他数据 A=mean(cyl_ch1,1); %求平均值
B=[0.002:0.002:5]; %构建时间序列
C=[B;A]; %构建时间数值序列
x=C(1,:); %将C 中的第一列赋值给x ，形成时间序列 y=C(2,:); %将C 中的第二列赋值给y ，形成被测量序列
%快速傅立叶变换
y=y-mean(y); %消去直流分量，使频谱更能体现有效信息
Fs=500; %得到原始数据female_1.mat时，仪器的采样频率。其实就是length(x)/(max(x)-min(x));
N=3000; %female_1.mat中的被测量个数，即采样个数。其实就是length(y);
z=fft(y,3000); %做10000点傅里叶变换
%频谱分析
f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N; %幅值，单位同被测变量y
Pyy=Mag.^2; %能量；对实数系列X ，有
X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
Power=Pyy/N;
subplot(3,1,1);
plot(f(1:N/2),Power(1:N/2),'b');
xlabel('频率(Hz)');
ylabel('day1功率');
title('day1功率谱图(频域图');
grid on;
%求功率谱day2
clear;
load 'E:\A\Database 2\male_day_2' %通过仪器测量的原始数据
A=cyl_ch1; %将测量数据赋给A ，A 可以为其他数据 A=mean(cyl_ch1,1); %求平均值
B=[0.002:0.002:5]; %构建时间序列
C=[B;A]; %构建时间数值序列
x=C(1,:); %将C 中的第一列赋值给x ，形成时间序列 y=C(2,:); %将C 中的第二列赋值给y ，形成被测量序列
%快速傅立叶变换
y=y-mean(y); %消去直流分量，使频谱更能体现有效信息
Fs=500; %得到原始数据female_1.mat时，仪器的采样频率。其实就是length(x)/(max(x)-min(x));
N=3000; %female_1.mat中的被测量个数，即采样个数。其实就是length(y);
z=fft(y,3000); %做10000点傅里叶变换
f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N; %幅值，单位同被测变量y
Pyy=Mag.^2;
Power=Pyy/N;
subplot(3,1,2);
plot(f(1:N/2),Power(1:N/2),'r');
xlabel('频率(Hz)');
ylabel('day2功率');
title('day2功率谱图(频域图');
grid on;
%求功率谱day3
clear;
load 'E:\A\Database 2\male_day_3' %通过仪器测量的原始数据
A=cyl_ch1; %将测量数据赋给A ，A 可以为其他数据 A=mean(cyl_ch1,1); %求平均值
B=[0.002:0.002:5]; %构建时间序列
C=[B;A]; %构建时间数值序列
x=C(1,:); %将C 中的第一列赋值给x ，形成时间序列 y=C(2,:); %将C 中的第二列赋值给y ，形成被测量序列
%快速傅立叶变换
y=y-mean(y); %消去直流分量，使频谱更能体现有效信息
Fs=500; %得到原始数据female_1.mat时，仪器的采样频率。其实就是length(x)/(max(x)-min(x));
N=3000; %female_1.mat中的被测量个数，即采样个数。其实就是length(y);
z=fft(y,3000); %做10000点傅里叶变换
%频谱分析
f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N; %幅值，单位同被测变量y
Pyy=Mag.^2; %能量；对实数系列X ，有
X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
Power=Pyy/N;
subplot(3,1,3);
plot(f(1:N/2),Power(1:N/2),'y');
xlabel('频率(Hz)');
ylabel('day3功率');
title('day3功率谱图(频域图');
grid on;
day123的cly 的无偏估计
clc;clear;
load 'E:\A\Database 2\male_day_1' %载入仪器测量的female1的原始数据，存储在.mat 文件中
A=mean(cyl_ch1,1); %求lat 动作在5秒的电位均值
%进行n 小波变换，确定阈值和去噪
s1=wden(A,'rigrsure','s','one',5,'sym8'); %自适应阈值选择采用无偏估计
%plot(s1);
%求取原始数据N 层小波分析后的无偏估计值
N=3000;
u=sqrt(sum((s1-mean(s1)).^2)/N);
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf(' n层小波分析的无偏估计值 = %7.3f \n',u) %输出本次采样数无偏估计值
%day2的无偏估计
clear;
load 'E:\A\Database 2\male_day_2' %载入仪器测量的female1的原始数据，存储在.mat 文件中
A=mean(cyl_ch1,1); %求lat 动作在5秒的电位均值
%进行n 小波变换，确定阈值和去噪
s1=wden(A,'rigrsure','s','one',5,'sym8'); %自适应阈值选择采用无偏估计
%plot(s1);
%求取原始数据N 层小波分析后的无偏估计值
N=3000;
u=sqrt(sum((s1-mean(s1)).^2)/N);
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf(' n层小波分析的无偏估计值 = %7.3f \n',u) %输出本次采样数无偏估计值
%day3的无偏估计
clear;
load 'E:\A\Database 2\male_day_3' %载入仪器测量的female1的原始数据，存储在.mat 文件中
A=mean(cyl_ch1,1); %求lat 动作在5秒的电位均值
%进行n 小波变换，确定阈值和去噪
s1=wden(A,'rigrsure','s','one',5,'sym8'); %自适应阈值选择采用无偏估计
%plot(s1);
%求取原始数据N 层小波分析后的无偏估计值
N=3000;
u=sqrt(sum((s1-mean(s1)).^2)/N);
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf(' n层小波分析的无偏估计值 = %7.3f \n',u) %输出本次采样数无偏估计值
day123的拟合曲线
%运行以下程序:
clc;clear;
A=[0.014 0.009 0.003 0.001 0.001 0.001 0.001 0.001];%小波变换后的无偏估计值 B=[1:8]; %构建时间序列
xdata=B; %将C 中的第一列赋值给x ，形成时间序列 ydata=A; %将C 中的第二列赋值给y ，形成被测量序列 x0=[1 1 1 1]'; %初始点选为全1向量
x=lsqcurvefit('cf',x0,xdata,ydata)
plot(xdata,ydata,'bo');
xi=0.002:0.002:8;
y=cf(x,xi);
grid on;
hold on;
plot(xi,y,'b')
xlabel('x')
ylabel('y')
clear;
A=[0.015 0.011 0.003 0.001 0.001 0.001 0.001 0.001];%小波变换后的无偏估计值 B=[1:8]; %构建时间序列
xdata=B; %将C 中的第一列赋值给x ，形成时间序列 ydata=A; %将C 中的第二列赋值给y ，形成被测量序列 x0=[1 1 1 1]'; %初始点选为全1向量
x=lsqcurvefit('cf',x0,xdata,ydata)
plot(xdata,ydata,'rx');
xi=0.002:0.002:8;
y=cf(x,xi);
grid on;
hold on;
plot(xi,y,'r')
xlabel('x')
ylabel('y')
title('lsqcurvefit函数曲线拟合')
legend('原始数据点',' 原始拟合曲线',' 待识别数据点',' 待识别拟合曲线')
clc;clear;
A=[0.024 0.016 0.007 0.001 0.001 0.001 0.001 0.001];%小波变换后的无偏估计值 B=[1:8]; %构建时间序列
xdata=B; %将C 中的第一列赋值给x ，形成时间序列 ydata=A; %将C 中的第二列赋值给y ，形成被测量序列 x0=[1 1 1 1]'; %初始点选为全1向量
x=lsqcurvefit('cf',x0,xdata,ydata)
plot(xdata,ydata,'gx');
xi=0.002:0.002:8;
y=cf(x,xi);
grid on;
hold on;
plot(xi,y,'g')
xlabel('x')
ylabel('y')
legend('day1原始数据点','day1原始拟合曲线','day2待识别数据点','day2待识别拟合曲线','day3待识别数据点','day2待识别拟合曲线')