肌电信号的时域和频域分析.docx
- 文档编号:24594203
- 上传时间:2023-05-29
- 格式:DOCX
- 页数:33
- 大小:171.54KB
肌电信号的时域和频域分析.docx
《肌电信号的时域和频域分析.docx》由会员分享,可在线阅读,更多相关《肌电信号的时域和频域分析.docx(33页珍藏版)》请在冰豆网上搜索。
肌电信号的时域和频域分析
肌电信号的时域和频域分析
摘要:
肌电信号是产生肌肉力的电信号根源,它是肌肉中很多运动单元动作电
位在时间和空间上的叠加,反映了神经,肌肉的功能状态,在基础医学研究、
临床诊断和康复工程中有广泛的应用。
其种类重要有两种:
一,临床肌电图检查多采用针电极插入肌肉检测肌
电图,其优点是干扰小,定位性好,易识别,但由于它是一种有创伤的检测
方法,其应用收到了一定的限制。
二,表面肌电则是从人体皮肤表面通过电
极记录下来的神经肌肉活动时发放的生物电信号,属于无创伤性,操作简单,
病人易接受,有着广泛的应用前景。
本次设计基于matlab用小波变换对肌电信号进行消噪处理,分别选用20N的肌电信号数据和50N的肌电数据进行对比,最后在GUI界面上完成相应的功能处理。
关键字:
肌电信号Matlab小波去噪GUI
第1章绪论
肌电信号是产生肌肉力的电信号根源,它是肌肉中很多运动单元动作电位在时间和空间上的叠加,反映了神经,肌肉的功能状态,在基础医学研究、临床诊断和康复工程中有广泛的应用。
其种类重要有两种:
一,临床肌电图检查多采用针电极插入肌肉检测肌电图,其优点是干扰小,定位性好,易识别,但由于它是一种有创伤的检测方法,其应用收到了一定的限制。
二,表面肌电则是从人体皮肤表面通过电极记录下来的神经肌肉活动时发放的生物电信号,属于无创伤性,操作简单,病人易接受,有着广泛的应用前景。
肌电信号本身是一种较微弱的电信号。
检测和记录表面肌电信号,需要考虑的主要问题是尽量消除噪声和干扰的影响,提高信号的保真度[1] 。
第2章肌电信号的时域分析
2.1肌电信号时域图的显示及比较
肌电信号采用两个不同的数据进行比较,通过比较时域图及其特性来进行分析[2]。
其图像如下所示:
如上图所示:
肌电数据分别是同一个体在20N的力和50N的力所反映的图像。
可以看出在不同作用力时,其图像的差别很大。
2.2时域参数
2.2.1均值
对于一个随机变量来说,均值是一个很重要的数值特征。
粗略的说,就是来描述一个群体的平均水平。
其严格的数学定义非常的简单,就是一个随机变量关于概率测度的积分。
这样的积分在测度轮或者实分析里是没有什么直观的解释的。
而在概率论里却成为了一个群体的主要指标。
在此处,均值表示肌电信号的平均水平。
2.2.2标准差
标准差(StandardDeviation),也称均方差(meansquareerror),是各数据偏离平均数的距离的平均数,它是离均差平方和平均后的方根,用σ表示。
标准差是方差的算术平方根。
标准差能反映一个数据集的离散程度。
平均数相同的,标准差未必相同。
2.2.3方差
方差是各个数据与平均数之差的平方的平均数。
在概率论和数理统计中,方(英文Variance)用来度量随机变量和其数学期望(即均值)之间的偏离程度。
在许多实际问题中,研究随机变量和均值之间的偏离程度有着很重要的意义。
选取一个信号,其执行结果如下所示:
其部分程序代码如下所示:
clear;
clc;
s=load('E:
\肌电信号数据\EMG\EMG3.txt');
%s=load('E:
\肌电信号数据\EMG\bs1.txt');
a=s(:
7);
t=s(:
1);
fprintf('\n数据基本信息:
\n');
fprintf(' 均值=%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)));
第3章肌电信号的时域分析
3.1傅里叶变换
傅里叶是离散傅立叶变换的快速算法,可以将一个信号变换到频域。
有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后就很容易看出特征了。
另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。
一个模拟信号,经过ADC采样之后,就变成了数字信号。
采样得到的数字信号,就可以做FFT变换了。
N个采样点,经过FFT之后,就可以得到N个点的FFT结果。
本次设计选取N=20000来进行肌电信号的频域分析[3]。
通过傅里叶变换后,两个不同信号的幅频图如下所示:
肌电相频图如下所示:
3.2功率谱分析
能量谱密度、功率谱密度函数表示信号的能量、功率密度随频率变化的情况。
通过研究功率谱密度,可以帮助了解信号的功率分布情况,确定信号的频带等。
功率密度谱虽然描述了随机信号的功率在各个不同频率上的分布,但因为它仅与幅度频谱有关,没有相位信息,所以从已知功率谱还难以完整地恢复原来的功率信号。
通过执行相应程序后,其功率谱的显示图如下所示:
其部分程序代码如下所示:
globala;
globalt;
globals;
s=fft(a,2000);
[CL]=wavedec(a,3,'db5');
cA3=appcoef(C,L,'db5',3);
cD1=detcoef(C,L,1);
cD2=detcoef(C,L,2);
cD3=detcoef(C,L,3);
thr1=thselect(cD1,'rigrsure');
thr2=thselect(cD2,'rigrsure');
thr3=thselect(cD3,'rigrsure');
TR=[thr1,thr2,thr3];
SORH='s';
[XC,CXC,LXC,PERFO,PERF2]=wdencmp('lvd',a,...
'db5',3,TR,SORH);
y1=fft(XC,20000);
fs=2000;
N=length(y1);
mag=abs(s);
f=(0:
N-1)/N*fs;
power1=(mag.^2)/2000;plot(handles.axes1,power1)
第4章小波去噪分析
小波分析方法是一种窗口大小(即窗口面积)固定但其形状可改变,时间窗和频率都可改变的时频局域化分析方法,即在低频部分具有较高的频率分辨率和较低的时间分辨率,在高频部分具有较高的时间分辨率和较低的频率分辨率,所以被称为数学显微镜。
正是这种特性,是小波变换具有对信号的自适应性[4]。
4.1常用的小波分析方法
4.1.1小波消噪
一维的小波消噪过程主要分为以下三个步骤:
①一维信号的分解:
获得各尺度上的细节分量和近似分量;②细节分量的阈值化:
对1到N尺度的上的每层细节分量选取合适的阈值,进行阈值化处理,保留有用信息,去除噪声成分;③一维信号的重构:
根据小波分解的第N层近似分量和经阈值化处理后的第1层到第N层细节分量,进行一维信号的小波重构。
本文还利用小波的多分辨率特性将心电信号分解为多层,然后进行有用信号的重构,将一些噪声信号去除。
4.1.2小波函数的选取
同一信号,选取不同的小波函数进行处理,将得到不同的效果,所以小波函数的选取显得尤为重要。
对于心电信号滤波来说,选择支集长度较短的小波可提高处理的实时性;选取对称性的小波可满足相移为基本线性,使心电信号不失真;选取正则性的小波可使重构以后的信号比较平滑。
经过对比降噪效果,在本文中用db5小波基进行消噪处理。
4.1.2阈值选取
阈值的确定是小波收缩消噪最关键的一步,阈值过小,则方差偏大,数据欠平滑;阈值过大,会使数据过平滑,信号的奇异性可能丧失。
对小波系数进行阈值操作过程中,有两种方式,其一对每一个小波系数进行阈值操作,其二是成块习俗进行阈值操作。
由信号的奇异性理论,心电信号里的噪声具有负的奇异性,其幅度和稠密度随尺度的增大而减小,而信号则相反,因此阈值的选取不能单一。
阙值函数体现了对小波系数的不同处理策略,主要分为软阈值函数和硬阈值函数它们的基本思想都是去除小幅值的系数,对大幅值的系数进行保留。
它们的去噪效果各有特点(可以从后面的具体例子中看出):
硬阈值法可以很好的保留信号的边缘等局部特征,但去噪的结果具有较大的方差,会出现伪Gibbs现象等视觉失真;软闭值法去噪结果相对平滑,但有较大的偏差,可能会造成边缘模糊等失真现象。
4.2小波去噪的具体分析
为了确定和比较本次设计中小波去噪的效果,先对小波进行一层分解,二层分解,三层分解,四层分解,让它们单独对同一信号进行消噪处理。
其结果图如下所示:
通过比较上图可以看出,第四幅图像的去噪效果比较明显,因次在本次设计中选择小波四层分解对信号进行去噪处理。
分别对20N的肌电信号和50N的肌电信号进行小波去噪,去噪前后的效果图如下所示:
左边是20N小波去噪前后的图,右边是50N作用前后的效果图。
通过前后图像的比较可以看出,当选用四层小波分解进行滤波时,已经能去掉大部分噪声的干扰,保存有用的信号。
小波去噪部分程序代码如下:
globala;
globalt;
globals;
globaly;
%globalIR;
%globalSORH;
M=length(a);
N=length(y);
p=size(a);
s=a(1:
20000);
[CL]=wavedec(a,4,'db5');
cA3=appcoef(C,L,'db5',4);
cD1=detcoef(C,L,1);
cD2=detcoef(C,L,2);
cD3=detcoef(C,L,3);
cD4=detcoef(C,L,4);
thr1=thselect(cD1,'rigrsure');
thr2=thselect(cD2,'rigrsure');
thr3=thselect(cD3,'rigrsure');
thr4=thselect(cD4,'rigrsure');
TR=[thr1,thr2,thr3,thr4];
SORH='s';
[XC,CXC,LXC,PERFO,PERF2]=wdencmp('lvd',a,...
'db5',4,TR,SORH);
L=p
(2);
x=a;
h=XC;
F=0;
M=0;
forii=1:
L
m(ii)=(x(ii)-y(ii))^2;
t(ii)=y(ii)^2;
f(ii)=t(ii)/m(ii);
F=F+f(ii);
M=M+m(ii);
end;
SNR=10*log10(F);
MSE=M/N;
SM=SNR/MSE;
%K=length(d);
%t1=(0:
K-1)/2000;
plot(handles.axes5,XC(1:
20000));
小波去噪前后的幅频和相频图如下:
表面肌电信号一般只有微伏级电压,信号中往往夹带着低频(接近直流)和高频的干扰信号,真正有用的肌电信号大致在10Hz~500Hz之间。
由上面的对比图可知,消除了噪声信号,有用信号得以保存。
其部分程序代码如下所示:
globala;
globalt;
%globalh;
[CL]=wavedec(a,3,'db5');
cA3=appcoef(C,L,'db5',3);
cD1=detcoef(C,L,1);
cD2=detcoef(C,L,2);
cD3=detcoef(C,L,3);
thr1=thselect(cD1,'rigrsure');
thr2=thselect(cD2,'rigrsure');
thr3=thselect(cD3,'rigrsure');
TR=[thr1,thr2,thr3];
SORH='s';
[XC,CXC,LXC,PERFO,PERF2]=wdencmp('lvd',a,...
'db5',3,TR,SORH);
y1=fft(XC,20000);
fs=2000;
N=length(y1);
mag1=abs(y1);
f=(0:
N-1)/N*fs;
plot(handles.axes5,f,mag1);
去噪前后的功率谱图如下所示:
部分程序代码如下:
globala;
globalt;
globals;
s=fft(a,2000);
[CL]=wavedec(a,3,'db5');
cA3=appcoef(C,L,'db5',3);
cD1=detcoef(C,L,1);
cD2=detcoef(C,L,2);
cD3=detcoef(C,L,3);
thr1=thselect(cD1,'rigrsure');
thr2=thselect(cD2,'rigrsure');
thr3=thselect(cD3,'rigrsure');
TR=[thr1,thr2,thr3];
SORH='s';
[XC,CXC,LXC,PERFO,PERF2]=wdencmp('lvd',a,...
'db5',3,TR,SORH);
y1=fft(XC,2000);
fs=2000;
N=length(y1);
mag=abs(y1);
f=(0:
N-1)/N*fs;
power2=(mag.^2)/2000;
plot(handles.axes5,power2);
第5章GUI界面的设计
GUI是由窗口、光标、按键、菜单、文字说明等对象(Objects)构成的一个用户界面。
用户通过一定的方法(如鼠标或键盘)选择、激活这些图形对象,使计算机产生某种动作或变化,比如实现计算、绘图等。
一个好的GUI能够使程序更加容易的使用。
它提供用户一个常见的界面,还提供一些控件,例如,按钮,列表框,滑块,菜单等。
用户图形界面应当是易理解且操作是可以预告的,所以当用户进行某一项操作,它知道如何去做。
例如,当鼠标在一个按钮上发生了单击事件,用户图形界面初始化它的操作,并在按钮的标签上对这个操作进行描述。
本次设计的GUI界面整体图如下所示:
程序运行之后产生的界面图如下所示:
本次设计模块涵盖多个功能,都是通过按钮来进行相应的信号处理。
只要M文件在当前目录或在Matlab搜索路径上,在Matlab命令窗口输入对应的M文件,就能打如上图所示的图形用户界面。
在此界面上可以进行相关的操作。
例如点击open按扭就会打开如下界面:
Open按钮的回调函数如下填写:
[filename,filepath]=uigetfile('*.txt','选择文件');%选择数据文件
str=[filepathfilename];
其他的按钮都按照此方法去激活,以实现相应的功能。
学习心得
通过这次设计,我对MATLAB有了更深入的理解,学会了在GUI界面通过按钮来实现相应的功能,对以前未接触的小波变换也有了大概了解。
在设计的过程中,我也认识到了自己所学知识的不足。
这也让我再次认识到知识是无尽的,只有不断的充实自己、完善自己的知识理论体系,才能够更好的胜任自己以后的工作。
设计过程中知识的不足也让我更加坚定了终身学习的决心。
参考文献
1.谢平、王娜、林洪斌等主编,信号处理原理及应用。
北京:
机械工业出版社,
2008.10
2.聂祥飞、王海宝、谭泽富主编,Matlab程序设计及其在信号处理中的应用。
成都:
西南交通大学出版社,2005
3.吴大正、高西全等主编,Matlab及在电子信息课程中的应用。
北京:
电子工业出版社,2006.3
4.李培芳、孙晖、李江主编,信号与系统分析基础。
北京:
清华大学出版社,2006.12
程序附录
functionvarargout=a1fig(varargin)
gui_Singleton=1;
gui_State=struct('gui_Name', mfilename,...
'gui_Singleton', gui_Singleton,...
'gui_OpeningFcn',@a1fig_OpeningFcn,...
'gui_OutputFcn', @a1fig_OutputFcn,...
'gui_LayoutFcn', [],...
'gui_Callback', []);
ifnargin&&ischar(varargin{1})
gui_State.gui_Callback=str2func(varargin{1});
end
ifnargout
[varargout{1:
nargout}]=gui_mainfcn(gui_State,varargin{:
});
else
gui_mainfcn(gui_State,varargin{:
});
end
%Endinitializationcode-DONOTEDIT
%---Executesjustbeforea1figismadevisible.
functiona1fig_OpeningFcn(hObject,eventdata,handles,varargin)
%Thisfunctionhasnooutputargs,seeOutputFcn.
%hObject handletofigure
%eventdata reserved-tobedefinedinafutureversionofMATLAB
%handles structurewithhandlesanduserdata(seeGUIDATA)
%varargin commandlineargumentstoa1fig(seeVARARGIN)
%Choosedefaultcommandlineoutputfora1fig
handles.output=hObject;
%Updatehandlesstructure
guidata(hObject,handles);
%UIWAITmakesa1figwaitforuserresponse(seeUIRESUME)
%uiwait(handles.figure1);
%---Outputsfromthisfunctionarereturnedtothecommandline.
functionvarargout=a1fig_OutputFcn(hObject,eventdata,handles)
%varargout cellarrayforreturningoutputargs(seeVARARGOUT);
%hObject handletofigure
%eventdata reserved-tobedefinedinafutureversionofMATLAB
%handles structurewithhandlesanduserdata(seeGUIDATA)
%Getdefaultcommandlineoutputfromhandlesstructure
varargout{1}=handles.output;
%---Executesonbuttonpressinpushbutton1.
functionpushbutton1_Callback(hObject,eventdata,handles)
%hObject handletopushbutton1(seeGCBO)
%eventdata reserved-tobedefinedinafutureversionofMATLAB
%handles structurewithhandlesanduserdata(seeGUIDATA
globals;
globala;
globalt;
[filename,filepath]=uigetfile('*.txt','选择文件');%选择数据文件
str=[filepathfilename];
s=load(str);
a=s(:
7);
t=s(:
1);
%---Executesonbuttonpressinpushbutton2.
functionpushbutton2_Callback(hObject,eventdata,handles)
%hObject handletopushbutton2(seeGCBO)
%eventdata reserved-tobedefinedinafutureversionofMATLAB
%handles structurewithhandlesanduserdata(seeGUIDATA)
try
delete(allchild(handles.axes1));
end
%---Executesonbuttonpressinpushbutton3.
functionpushbutton3_Callback(hObject,eventdata,handles)
%hObject handletopushbutton3(seeGCBO)
%eventdata reserved-tobedefinedinafutureversionofMATLAB
%handles structurewithhandlesanduserdata(seeGUIDATA)
globala;
globalt;
plot(handles.axes1,t,a);title('信号1');
%---Executesonbuttonpressinpushbutton4.
functionpushbutton4_Callback(hObject,eventdata,handles)
%hObject handletopushbutton4(seeGCBO)
%eventdata reserved-tobedefinedinafutureversionofMATLAB
%handles structurewithhandlesanduserdata(seeGUIDATA)
globala;
globalt;
globaly;
%a=EMG3(1:
points,2);
y=fft(a,20000);%对信号进行快速Fourier变换
fs=2000;
N=length(y);
mag=abs(y);%求得Fourier变换后的振幅
f=(0:
N-1)/N*fs;%频率序列
plot(handles.axes1,f,mag);%绘制图形
%---Executesonbuttonpressinpushbutton5.
functionpushbutton5_Callback(hObject,eve
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 电信号 时域 分析