数字信号处理期末报告——基于matlab的消除肌电信号工频干扰的陷波器的设计.docx

数字信号处理——基于matlab的消除肌电信号工频干扰的陷波器的设计一、课题题目基于matlab的消除肌电信号工频干扰的陷波器的设计二、课题背景表面肌电信号（sEMG）是由肌肉兴奋时所募集的运动单位产生的一个个动作电位系列在皮肤表面电极叠加而成的，是一种非平稳的微弱信号肌电信号变化与肌肉活动的水平和功能状态有关，并且它的检测具有非损伤性和多位点同步测量的优点，因此常被作为肌肉功能评价和肌肉活动控制研究的重要生理学指标而肌电信号信号弱、噪声强、频率范围低、随机性强且非平稳的特点要求我们在表面机电检测中，需要考虑的主要问题就是消除噪声干扰，提高信号的保真度表面信号检测中的噪声主要来自于电子器件的固有噪声和EMG信号本身在0~20Hz范围内受发放率的影响而在幅度上呈现出的固有的不稳定性，也会受电台、电视发射，荧光灯和电力线造成的环境噪声的干扰这些噪声或者干扰都可以通过合理的电路设计以及检测前的皮肤处理等消除或减弱影响但是，市电产生的电磁场会通过人体分布电容或者连接导线将50Hz及其高次谐波干扰引入到表面肌电信号中，而表面肌电信号的频率也主要集中在50~150Hz，这种50Hz工频干扰恰巧落在了信号的主要频带内，（且幅度可达EMG的1~3个数量级），无法轻易消除。

三、理论基础：IIR数字滤波器的系统函数可以表示为：无限冲激响应数字滤波器的设计方法： 给定数字滤波器技术指标 -> 转换成模拟滤波器技术指标 -> 设计模拟低通滤波器->转化成其他类型的模拟滤波器，得到H(s) -> 转化得到数字滤波器H(z)而这里我们选择了butterworth（通带最平坦）滤波器：，其中C为待定常数，N为待定的滤波器阶次通过通带截止频率Wp、阻带截止频率Ws、通带允许的最大衰减、阻带内应达到的最小衰减的计算（），可以得到滤波器的阶次N和3dB频率Wn当需要设计的是一个带阻滤波器时，Wp和Ws会由一个1x2的向量表示，需要将其进行归一化（即除以fs/2，fs为采样频率），便于计算机识别四、计算机仿真：（1）仿真实验1：原信号：；干扰信号：；信号长度N=200;采样频率fs=200Hz；陷波器信息：Wp=(40~60)/100;Ws=(49~51)/100;Rp=3dB;Rs=20dB经过实验可得以下图像：通过快速傅里叶变换FFT后的图，可以很明显地看到混合后的信号在40、50、60Hz的地方都有明显的峰值（由于FFT结果的对称性，通常我们只使用前半部分的结果），而经过滤波器后的时域对比可以看出，50Hz处的信号完全被滤去，但其他频率对应的信号因为过渡带的存在也有所衰减，但通过陷波器的信号在整体走向上与原信号保持了高度一致。

2）仿真实验2：原信号：随机产生，值在0~3之间；干扰信号：；信号长度N=200;采样频率fs=200Hz；陷波器信息：Wp=(40~60)/100;Ws=(49~51)/100;Rp=3dB;Rs=20dB考虑到原信号本身在50Hz的地方有可能会有成分，而通过陷波器后会将信号本身的成分消除，我们在50Hz陷波器的基础上引入了谱内插法这种方法需要建立在假设信号的功率谱在工频及其谐波位置处与其相邻的频率成分为连续变化的基础上，其思路可以表示为对信号进行傅里叶变换，计算平均功率谱P根据需要估算的频率两边的点、对应的功率值，内插得到，其中是频率分辨率信号傅里叶变换后的滤去的部分的幅值用替换，相位保持不变进行傅里叶反变换，即得到消除工频干扰后的信号经过实验可得以下图像：由于该实验中数据变较大，而且可以看出过渡带也较大，50Hz周围很大一段频率的信号都被削弱了，这时内插法最好使用数据拟合多项式的方式，用高次方程来近似原来的函数，已得到所需要的点对应的函数值但是这里我们只用了线性内插，通过两边的点直接对中间的点做插值，所以效果其实不是很好，这样做反变换后得到的也与未插值的变化不大五、真实肌电信号处理：信号信息：胫骨前肌肌电信号N=10000;fs=1000Hz陷波器信息：Wp=(1~100)/500;Ws=(49~51)/500;Rp=3dB;Rs=2dB通过对肌电信号进行滤波和谱内插的不同处理，可以得到比较明显的对比图。

可以看出，谱内插法后呈现出的信号变化更接近原信号六、课题结论：在matlab平台上，可以通过直接程序法快捷有效地完成滤波器的设计，从滤波前后的波形和频谱分析，50Hz陷波器可以很好地消除工频干扰在处理肌电信号的工频干扰时，对于50Hz陷波器，由于他在消除工频干扰的同时去掉了emg的信号成分，所以只能运用于简单的生物反馈，而谱内插法不但可以消除工频干扰，并且多EMG信号几乎没有衰减，是抑制工频干扰的有效措施附一： 参考文献：[1]梁津国等.基于自适应滤波器的表面肌电信号 消噪方法研究[J].中国医学物理杂志，2008，25卷3期：679~681[2]雷培源等.基于三级滤波器的表面肌电信号降噪处理[J].北京生物医学工程，2011，30卷1期：62~66[3]陈胤等.表面肌电信号去噪研究[J].Computer Era，2008，6期：22~24[4]邢国泉.消除50Hz工频干扰数字滤波器的设计[J].医疗卫生装备，2008，29卷12期[5]丁祥峰等.表面肌电检测中消除工频干扰的方法[J].北京生物医学工程，2006，25卷1期：63~66[6]戚仕涛等.基于MATLAB的工频干扰陷波器设计[J].医疗设备信息，2005，20卷3期：7~9附二：实验部分代码clcclearN=512;fs=200;n=0:199;t=n/fs;for i=1:200 x(i)=rand()*3;end%x=sin(2*pi*40*t)+sin(2*pi*60*t);drift=sin(2*pi*50*t);x1=x+drift;figure(1)subplot(211)plot(x);title('原信号')subplot(212)plot(x1);title('加入50Hz的正弦信号干扰后的信号')axis([-inf inf 0 5])X=fft(x,N);Y=X.*conj(X)/N;X1=fft(x1,N);Y1=X1.*conj(X1)/N;f=(0:N/2-1)*fs/N;figure(2)subplot(211)plot(f(1:N/2),abs(Y(1:N/2)));title('原信号经过傅立叶变换')xlabel('频率f/Hz')ylabel('·幅值/db')axis([20 70 0 3])subplot(212)plot(f(1:N/2),abs(Y1(1:N/2)));title('加入干扰后信号的傅里叶变换结果')xlabel('频率f/Hz')ylabel('幅值/db')axis([20 70 0 3]) wp=[40 60]/100;ws=[49 51]/100;rp=3;rs=20;[n,wn]=buttord(wp,ws,rp,rs);[h]=butter(n,wn,'stop');figure(3)freqz(h,N,fs);title('陷波器信息') sf=filter(h,1,x1);figure(4)subplot(211)plot(x);title('原信号')subplot(212)plot(sf)title('滤波后的信号')axis([-inf inf 0 5])SF=fft(sf,N);SF1=SF.*conj(SF)/N;figure(5)plot(f(1:N/2),abs(SF1(1:N/2)));title('滤波后的频谱')xlabel('频率f/Hz')ylabel('幅值/db')axis([20 70 0 3]) a=SF(1:N/2);P=real(a).^2+imag(a).^2;figure(6)subplot(211)plot(f,P)title('滤波后的能量频谱')axis([30 70 0 300])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(116:144)=PT(1:29);subplot(212)plot(f,P)title('对能量频谱做线性内插')axis([30 70 0 300])SF(116:144)=sqrt(P(116:144)).*(cos(angle(SF(116:144)))+sin(angle(SF(116:144)))*i);SF(N/2+2:N)=abs(fliplr(SF(2:N/2))).*(cos(angle(SF(N/2+2:N)))+sin(angle(SF(N/2+2:N)))*i);answer=ifft(SF,N);figure(7)subplot(311),plot(t,x)subplot(312),plot(t,sf)subplot(313),plot(t,answer(1:200))。