基于B样条小波变换的EMG肌电信号去噪算法matlab仿真

目录

1.B样条小波变换

2.EMG信号

3.MATLAB程序

4.仿真结果

---

        肌电图（Electromyography, EMG）信号是研究神经肌肉系统功能的重要生物电信号，但由于采集过程中存在的各种噪声，如肌腱振动、电极接触噪声、电源干扰等，严重影响了信号的分析和应用。基于B样条小波变换（B-Spline Wavelet Transform, BSWT）的EMG信号去噪算法是一种结合了B样条函数平滑特性和小波变换多尺度分析优势的方法，能够有效地从复杂信号中提取肌电信号的有用成分，同时抑制噪声。

1.B样条小波变换

      B样条（Basis spline）是一种在计算机图形学和数值分析中广泛应用的平滑曲线，其特性在于局部支撑、连续性好且易于控制形状。B样条基函数定义在闭区间[0,1]上，通过递归定义得到，具有良好的局部性和全局光滑性。在小波分析中，B样条函数被用作小波基，生成的小波基具有更好的正交性和紧凑支持性。

       传统的离散小波变换（Discrete Wavelet Transform, DWT）通常基于Haar小波或Daubechies小波等，而B样条小波变换则是基于B样条函数构造的小波基进行的。B样条小波基的生成依赖于B样条函数的升阶和降阶操作，从而形成一个多分辨率分析框架。

       利用多尺度分析理论构造出的B样条小波函数同时具备紧支集，对称性和反对称性，很高的正则性，高阶消失矩及半正交性等优良性质。因此，B样条小波的应用范围也非常广泛。m阶B样条小波：

2.EMG信号

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

      2.噪声抑制：在不同尺度上，噪声和信号的分布特性不同。通常，高频尺度的小波系数主要包含噪声信息，而低频尺度的尺度系数则主要携带信号的有用信息。因此，可以对高频系数施加硬阈值或软阈值处理来去除噪声

3.MATLAB程序

%计算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);


%显示细节信号
figure
subplot(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');
up4094

4.仿真结果

·输入的信号：

   ·通过小波变换的方法去噪声，其仿真结果如下所示：

      从图可以看出，经全局阈值和分层阈值方法降噪后的信号都很好的保留了信号发展初期的高频特性。分层阈值法虽然损失了部分性能（与原信号的相似性：经计算，分层阈值降噪信号与原信号的标准差较全局阈值法较大），但降噪信号更光滑，而其信号发展初期的高频特性也几乎不受影响，最大限度地反映了原信号本身的特性。

    ·通过计算每一层的能量，可以得到如下的仿真结果：