MATLAB课程设计报告

PAGEPAGE19燕山大学课程设计说明书题目：基于matlab的肌电信号工频干扰的消除学院（系）：电气工程学院年级专业：09级测控仪表学号学生姓名：指导教师：教师职称：燕山大学课程设计（论文）任务书院（系）：电气工程学院基层教学单位：自动化仪表系学号学生姓名专业（班级）设计题目基于matlab的肌电信号工频干扰的消除设计技术参数采样频率：fs=1000HZ干扰频率：f=50HZ（100HZ,150HZ,200HZ,250HZ…）滤波器参数：通带最大波纹Rp阻带最小衰减Rs通带边界频率wp阻带边界频率ws最小阶数N截止频率Wc设计要求设计一个梳妆滤波器，用于滤除肌电采集信号中的50Hz工频及其谐波干扰，分析所涉及滤波器的频响，比较滤波前后的时域波形和频谱特性。肌电信号的采样频率为1000Hz。工作量1对原信号进行频谱分析；2设计滤波器并对其进行频谱分析；3用滤波器处理原始信号，并进行滤波分析；4设计显示界面。工作计划周一：分析处理原始信号：周（二三）：设计滤波器并对其进行频谱分析；周四：用滤波器处理原始信号，并进行滤波分析；周五：设计显示界面并整论文。参考资料谢平、王娜、林洪斌等主编，信号处理原理及应用。北京：机械工业出版社，宋爱国、刘文波等主编，测试信号分析与处理。北京：机械工业出版社，2005,9聂祥飞、王海宝、谭泽富主编，Matlab程序设计及其在信号处理中的应用。成都：西南交通大学出版社，2005吴大正、高西全等主编，Matlab及在电子信息课程中的应用。北京：电子工业出版社，2006.3李培芳、孙晖、李江主编，信号与系统分析基础。北京：清华大学出版社，2006.12指导教师签字谢平吴晓光基层教学单位主任签字谢平2012年6月燕山大学课程设计评审意见表指导教师评语：成绩：指导教师：谢平吴晓光2012年6答辩小组评语：成绩：组长：2012年6课程设计总成绩：答辩小组成员签字：2012年6摘要肌电信号是产生肌肉力的电信号根源，它是肌肉中很多运动单元动作电位在时间和空间上的叠加，反映了神经，肌肉的功能状态，在基础医学研究、临床诊断和康复工程中有广泛的应用。梳状滤波器它是由许多按一定频率间隔相同排列的通带和阻带，只让某些特定频率范围的信号通过。梳状滤波器其特性曲线象梳子一样，故称为梳状滤波器。本文主要介绍使用梳状滤波器消除在检测和记录肌电信号时噪声和干扰的影响，如50HZ的工频及其谐波，并使用MATLAB软件对原信号、滤波器、处理信号进行频谱分析，观察其频谱特性。关键词：梳状滤波器；MATLAB；频谱分析目录第一章设计过程………61.1滤波器设计结构图………………61.2设计滤波器步骤…………………6第二章MATLAB的应用……………………62.1MATLAB语言简介………………62.2一些函数及指令简介……………7第三章程序和仿真图……………………93.1滤波器设计程序及仿真图………93.2原始信号的频谱分析……………113.3信号滤波的程序设计及频谱分析………………12第4章学习心得…………16参考文献…………………17第一章设计过程1.1滤波器设计结构图椭圆滤波器设计结构图如图所示： 采样（1000HZ）梳状滤波器连续混合采样（1000HZ）梳状滤波器信号图4.1结构框图1.2设计数字滤波器的步骤1.确定数字滤波器的性能指标：Wp、Ws、Rp、Rs、N、Wc。2.将数字滤波器的性能指标转化成相应模拟滤波器的性能指标。3.设计满足指标要求的模拟滤波器Ha(s)。4.通过变换将模拟滤波器转换成数字滤波器4.数字滤波器的MATLAB实现和频谱分析第二章MATLAB的应用2.1MATLAB语言简介MATlAB是在20世纪80年代初期，由美国的MathWorks软件开发公司正式推出的一种数学工具软件。它以矩阵运算为基础，把计算、可视化、程序设计有机地融合到了一个简单易学的交互式工作环境中，其出色的数值计算能力和强大的图形处理功能得到了广大科研人员和工程技术人员的垂青。它拥有功能全面的函数库，他把大量的函数封装起来，让用户脱离了繁琐复杂的程序设计过程，只需要将精力集中到模型的建立和工作上即可，大大提高了工作效率。利用MATLAB可以实现科学计算、符号运算、算法研究、数学建模和仿真、数据分析和可视化、科学工程绘图以及图形用户界面设计等强大功能。更重要的是，MATLAB提供了图形用户界面（GUI)设计功能，用户可以自行设计别具风格的人机交互界面。此外，MATLAB提供了与C语言或Fortran语言的交互操作，是用户能通过程序接口很方便地进行程序间的调用和转换。2.2一些函数及指令简介MATLAB的信号处理工具箱提供了设计椭圆滤波器的函数：ellipord函数、ellip函数、filter函数、freqz函数、fft函数和abs函数。1.ellipord函数的功能是求滤波器的最小阶数，其调用方式为：[n,wp]=ellipord(wp,ws,rp,rs)n-椭圆数字滤波器最小阶数；wp-椭圆滤波器通带截止角频率；ws-椭圆滤波器阻带起始角频率；rp-通带最大衰减量也即通带波纹（dB）；rs-阻带最小衰减量(dB)；这里wp、ws都是归一化频率，即0wp（或ws）1,1对应弧度。2.ellip函数的功能是用来设计椭圆滤波器，其调用方式为：（1）[b,a]=ellip(n,rp,rs,wp)rp用来指定通带内波纹的最大衰减；rs用来指定阻带内波纹的最小衰减；在这里，wp的取值范围是（0.0,1.0），其中1对应于0.5Fs，Fs为采样频率。（2）[b,a]=ellip(n,rp,rs,wp，'high')：设计高通椭圆滤波器。上述返回的向量b和a的维数都是n+1，而不是n。filter函数功能：利用IIR滤波器和FIR滤波器对数据进行滤波，其调用方式为：y=filter(b,a,x)：计算输入信号X的滤波输出，向量b和a分别是所采用的滤波器的分子系数向量和滤波器的分母系数向量。[y,zf]=filter(b,a,x,zi)：参数zi指定滤波器的初始条件值，其大小为zi=max（length（a），length（b））-1。y=filter(b,a,x,[],dim)或y=filter(b,a,x,zi,dim)：参数dim指定滤波的维数。filter采用数字滤波器对数据进行滤波，其实现采用移位直接Ⅱ型结构，因而适用于IIR和FIR滤波器。滤波器的系统函数为即滤波器系数a=[a0a1a2...an],b=[b0b1...bm],输入序列矢量为x。这里，标准形式为a0=1，如果输入矢量a时，a0≠1，则MATLAB将自动进行归一化系数的操作；如果a0=0，则给出出错信息。4.freqz函数功能：求解数字滤波器频率响应，其调用方式为：（1）[h,w]=freqz(b,a,n)：返回数字滤波器的n点复频率响应，输入参数b和a分别是滤波器系数的分子和分母向量；输出参数h是复频率响应，w是频率点。输入参数n的默认值为512。（2）h=freqz(b,a,w)：计算由向量w（单位为rad/sample，范围为[0,]）指定的频率点的复频率响应。（3）[h,f]=freqz(b,a,n,Fs)：同时输出实际频率点。用于对H(ejω)在[0,Fs/2]上等间隔采样n点，采样点频率及相应频响值分别记录在f和h中。由用户指定FS（以HZ为单位）值。（4）h=freqz(b,a,f,Fs)：计算由向量f指定的频率点的复频率响应。5.Abs函数调用方式为：Y=abs（x）：计算x各元素的绝对值。当x为一个复数的时候，则计算x的复数模。第三章程序和仿真图3.1滤波器设计程序及仿真图ws2=150.1;rp=0.3;rs=40;Fs=1000;ws2=150.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H3,W3]=freqz(num,den);%150HZ滤波wp1=195;wp2=205;ws1=199.9;ws2=200.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H4,W4]=freqz(num,den);%200HZ滤波wp1=245;wp2=255;ws1=249.9;ws2=250.1;wp1=45;wp2=55;ws1=49.9;ws2=50.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H1,W1]=freqz(num,den);%50HZ滤波wp1=95;wp2=105;ws1=99.9;ws2=100.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H2,W2]=freqz(num,den);%100HZ滤波wp1=145;wp2=155;ws1=149.9rp=0.3;rs=40;Fs=1000;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H5,W5]=freqz(num,den);%250HZ滤波wp1=295;wp2=305;ws1=299.9;ws2=300.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H6,W6]=freqz(num,den);%300HZ滤波wp1=345;wp2=355;ws1=349.9;ws2=350.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H7,W7]=freqz(num,den);%350HZ滤波wp1=395;wp2=405;ws1=399.9;ws2=400.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H8,W8]=freqz(num,den);%400HZ滤波wp1=445;wp2=455;ws1=449.9;ws2=450.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H9,W9]=freqz(num,den);%450HZ滤波H=H1+H2+H3+H4+H5+H6+H7+H8+H9;W=W1+W2+W3+W4+W5+W6+W7+W8+W9;plot(W*Fs/(18*pi),abs(H));gridon;xlabel('频率/Hz');ylabel('幅值');%梳状滤波器频谱显示滤波器仿真图3.2原始信号的频谱分析X=[原始信号，肌电信号];N=length(x);y=fftshift(fft(x));df=fs/(N-1);n=-(N-1)/2:(N-1)/2;q=n*df;subplot(2,1,1);plot(q,abs(y));title('幅频图像')subplot(2,1,2);plot(q,angle(y));title('相频图像')原始信号仿真图3.3信号滤波的程序设计及频谱分析x=load(‘EMG.txt’);wp1=45;wp2=55;ws1=49.9;ws2=50.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H1,W1]=freqz(num,den);sf1=filter(num,den,x);wp1=95;wp2=105;ws1=99.9;ws2=100.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H2,W2]=freqz(num,den);sf2=filter(num,den,sf1);wp1=145;wp2=155;ws1=149.9;ws2=150.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H3,W3]=freqz(num,den);sf3=filter(num,den,sf2);wp1=195;wp2=205;ws1=199.9;ws2=200.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H4,W4]=freqz(num,den);sf4=filter(num,den,sf3);wp1=245;wp2=255;ws1=249.9;ws2=250.1;rp=0.3;rs=40;Fs=1000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);[N,wc]=ellipord(wp,ws,rp,rs,'z');[num,den]=ellip(N,rp,rs,wc,'stop');[H5,W5]=freqz(num,den);sf5=filter(num,den,sf4);N=length(sf5);y=fftshift(fft(sf5));%fftshift的作用正是让正半轴部分和负半轴部分的图像分别%关于各自的中心对称。因为直接用fft得出的数据与频率不%是对应的df=1000/(N-1);n=-(N-1)/2:(N-1)/2;q=n*df;subplot(2,1,1);plot(q,abs(y));title('幅频图像')subplot(2,1,2);plot(q,angle(y));title('相频图像')滤波信号仿真图第四章学习心得经过一周的课程设计，我学到了很多东西。对于以前不理解的知识，通过试验的学习得到