肌电信号的时域和频域分析

肌电信号的时域和频域分析摘要：肌电信号是产生肌肉力的电信号根源，它是肌肉中很多运动单元动作电位在时间和空间上的叠加，反映了神经，肌肉的功能状态，在基础医学研究、临床诊断和康复工程中有广泛的应用。其种类重要有两种：一，临床肌电图检查多采用针电极插入肌肉检测肌电图，其优点是干扰小，定位性好，易识别，但由于它是一种有创伤的检测方法，其应用收到了一定的限制。二，表面肌电则是从人体皮肤表面通过电极记录下来的神经肌肉活动时发放的生物电信号，属于无创伤性，操作简单，病人易接受，有着广泛的应用前景。本次设计基于matlab用小波变换对肌电信号进行消噪处理，分别选用20N的肌电信号数据和50N的肌电数据进行对比，最后在GUI界面上完成相应的功能处理。关键字：肌电信号Matlab小波去噪GUI第一章绪论肌电信号是产生肌肉力的电信号根源，它是肌肉中很多运动单元动作电位在时间和空间上的叠加，反映了神经，肌肉的功能状态，在基础医学研究、临床诊断和康复工程中有广泛的应用。其种类重要有两种：一，临床肌电图检查多采用针电极插入肌肉检测肌电图，其优点是干扰小，定位性好，易识别，但由于它是一种有创伤的检测方法，其应用收到了一定的限制。二，表面肌电则是从人体皮肤表面通过电极记录下来的神经肌肉活动时发放的生物电信号，属于无创伤性，操作简单，病人易接受，有着广泛的应用前景。肌电信号本身是一种较微弱的电信号。检测和记录表面肌电信号，需要考虑的主要问题是尽量消除噪声和干扰的影响,提高信号的保真度[1]。2第二章肌电信号的时域分析2.1肌电信号时域图的显示及比较肌电信号采用两个不同的数据进行比较，通过比较时域图及其特性来进行分析[2]。其图像如下所示：如上图所示：肌电数据分别是同一个体在20N的力和50N的力所反映的图像。可以看出在不同作用力时，其图像的差别很大。2.2时域参数2.2.1均值对于一个随机变量来说，均值是一个很重要的数值特征。粗略的说，就是来描述一个群体的平均水平。其严格的数学定义非常的简单，就是一个随机变量关于概率测度的积分。这样的积分在测度轮或者实分析里是没有什么直观的解释的。而在概率论里却成为了一个群体的主要指标。在此处，均值表示肌电信号的平均水平。2.2.2标准差标准差（StandardDeviation），也称均方差（meansquareerror），是各数据偏离平均数的距离的平均数，它是离均差平方和平均后的方根，用σ表示。3标准差是方差的算术平方根。标准差能反映一个数据集的离散程度。平均数相同的，标准差未必相同。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)));4第三章肌电信号的时域分析3.1傅里叶变换傅里叶是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后就很容易看出特征了。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。本次设计选取N=20000来进行肌电信号的频域分析[3]。通过傅里叶变换后，两个不同信号的幅频图如下所示：肌电相频图如下所示：53.2功率谱分析能量谱密度、功率谱密度函数表示信号的能量、功率密度随频率变化的情况。通过研究功率谱密度，可以帮助了解信号的功率分布情况，确定信号的频带等。功率密度谱虽然描述了随机信号的功率在各个不同频率上的分布，但因为它仅与幅度频谱有关，没有相位信息，所以从已知功率谱还难以完整地恢复原来的功率信号。通过执行相应程序后，其功率谱的显示图如下所示：其部分程序代码如下所示：globala;globalt;6globals;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)7第四章小波去噪分析小波分析方法是一种窗口大小（即窗口面积）固定但其形状可改变，时间窗和频率都可改变的时频局域化分析方法，即在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，所以被称为数学显微镜。正是这种特性，是小波变换具有对信号的自适应性[4]。4.1常用的小波分析方法4.1.1小波消噪一维的小波消噪过程主要分为以下三个步骤：①一维信号的分解：获得各尺度上的细节分量和近似分量；②细节分量的阈值化：对1到N尺度的上的每层细节分量选取合适的阈值，进行阈值化处理，保留有用信息，去除噪声成分；③一维信号的重构：根据小波分解的第N层近似分量和经阈值化处理后的第1层到第N层细节分量，进行一维信号的小波重构。本文还利用小波的多分辨率特性将心电信号分解为多层，然后进行有用信号的重构，将一些噪声信号去除。4.1.2小波函数的选取同一信号，选取不同的小波函数进行处理，将得到不同的效果，所以小波函数的选取显得尤为重要。对于心电信号滤波来说，选择支集长度较短的小波可提高处理的实时性；选取对称性的小波可满足相移为基本线性，使心电信号不失真；选取正则性的小波可使重构以后的信号比较平滑。经过对比降噪效果，在本文中用db5小波基进行消噪处理。4.1.2阈值选取阈值的确定是小波收缩消噪最关键的一步，阈值过小，则方差偏大，数据欠平滑；阈值过大，会使数据过平滑，信号的奇异性可能丧失。对小波系数进行阈值操作过程中，有两种方式，其一对每一个小波系数进行阈值操作，其二是成块习俗进行阈值操作。由信号的奇异性理论，心电信号里的噪声具有负的奇异性，8其幅度和稠密度随尺度的增大而减小，而信号则相反，因此阈值的选取不能单一。阙值函数体现了对小波系数的不同处理策略，主要分为软阈值函数和硬阈值函数它们的基本思想都是去除小幅值的系数，对大幅值的系数进行保留。它们的去噪效果各有特点(可以从后面的具体例子中看出)：硬阈值法可以很好的保留信号的边缘等局部特征，但去噪的结果具有较大的方差，会出现伪Gibbs现象等视觉失真；软闭值法去噪结果相对平滑，但有较大的偏差，可能会造成边缘模糊等失真现象。4.2小波去噪的具体分析为了确定和比较本次设计中小波去噪的效果，先对小波进行一层分解，二层分解，三层分解，四层分解，让它们单独对同一信号进行消噪处理。其结果图如下所示：9通过比较上图可以看出，第四幅图像的去噪效果比较明显，因次在本次设计中选择小波四层分解对信号进行去噪处理。分别对20N的肌电信号和50N的肌电信号进行小波去噪，去噪前后的效果图如下所示：10左边是20N小波去噪前后的图，右边是50N作用前后的效果图。通过前后图像的比较可以看出，当选用四层小波分解进行滤波时，已经能去掉大部分噪声的干扰，保存有用的信号。小波去噪部分程序代码如下:globala;globalt;globals;globaly;%globalIR;%globalSORH;M=length(a);11N=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:Lm(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);12MSE=M/N;SM=SNR/MSE;%K=length(d);%t1=(0:K-1)/2000;plot(handles.axes5,XC(1:20000));小波去噪前后的幅频和相频图如下：13表面肌电信号一般只有微伏级电压,信号中往往夹带着低频(接近直流)和高频的干扰信号,真正有用的肌电信号大致在10Hz~500Hz之间。由上面的对比图可知，消除了噪声信号，有用信号得以保存。其部分程序代码如下所示：globala;globalt;%globalh;[CL]=wavedec(a,3,'db5');cA3=appcoef(C,L,'db5',3);cD1=detcoef(C,L,1);14cD2=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);去噪前后的功率谱图如下所示：15部分程序代码如下：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');16thr2=thselect(cD2,'rigrsure');thr3=thselect(cD