数字信号处理期末报告——基于matlab的消除肌电信号工频干扰的陷波器的设计.docx
数字信号处理基于mat1ab的消除肌电信号工频干扰的陷波器的设计一、课题题目基于mat1ab的消除肌电信号工频干扰的陷波器的设计二、课题背景表面肌电信号(sEMG)是由肌肉兴奋时所募集的运动单位产生的一个个动作电位系列在皮肤表面电极叠加而成的，是一种非平稳的微弱信号。肌电信号变化与肌肉活动的水平和功能状态有关，并且它的检测具有非损伤性和多位点同步测量的优点，因此常被作为肌肉功能评价和肌肉活动控制研究的重要生理学指标。而肌电信号信号弱、噪声强、频率范围低、随机性强且非平稳的特点要求我们在表面机电检测中，需要考虑的主要问题就是消除噪声干扰，提高信号的保真度。表面信号检测中的噪声主要来自于电子器件的固有噪声和EMG信号本身在020Hz范围内受发放率的影响而在幅度上呈现出的固有的不稳定性，也会受电台、电视发射，荧光灯和电力线造成的环境噪声的干扰。这些噪声或者干扰都可以通过合理的电路设计以及检测前的皮肤处理等消除或减弱影响。但是，市电产生的电磁场会通过人体分布电容或者连接导线将50HZ及其高次谐波干扰引入到表面肌电信号中，而表面肌电信号的频率也主要集中在50150Hz,这种50HZ工频干扰恰巧落在了信号的主要频带内，(且幅度可达EMG的13个数量级)，无法轻易消除。三、理论基础：IIR数字滤波器的系统函数可以表示为：MNbkZ4A=I无限冲激响应数字滤波器的设计方法：给定数字滤波器技术指标->转换成模拟滤波器技术指标->设计模拟低通滤波器->转化成其他类型的模拟滤波器，得到H(S)->转化得到数字滤波器H(z)o而这里我们选择了butterworth(通带最平坦)滤波器：I"C)F=4F其中C为待定常数，N为待定的滤波器阶次通过通带裁止频率Wp、阻带截止频率Ws、通带允许的最大衰减%、阻带内CCp=20g”.)J-201gH(e网)PH(esp),=20IgI"""=-20IgH(%)应达到的最小衰减%的计算("(*),可以得到滤波器的阶次N和3dB频率Wno当需要设计的是一个带阻滤波器时，Wp和Ws会由一个1x2的向量表示，需要将其进行归一化(即除以fs2,fs为采样频率)，便于计算机识别。四、计算机仿真：(1)仿真实验1:原信号.X=sin(8Om)+sin(120R).干扰信号.drift=Sin(IoO加).信号长度N=200;采样频率fs=200Hz;陷波器信息:Wp=(40"60)100;Ws=(49"51)100;Rp=3dB;Rs=20dB经过实验可得以下图像：原信号210-1-2(BP)-s-6-2。OCOS史守P)seqd02040608010012014016018020050-1OO-150陷淞器信息0.10.20.30.40.50.60.70.80.9Norma1izedFrequency(XKradsamp1)-1OOO0.10.20.30.4.60.60.7O.0.9Norma1izedFrequency(rad/samp1e)50O-50通过快速傅里叶变换FFT后的图，可以很明显地看到混合后的信号在40、50、60Hz的地方都有明显的峰值(由于FFT结果的对称性，通常我们只使用前半部分的结果)，而经过滤波器后的时域对比可以看出，50HZ处的信号完全被滤去，但其他频率对应的信号因为过渡带的存在也有所衰减，但通过陷波器的信号在整体走向上与原信号保持了高度一致。(2)仿真实验2:原信号：随机产生，值在03之间；干扰信号：drift=Sin(KXto).信号长度N=200;采样频率fs=200Hz;陷波器信息：Wp=(40"60)/100；Ws=(4951)/100；Rp3dB；Rs=20dB考虑到原信号本身在50HZ的地方有可能会有成分，而通过陷波器后会将信号本身的成分消除，我们在50HZ陷波器的基础上引入了谱内插法。这种方法需要建立在假设信号的功率谱在工频及其谐波位置处与其相邻的频率成分为连续变化的基础上，其思路可以表示为对信号进行傅里叶变换，计算平均功率谱P根据需要估算的频率“。两边的点口、/一1口对应的功率值，内插得到P(Go),其中"G是频率分辨率。信号傅里叶变换后的滤去的部分的幅值用Jp(豌)替换，相位保持不变。进行傅里叶反变换，即得到消除工频干扰后的信号。经过实验可得以下图像：后佰与力入印IX的正G玄佰亭干扰后门口E鼻S>r-1r-r-r-r-r-1.4-;此向岬刑郴20-4060SDIOD1201401BO180200jr-g.1OD40eoOToO-120QdOTBOQSo2005后的EW原信号经过得或叶交换303540455055606670算率创2加入干扰后僖号的傅里叶变换结果30354045505560S570250O20406080100120140160180200由于该实验中数据变较大，而且可以看出过渡带也较大，50HZ周围很大一段频率的信号都被削弱了，这时内插法最好使用数据拟合多项式的方式，用高次方程来近似原来的函数，已得到所需要的点对应的函数值。但是这里我们只用了线性内插，通过两边的点直接对中间的点做插值，所以效果其实不是很好，这样做反变换后得到的也与未插值的变化不大。五、真实肌电信号处理：信号信息：胫骨前肌肌电信号N=10000Jfs=IOOOHz陷波器信息：WP=(I100)/500；WS=(49"51)500；Rp=3dB;Rs=2dB信号的时域图信号的频域图4°050100150200250300350频率f/Hz400450500-500陷波器信息0.20.30.40.50.60.70.8Norma1izedFrequency(×radsamp1)0.9200SB8p)8詈d100-10000.10.20.30.40.50.60.70.80.91Norma1izedFrequency(rad/samp1e)流波后的频谱频率f/Hz原信号通过对肌电信号进行滤波和谱内插的不同处理，可以得到比较明显的对比图。可以看出，谱内插法后呈现出的信号变化更接近原信号。六、课题结论：在mat1ab平台上，可以通过直接程序法快捷有效地完成滤波器的设计，从滤波前后的波形和频谱分析，50Hz陷波器可以很好地消除工频干扰。在处理肌电信号的工频干扰时，对于50HZ陷波器，由于他在消除工频干扰的同时去掉了emg的信号成分，所以只能运用于简单的生物反馈，而谱内插法不但可以消除工频干扰，并且多EMG信号几乎没有衰减，是抑制工频干扰的有效措施O附一：参考文献：1梁津国等.基于自适应滤波器的表面肌电信号消噪方法研究J.中国医学物理杂志，2008,25卷3期：6796812雷培源等.基于三级滤波器的表面肌电信号降噪处理J.北京生物医学工程，2011,30卷1期：62663陈胤等.表面肌电信号去噪研究J.COn1PUterEra,2008,6期:22244邢国泉.消除50Hz工频干扰数字滤波器的设计J.医疗卫生装备，2008,29卷12期5丁祥峰等.表面肌电检测中消除工频干扰的方法J.北京生物医学工程,2006,25卷1期：63666戚仕涛等.基于MAT1AB的工频干扰陷波器设计J.医疗设备信息，2005,20卷3期：79附二：实验部分代码c1cc1earN=512;fs=200;n=0:199;t=nfs;fori=1:200x(i)=rand()*3;end%x=sin(2*pi*40*t)+sin(2*pi*60*t);drift=sin(2*pi*50*t);x1=x+drift;subp1ot(211)p1ot(x);tit1e(,原信号，)subp1ot(212)p1ot(x1);tit1e('加入50HZ的正弦信号干扰后的信号Daxis(-infinf05)X=fft(x,N);Y=X.*conj(X)/N;X1=fft(x1,N);Y1=X1.*conj(X1)/N;f=(0:N/2-1)*fsN;subp1ot(211)p1ot(f(1:N/2),abs(Y(1:N/2);tit1e('原信号经过傅立叶变换，)XIabeI(,频率f/Hz，)y1abe1('幅值db')axis(207003)subp1ot(212)p1ot(f(1:N/2),abs(Y1(1:N/2);tit1e('加入干扰后信号的傅里叶变换结果)XIabe1('频率f/Hz，)y1abe1(')1ffidb*)axis(207003J)wp=4060100;ws=4951)/100;rp=3;rs=20;n,wn=buttord(wp,ws,rp,rs);h=butter(n,wn,'stop');freqz(h,N,fs);tit1e(陷波器信息Dsf=fi1ter(h,1,x1);subp1ot(211)p1ot(x);tit1e(,原信号，)subp1ot(212)p1ot(sf)tit1e(，滤波后的信号，)axis(-infinf05)SF=fft(sf,N);SF1=SF.*conj(SF)/N;p1ot(f(1:N/2),abs(SF1(1:N/2);tit1e(，滤波后的频谱，)XIabe1('频率f/Hz，)y1abe1。幅值dk)axis(207003J)a=SF(1:N/2);P=rea1(a).2+imag(a).2;subp1ot(211)p1ot(f,P)tit1e(，滤波后的能量频谱，)axis(30700300)ft(1:115)=f(1:115);ft(116:227)=f(145:256);pt(1:115)=P(1:115);pt(116:227)=P(145:256);PT=interp1(ft,pt,f(116:144);P(II6:144)=PT(I:29);subp1ot(212)p1ot(f,P)tit1e(,对能量频谱做线性内插,)axis(30700300)SF(116:144)=sqrt(P(116:144).*(cos(ang1e(SF(116:144)+sin(ang1e(SF(116:144)*i);SF(N/2+2:N)=abs(fIip1r(SF(2:N/2).*(cos(ang1e(SF(N2+2:N)+sin(ang1e(SF(N/2+2:N)*i);answer=ifft(SF,N);subp1ot(311),p1ot(t,x)subp1ot(312),p1ot(t,sf)subp1ot(313),p1ot(t,answer(1:200)
