一、理论基础

肌电信号(EMG)是一种非常微弱的生物电信号，它与神经肌肉活动密切相关，其中蕴涵着很多与肢体运动相关联的信息。肌电信号可通过表面肌电拾取电极或针式肌电拾取电极加以引导、记录，通过表面电极拾取的肌电信号称之为表面肌电信号，而通过针电极拾取的肌电信号称之为针肌电信号。目前，肌电信号已经深入应用到临床医学、运动医学、生物医学工程等领域。特征提取是肌电信号分析的基础，肌电分析被广泛地应用于肌肉疲劳、重症肌无力、肌强直和肌萎缩等各种肌肉疾病的临床诊断；肌电图也是运动员训练中动作和强度分析的依据；另外，还可利用人体表面肌电的某些特征进行模式分类，进而驱动假肢或其它运动机械(如机器人)的各种动作，实现仿生控制。由此可见，提取有效的肌电信号特征对肌电信号分析意义重大。随着检测技术、信息处理技术及计算机技术的发展，使肌电信号的定量分析及更深入的研究成为可能。

二、案例背景

本项目的目的是开发利用小波分析处理针肌电信号的信号处理方法，这将有助于神经肌肉疾病患者的诊断。

使用小波分析的方法分析肌电图信号，从而实现病人数据的分析。

临床肌电图是评估神经肌肉疾病患者的标准方法。目前没有自动诊断患者的方法，因此诊断必须由临床医生通过目视检查信号或在扬声器上收听其属性来确定。我们已经研究了小波变换作为一种从信号中获得比仅时域频域（FFT）分析更多信息的手段。组合的时频信息为创建患者诊断方法提供了真正的可能性。我们以前发表过几篇关于这个问题的文章。

此段说明，采用传统的快速傅立叶变换方法已经无法达到所要的目的。

有许多不同的小波可用。该项目将首先使用B-sline小波，因为它在时间局部化方面具有理想的特性，但许多包括其他小波。到目前为止，我们一直在使用复杂的Labview程序在PC上收集肌电图数据。新的EMG机器将允许在机器本身上收集数据，但我们需要一个带有医生界面的Labview程序来控制信号向患者的呈现，以便在最大自主收缩（MVC）和大约50%的MVC时收集信号。我们希望对两种信号的比较分析将揭示患者和正常受试者之间的差异。

这段的要求，就是使用B样条小波分析/或其余小波分析的方法来分析肌电图数据；

该项目将从为医生设置控制程序开始。一旦在阿伯丁收集到数据，这些数据将被匿名并发送到邓迪，在那里将使用小波分析方法进行分析。学生将开发分析方法。分析程序将用Matlab编写。该项目将适合具有编程经验和信号处理知识的人员。小波理论需要良好的数学技能才能理解使用MATLAB编写相关代码来实现分析。

三、matlab程序

clc;clear;close all;warning off;[signals,Fs]=wavread('data\triceps1.wav');%由于数据采样率太高，对于数据较大的数据，只选择其中一段if length(signals) > 1000000   signal = signals(1:1000000);else   signal = signals;endclear signals;figure;plot(signal);title('Initial data');save mat\signal.mat signal save mat\Fs.mat Fs clear signal%Step1:利用小波变换做预处理，去噪%Step1:利用小波变换做预处理，去噪%Step1:利用小波变换做预处理，去噪%画出原始信号 load mat\signal.mats = signal';figure;subplot(221); plot(s); title('原始信号');grid; %用db1小波对原始信号进行3层分解并提取系数 [c,l]=wavedec(s,3,'db1'); ca3=appcoef(c,l,'db1',3); cd3=detcoef(c,l,3); cd2=detcoef(c,l,2); cd1=detcoef(c,l,1); %对信号进行强制性消噪处理并图示结果 cdd3=zeros(1,length(cd3)); cdd2=zeros(1,length(cd2)); cdd1=zeros(1,length(cd1)); c1=[ca3 cdd3 cdd2 cdd1]; s1=waverec(c1,l,'db1'); subplot(2,2,2); plot(s1); title('强制消噪后的信号');grid; %用ddencmp函数获得信号的默认阈值,使用wdencmp命令函数实现消噪过程 [thr,sorh,keepapp]=ddencmp('den','wv',s); s2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp); subplot(2,2,3); plot(s2);title('默认阈值消噪后的信号');grid; %用给定的软阈值进行消噪处理 cd1soft=wthresh(cd1,'s',1.465); cd2soft=wthresh(cd2,'s',1.823); cd3soft=wthresh(cd3,'s',2.768); c2=[ca3 cd3soft cd2soft cd1soft]; s3=waverec(c2,l,'db1'); subplot(2,2,4); plot(s3); title('给定软阈值消噪后的信号');grid set(gcf,'color','w');save mat\delete.mat s s1 s2 s3clear all%Step2:直接对信号进行功率谱的计算分析%Step2:直接对信号进行功率谱的计算分析%Step2:直接对信号进行功率谱的计算分析load mat\delete.matload mat\Fs.matsignals = s2;figure;plotspec(signals,1/Fs)set(gcf,'color','w');save mat\signals.mat signalsclear all%Step3:对该信号进行小波变换%Step3:对该信号进行小波变换%Step3:对该信号进行小波变换%db10小波进行4层分解%一维小波分解load mat\delete.matload mat\Fs.matload mat\signals.mattype    = 'db5';n       = 7;signals = s2;[c,l]   = wavedec(signals,n,type);%分解a7 = wrcoef('a',c,l,type,7);a6 = wrcoef('a',c,l,type,6);a5 = wrcoef('a',c,l,type,5);a4 = wrcoef('a',c,l,type,4);a3 = wrcoef('a',c,l,type,3);a2 = wrcoef('a',c,l,type,2);a1 = wrcoef('a',c,l,type,1);%重构第1～4层细节系数d7 = wrcoef('d',c,l,type,7);d6 = wrcoef('d',c,l,type,6);d5 = wrcoef('d',c,l,type,5);d4 = wrcoef('d',c,l,type,4);d3 = wrcoef('d',c,l,type,3);d2 = wrcoef('d',c,l,type,2);d1 = wrcoef('d',c,l,type,1);%显示细节信号figuresubplot(5,3,1);plot(d7);title('d7');subplot(5,3,2);plot(d6);title('d6');subplot(5,3,3);plot(d5);title('d5');subplot(5,3,4);plot(d4);title('d4');subplot(5,3,5);plot(d3);title('d3');subplot(5,3,6);plot(d2);title('d2');subplot(5,3,7);plot(d1);title('d1');subplot(5,3,9) ;plot(a7);title('a7');subplot(5,3,10);plot(a6);title('a6');subplot(5,3,11);plot(a5);title('a5');subplot(5,3,12);plot(a4);title('a4');subplot(5,3,13);plot(a3);title('a3');subplot(5,3,14);plot(a2);title('a2');subplot(5,3,15);plot(a1);title('a1');set(gcf,'color','w');save mat\d.mat d1 d2 d3 d4 d5 d6 d7clear allload mat\delete.matload mat\signals.matload mat\d.matload mat\Fs.mattype    = 'db5';n       = 7;%7层小波包分解T=wpdec(s2,n,type);%重构低频信号y1=wprcoef(T,[7,1]);figure;plotspec(y1,1/Fs)set(gcf,'color','w');%第1层的频谱图figure;plotspec(d1,1/Fs);title('d1');%第2层的频谱图figure;plotspec(d2,1/Fs);title('d2');set(gcf,'color','w');%第3层的频谱图figure;plotspec(d3,1/Fs);title('d3');set(gcf,'color','w');%第4层的频谱图figure;plotspec(d4,1/Fs);title('d4');set(gcf,'color','w');%第5层的频谱图figure;plotspec(d5,1/Fs);title('d5');set(gcf,'color','w');%第6层的频谱图figure;plotspec(d6,1/Fs);title('d6');set(gcf,'color','w');%第7层的频谱图figure;plotspec(d7,1/Fs);title('d7');set(gcf,'color','w');clear alltype    = 'db5';n       = 7;load mat\delete.mat%计算各层能量分布for i=1:n    wpt=wpdec(s2,i,type);       e=wenergy(wpt);    E=zeros(1,length(e));    for j=1:2^i        E(j)=sum(abs(wprcoef(wpt,[i,j-1])).^2);    end    figure(20);    subplot(7,1,i);    bar(e);    axis([0 length(e) 0 130]);    title(['第 ',num2str(i), ' 层']);endclear all;

clc;clear;close all;warning off;[signals,Fs]=wavread('data\EDC3.wav');%由于数据采样率太高，对于数据较大的数据，只选择其中一段if length(signals) > 1000000   signal = signals(1:1000000);else   signal = signals;endclear signals;x=signal;l0=length(x);%计算B样条小波函数的尺度特性曲线m     = 3;%B样条阶数N     = func_bspline(2*m,2*m+1);%求样条函数在整数点的值w     = 0:0.1*pi:10*pi;%求尺度函数mea_w = fun_bspline_wavelet(m,w,N);figure;subplot(131);plot(w/pi,abs(mea_w));title('尺度函数的幅频图');mea_ang=atan(imag(mea_w)./(real(mea_w)+eps));%计算相位角subplot(132),plot(w/pi,mea_ang);title('尺度函数的相频图');%求反傅立叶变换，求出尺度函数的时域函数。mea_x=ifft(mea_w);                                  mea_x=ifftshift(mea_x);                              subplot(1,3,3);plot(w/pi,real(mea_x));title('尺度函数的时域图');set(gcf,'color','w');x=x';%三次样条小波的滤波器系数ld=[0,0.0625,0.25,0.375,0.25,0.0625,0,0]; hd=[-0.00008,-0.01643,-0.10872,-0.59261,0.59261,0.10872,0.01643,0.00008];   %第1层小波分解ld=rot90(ld,2);hd=rot90(hd,2);a1=wconv1(x,ld);d1=wconv1(x,hd);l1=length(a1);a1=a1([l1-l0+1:l1]);d1=d1([l1-l0+1:l1]);a1=dyaddown(a1);d1=dyaddown(d1);l1=length(a1);%第2层小波分解a2=wconv1(a1,ld);d2=wconv1(d1,hd);l2=length(a2);a2=a2([l2-l1+1:l2]);d2=d2([l2-l1+1:l2]);a2=dyaddown(a2);d2=dyaddown(d2);l2=length(a2);%第3层小波分解a3=wconv1(a2,ld);d3=wconv1(d2,hd);l3=length(a3);a3=a3([l3-l2+1:l3]);d3=d3([l3-l2+1:l3]);a3=dyaddown(a3);d3=dyaddown(d3);l3=length(a3);%第4层小波分解a4=wconv1(a3,ld);d4=wconv1(d3,hd);l4=length(a4);a4=a4([l4-l3+1:l4]);d4=d4([l4-l3+1:l4]);a4=dyaddown(a4);d4=dyaddown(d4);l4=length(a4);%显示细节信号figuresubplot(421);plot(a4);title('a4');subplot(422);plot(a3);title('a3');subplot(423);plot(a2);title('a2');subplot(424);plot(a1);title('a1');subplot(425);plot(d4);title('d4');subplot(426);plot(d3);title('d3');subplot(427);plot(d2);title('d2');subplot(428);plot(d1);title('d1');set(gcf,'color','w');%第1层的频谱图figure;plotspec(d1,1/Fs);title('d1');set(gcf,'color','w');%第2层的频谱图figure;plotspec(d2,1/Fs);title('d2');set(gcf,'color','w');%第3层的频谱图figure;plotspec(d3,1/Fs);title('d3');set(gcf,'color','w');%第4层的频谱图figure;plotspec(d4,1/Fs);title('d4');set(gcf,'color','w');