陷波器设计完整版本

PAGEPAGEII数字信号处理课程设计报告书题目陷波器设计课程设计任务书班级电信111姓名小丢题目陷波器设计技术参数、设计要求、检测数据等在信号采集时，往往受到50Hz电源频率干扰，分别用FIR和IIR滤波器设计方法设计实现一个50Hz陷波器，要求陷波器阻带在以内，给定采样频率，（1）设计一个50Hz的陷波器，给出原理框图，用多频信号（含有50Hz）作为输入，输出为滤去50Hz单频干扰之后的信号；（2）要求设计FIR和IIR滤波器分别至少采用2种方案实现；（3）分析比较各种方法所设计陷波器性能。设计进度安排或工作计划2014.7.1~2014.7.2：熟悉课题，查询相关资料，完成方案选择。2014.7.3~2014.7.6：设计模块划分、实现及各模块调试、验证。2014.7.7~2014.7.8：设计整体实现、调试及验证，并开始撰写报告。2014.7.9~2014.7.10：设计完成，课程设计报告撰写并定稿，上交。其它认真阅读数字信号处理课程设计报告撰写规范；课题小组经协商好要指定组长并明确分工，形成良好团队工作氛围；基于课题基本要求，各小组课再细化、增加要求；课题小组每成员均需各自撰写一份课程设计报告。课题题目摘要随着数字技术的发展，数字滤波器在许多领域得到广泛的应用。它是通信、语言、图像、自动控制、雷达、航空航天、生物医学信号处理等领域中的一种基本处理部件，具有稳定性好、精度高、灵活性大等突出优点。在信号采集时，往往受到50Hz电源频率干扰，尤其是在供电系统不稳定、外界环境适应性差时严重影响要采集信号的正确判断。本设计研究一种在MATLAB语言环境下分别用IIR和FIR滤波器设计方法设计实现一个数字陷波器，并将设计的滤波器应用到混合的正弦信号，通过仿真测试，用两种方法设计的滤波器可以很好的消除50Hz的工频干扰，并分析比较了各种方法所设计的陷波器性能。在设计IIR数字陷波器过程中,是用椭圆数字陷波滤波器的设计方法，而FIR数字陷波器的设计主要用窗函数法、频率采样法及等波纹逼近法。FIR滤波器可以得到严格的线性相位，但它的传递函数的极点固定在原点，只能通过改变零点位置来改变性能，为了达到高的选择性，必须用较高的阶数，对于同样的滤波器设计指标，FIR滤波器要求的阶数可能比IIR滤波器高5～10倍。IIR滤波器的设计可以利用模拟滤波器的许多现成的设计公式、数据和表格，计算的工作量较小。关键词数字陷波器；50Hz工频干扰；IIR和FIR滤波器目录课程设计任务书 I摘要 II1 设计概述 11.1 设计背景 11.2 设计目的 11.3 设计内容及要求 12 设计方案及实现 22.1 总体方案设计 22.2 设计原理 32.2.1 数字陷波器原理 32.2.2 IIR数字陷波器原理 32.2.3 FIR数字陷波器原理 43 设计结果分析 93.1 IIR数字陷波器设计 93.2 FIR数字陷波器设计 113.2.1 用窗函数法设计陷波器 113.2.2 频率采样法设计陷波器 133.2.3 基于切比雪夫等波纹逼近法 144 总结 17PAGE22设计概述设计背景在我国采用的是50Hz频率的交流电，所以在平时需要对信号进行采集处理和分析时，常会存在50Hz的工频干扰，对我们的信号处理造成很大干扰，因此50Hz陷波器在日常成产生活中被广泛应用，其技术已基本成熟。

陷波器也称带阻滤波器，它能在保证其他频率的信号不损失的情况下，有效的抑制输入信号中某一频率信息。所以当电路中需要滤除存在的某一特定频率的干扰信号时，就经常用到陷波器。50Hz陷波器不仅在通信领域里被大量应用，还在自动控制、雷达、声纳、人造卫星、仪器仪表测量及计算机技术等领域有着广泛的应用。鉴于MATLAB是一个具有强大分析能力的数学工具，本设计采用基于MATLAB的程序设计的方法来实现IIR数字陷波器和FIR数字陷波器，并分别通过仿真测试，实现50Hz工频信号的消除。设计目的用IIR和FIR滤波器分别设计实现一个数字陷波器，并将设计的滤波器应用到混合的正弦信号，通过仿真测试，消除50Hz的工频干扰，分析比较各种方法所设计陷波器性能。设计内容及要求深入了解利用MATLAB设计IIR和FIR数字滤波器的基本方法，在课程设计的过程中掌握程序编译及软件设计的基本方法，提高自己对于新知识的学习能力及进行实际操作的能力，锻炼自己通过网络及各种资料解决实际问题的能力。分别用FIR和IIR滤波器设计方法设计实现一个50Hz陷波器，要求陷波器阻带在以内，给定采样频率，（1）设计一个50Hz的陷波器，给出原理框图，用多频信号（含有50Hz）作为输入，输出为滤去50Hz单频干扰之后的信号；（2）要求设计FIR和IIR滤波器分别至少采用2种方案实现；（3）分析比较各种方法所设计陷波器性能。

设计方案及实现总体方案设计正弦信号正弦信号50Hz干扰信号混合信号源IIR滤波器滤波FIR滤波器滤波频谱分析滤波后信号频谱分析Butterwoth、Chebyshev、Ellipse滤波器窗函数、频率采样、等波纹逼近法图2-1数字陷波器总体方案设计原理框图设计原理数字陷波器原理陷波滤波器是一种特殊的带阻滤波器，当带阻滤波器的阻带很窄时就成了陷波滤波器。其阻带在理想情况下只有一个频率点，因此也被称为点阻滤波器。陷波器设计的目的是滤去一定范围的频率干扰，但在滤波过程中要求基本不改变其他频率成份，因此，要求陷波器的频率响应为一个理想的定点陷波器的频率响应要在消除的信号频率点处，其值为零；而在其他频率处，其值为l，如图2-2所示。图2-2陷波器的频率响应IIR数字陷波器原理IIR滤波器设计思想是：利用已有的模拟滤波器设计理论，首先根据设计指标设计一个合适的模拟滤波器，然后再通过脉冲响应不变法或双线性变换法，完成从模拟到数字的变换。常用的模拟滤波器有巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Ellipse)滤波器等，这些滤波器各有特点，供不同设计要求选用。滤波器的模拟数字变换，通常是复变函数的映射变换，也必须满足一定的要求。相比于巴特沃斯和切比雪夫I型、II型滤波器，椭圆滤波器可以获得对理想滤波器幅频响应的最好逼近，是一种性价比最高的滤波器，而且椭圆滤波器可以获得比较低的滤波器阶数，故本设计采用椭圆滤波器设计IIR数字陷波器。Matlab中ellipord和ellip函数默认采用双线性变换法[1]。由于数字滤波器传输函数只与频域的相对值有关，故在设计时可先将滤波器设计指标进行归一化处理。IIR数字滤波器的步骤：（1）将设计指标归一化处理。由于要求设计的陷波器阻带在以内，故本设计的通阻带截止频率定为：wp1=48;wp2=52;ws1=49.9;ws2=50.1;经过归一化处理后，得到：wp=[wp1,wp2]/(fs/2);ws=[ws1,ws2]/(fs/2);（2）根据归一化频率，确定最小阶数N和频率参数Wn。可供选用的阶数选择函数有:buttord，cheblord，cheb2ord，ellipord等。[N,wc]=ellipord(wp,ws,rp,rs,’z’);（3）运用最小阶数N设计模拟低通滤波器原型。根据最小阶数直接设计模拟低通滤波器原型，用到的函数有：butter,chebyl,cheby2和ellip。本设计采用椭圆滤波器，[B,A]=ellip(N,Wc,'type')设计'type'型椭圆(Elliptic)滤波器filter。N为滤波器阶数，Wc为截止频率，type决定滤波器类型，ftype=stop，设计带阻IIR滤波器[2]。（4）再用freqz函数验证设计结果。FIR数字陷波器原理一般FIR滤波器的设计方法有窗函数法、频率抽样法、切比雪夫等波逼近法等。（1）窗函数法设计陷波器原理：50Hz陷波器由截止频率为48Hz的FIR低通滤波器和截止频率为52Hz的FIR高通滤波器组成，而高通滤波器又是由全通滤波器减去截止频率为52Hz的低通FIR滤波器所得，所以FIR低通滤波器的设计是关键。使用窗函数设计FIR低通滤波器的基本思想为：首先选择一个适当的理想低通滤波器(它总是具有—个非因果、无限持续时间脉冲响应)，然后截取(或加窗)它的脉冲响应得到线性相位和因果FIR滤波器，因此这种方法的重点是选择—个合适的窗函数和理想低通滤波器。表示理想的低通滤波器，在通带上具有单位增益和线性相位，在阻带上具有零响应，表示形式为：（1-1）其中是截止频率，这个滤波器的脉冲响应具有无限持续时间，（1-2）注意关于对称，这对线性相位滤波器非常有用的。为了从得到一个FIR滤波器，必须同时在两边截取，要得到一个长度为M的因果线性相位滤波器，必须有（1-3）这种操作叫做“加窗”，一般地讲，可看作是与窗函数的乘积。（1-4）根据上述的不同定义，可以得到不同的窗结构。例如，式（1-9）为原先定义的矩形窗：（1-5）在频域中，因果FIR滤波器响应由和窗响应的周期卷积得到，如式（1-6）所示。 （1-6）Kaiser窗是最接近最优化窗结构的窗函数，因此选用其作为FIR滤波器的窗函数。通过Matlab调用“”可得到该窗函数。其中N为窗口函数的长度，值确定了窗函数的主瓣宽度和旁瓣衰减直接的关系，为返回的窗函数系数。改变N和就能改变滤波器的阻带衰减和过度带宽[3]。相反，根据对阻带衰减和过渡带带宽的要求，也能确定N和的值，它们之间有一套具体的计算公式。Kaiser窗的函数为：，其余。（1-7）式中为主瓣宽度参数，由下式计算：（1-8）可以控制窗的形状。一般加大，主辨加宽，旁辨幅度减小，典型取值为。（1-9）为零阶第一类贝塞尔函数，可用下面级数计算：（1-10）根据具体的情况N可由式N=ceil((rs-7.95)/2.285/Bt)+1求得，=7.865。选定窗函数后要构造理想低通滤波器，理想低通滤波器的脉冲响应是Sinc函数，如式(1-11)所示：（1-11）其中的主要参数是截止频率。本设计使用表示截止频率是,点数为n的理想低通滤波器。由于50Hz陷波器由截止频率为48Hz的FIR低通滤波器和通带频率为52Hz的FIR高通滤波器组成，而高通滤波器是由全通滤波器减去低通FIR滤波器所得[4]。所以理想陷波器可以表示成：（1-12）信号采集系统采样频率是400Hz，归一化后的。由于理想陷波器是非因果系统，需要使用已选的窗函数截断。利用下面的语句对该滤波器进行加Kaiser窗处理：（1-13）notch即为满足要求的50Hz陷波滤波器的冲击响应序列，它的数值即可作为线性相位50HzFIR陷波器的滤波系数。

（2）频率抽样法设计原理：周期序列的离散傅里叶级数的系数的值和的一个周期的Z变换在单位圆的N个均分点上的抽样值相等，这就实现了频域的抽样。通过有线的频率特性取样值去逼近理想的滤波特性，然后由有限的频率特性取样值（如系统冲击响应的DFT）取得系统系数。窗函数设计法是从时域出发，把理想的用一定形状的窗函数截取成有限长的，以来近似，从而使频率相应函数近似理想频率响应。频率取样法是从频域出发，对理想的频率相应进行等间隔取样[5]，以有限个频率响应去近似理想频率响应。即：，k=0,1,2,…,N-1（1-14）相比较而言，窗函数法思路为：；频率抽样法思路为：。对于应用系统函数得到插值公式：（1-15）当采样点数N已知后，便是常数，只要采样值确定，则系统函数就可确定，要求的FIR滤波器就设计出来了。对滤波系统的频率特性有：（1-16）其中是内插函数。重构的频率响应在采样点上严格等于，而在采样点之间，频率相应有加权的内插函数延伸叠加而成。设计步骤：由于设计的陷波器的过渡带要尽量窄，将过渡带宽度设置为。选取频率幅度特性的样本点。样本点数N的选取主要考虑过渡带宽，K为过渡带点数（1-17）频率取样法所得的FIR滤波器通带和阻带波动主要是由于过渡带的突变引起的。通过优化过渡带样点可以得到较好的通阻带特性，缺点是过渡带加宽，但可以通过增加样点数N来克服。一般1个过渡带样点可以使最小阻带衰减达-40dB，而2个过渡带样点可以使最小阻带衰减达-60dB，为了得到更好的设计效果，因此本设计取了2个过渡带样点。将和K=2带入式（1-17）可确定抽样点数N=750。（3）基于切比雪夫等波纹逼近法设计原理：窗口法设计和频率采样设计都存在某些缺陷。首先，在设计中不能将边缘频率和精确地给定。其次，不能够同时标定波纹因子δ1和δ2。最后，近似误差（即理想响应和实际响应之间的差）在频率区间上不是均匀分布的。而最优等纹波设计法能解决上面三个问题。等波纹切比雪夫逼近准则就是通过通带和阻带使用不同的加权函数，实现在不同频段（通常指的是通带和阻带）的加权误差最大值相同，从而实现其最大误差在满足性能指标的条件下达到最小值，即使得和之间的最大绝对误差最小。等波纹切比雪夫逼近是采用加权逼近误差，它可以表示为：（1-18）其中，为逼近误差加权函数在误差要求高的频段上，可以取较大的加权值，否则，应当取较小的加权值。尽管按照FIR数字滤波器单位取样响应的对称性和的奇、偶性，FIR数字滤波器可以分为4种类型，但滤波器的频率响应可以写成统一的形式：（1-19）其中，为幅度函数，且是一个纯实数，表达式也可以写成统一的形式：（1-20）其中，为的固定函数，为M个余弦函数的线性组合。设计步骤：利用数字信号处理工具箱中的remezord和remez函数可以实现FIRDF的最优化设计。在此先介绍这两个函数：功能:利用remezord函数可以通过估算得到滤波器的近似阶数M，归一化频率带边界fo，频带内幅值mo各个频带内的加权系数w。输入参数f为频带边缘频率，m为各个频带所期望的幅度值，rip是各个频带允许的最大波动，fs为采样频率。功能:利用remez函数可以得到最优化设计FIR滤波器的系数，输入参数M是滤波器的阶数，f，m，w参数含义说明同(1)。实际设计中，由于remezord函数可跑高估或低估滤波器的阶数，因此在得到滤波器系数后，必须检查其阻带最小衰减是否满足设计要求。如果此时的技术指标不能满足设计要求，则必须提高滤波器的阶数到等。设计结果分析IIR数字陷波器设计图3-1IIR陷波器幅频特性和衰减特性曲线图中，利用椭圆带阻滤波器无限缩短阻带宽度趋近陷波器，从幅频特性曲线中可以看出，输入信号频率小于48Hz或大于52Hz信号被通过，而频率在49.5~51.5Hz之间的信号被截止，通阻带内信号产生的跃变波形抖动是椭圆滤波器的通阻带波纹抖动引起的。从衰减特性曲线中可以看出，通带衰减为-0.9dB，阻带衰减为-100dB,满足滤波要求。图3-2滤波前后信号波形图由图，输入信号为分别含有频率为48Hz、50Hz、52Hz的正弦信号，其时域波形为含有3段不同包络的正弦信号，滤波以后不含频率为50HZ信号的包络了。图3-3输入输出信号频谱图由图易知，对于分别含有频率为48Hz、50Hz、52Hz的输入信号，在其对应的频谱有三根谱线，经椭圆滤波器滤波后，48HZ和52HZ的谱线幅度均变为原来的95%，50HZ的谱线长度几乎为0，而参数要求的通带衰减由公式计算得到已达到-0.9dB和阻带衰减-100dB均能满足设计要求。还可以看出，用椭圆滤波器设计的陷波器，旁瓣峰值很小，过渡带的宽度也尽量设计到很小。FIR数字陷波器设计用窗函数法设计陷波器图3-4FIR陷波器幅频特性和衰减特性曲线从图3-4可以看出，利用Kaiser窗设计的FIR滤波器无限缩短阻带宽度趋近陷波器，从幅频特性曲线中可以看出，输入信号频率小于48Hz或大于52Hz信号被通过，而频率在50Hz附近很小的范围内的信号被截止。从衰减特性曲线中可以看出，在陷波器中心频率50Hz处衰减接近-100dB。

图3-5频率为48、50、52Hz混合正弦信号滤波前后波形从图3-5可以知道当输入信号为分别含有频率为48Hz、50Hz、52Hz的混合正弦信号，其时域波形为含有3段不同包络的正弦信号，滤波以后不含频率为50Hz信号的包络了。但是由于滤波器的阶数比较高，因此有比较大的“群延时”，即每个点都向后延迟了相同的一段时间，“群延时”正比于滤波器的阶数，阶数越小滤波器的延迟越小，但滤波性能会随之下降，是每个FIR滤波器不可避免的问题。故本设计在首先达到滤波性能要求的前提下，尽量减小阶数。图3-6混合信号滤波前后频谱图从频谱图可以看出对于分别含有频率为48Hz、50Hz、52Hz的输入信号，其对应的幅频特性有三根谱线，经FIR滤波器滤波后，48Hz和52Hz的谱线均变为原来的87.5%，通带衰减也小于3dB，50HZ的谱线长度为15，变为原来的0.75%，阻带衰减达到-42.5dB满足要求。频率采样法设计陷波器图3-7频率采样法设计FIR陷波器幅频特图3-8频率采样法设计FIR陷波器损耗特性曲线从图可以看出，利用频率采样法设计的FIR滤波器也能满足设计要求，从幅频特性曲线中可以看出，输入信号频率小于48Hz或大于52Hz信号均能被通过，而频率在50Hz附近很小的范围内的信号被截止。从衰减特性曲线中可以看出，在陷波器中心频率50Hz处衰减接近-320dB。图3-948,、50、52Hz混合信号滤波前后频谱图从频谱图可以看出，当输入信号频率小于48Hz或大于52Hz通过滤波器后均没有衰减，而在50Hz信号幅度几乎为0，满足设计要求。基于切比雪夫等波纹逼近法图3-10频率采样法设计FIR陷波器幅频特性图3-11等波纹逼近法设计FIR陷波器幅频及损耗特性曲线从图可以看出，在50Hz处，信号的衰减达到-70dB完全符合设计要求，在48Hz和52Hz处信号也能完整的被保留下来，没有产生衰减。图3-1248,、50、52Hz混合信号滤波前后频谱图从频谱图可以看出，当输入信号频率小于48Hz或大于52Hz通过滤波器后均没有衰减，而在50Hz信号幅度几乎为0，满足设计要求。结论：频率采样法是直接在频域上采样，保证设计的和希望的滤波器幅度值相等，而在采样点之间是用内插函数和相乘的线性组合形成的，这样使得频域不连续点附近的误差大，且边界频率不易受到控制。而窗函数法中是用窗函数直接截取希望设计的滤波器的一段，作为滤波器的。其频率响应有一定的误差，并且由于吉布斯效应，使得过渡带附近的通带内具有较大的上冲，而阻带衰减过小,而且窗函数法在阶数较低时，阻带特性不满足设计要求，只有当滤波器阶数较高时，使用凯塞窗基本可以达到阻带衰耗要求。以上两种方法存在的共同缺陷是在边界频率的控制上不如等波纹设计精确。利用切比雪夫等波纹逼近方法设计的滤波器将所有的抽样值皆作为变量，在某个优化准则下，通过计算机的迭代运算，以得到最优结果，能够获得较好的通带和阻带性能，并能准确地指定通带和阻带边缘。由于该滤波器在通带和阻带的误差是均匀分布的，因此其频率响应在通带和阻带内显示出等波纹性，阶次可以比较低，但是等波纹逼近法设计的滤波器具有较大的通带和阻带波纹，而且硬件实现难度较大。IIR滤波器和FIR滤波器比较：1.在相同技术指标下，IIR滤波器由于存在着输出对输入的反馈，因而可用比FIR滤波器较少的阶数来满足指标的要求，这样一来所用的存储单元少，运算次数少，较为经济。例如用频率抽样法设计阻带衰减为-20db的FIR滤波器，其阶数要33阶才能达到，而如果用双线性变换法设计只需4-5阶的切贝雪夫滤波器，即可达到指标要求，所以FIR滤波器的阶数要高5-10倍左右。2.FIR滤波器可得到严格的线性相位，而IIR滤波器则做不到这一点，IIR滤波器选择性愈好，则相位的非线性愈严重，困而，如果IIR滤波器要得到线性相位，又要满足幅度滤波的技术要求，必须加全通网络进行相位校正，这同样会大大增加滤波器的阶数，从这一点上看，FIR滤波器又优于IIR滤波器。3.FIR滤波器主要采用非递归结构，因而从理论上以及时性从实际的有限精度的运算中，都是稳定的。有限精度运算误差也较小，IIR滤波器必须采用递归的结构，极点必须在Z平面单位圆内，才能稳定，这种结构，运算中的四舍五入处理，有时会引起寄生振荡。4.FIR滤波器，由于冲激响应是有限长的，因而可以用快速傅里叶变换算法，这样运算速度可以快得多，IIR滤波器则不能这样运算。5.从设计上看，IIR滤波器可以利用模拟滤波器设计的现成闭合公式、数据和表格，因而计算工作量较小，对计算工具要求不高。FIR滤波器则一般没有现成的设计公式，窗函数法只给出窗函数的计算工式，但计算通带、阻带衰衰减仍无显示表达式。一般FIR滤波器设计只有计算机程序可资利用，因而要借助于计算机。6.IIR滤波器主要是设计规格化的、频率特性为分段常数的标准低通、高通、带通、带阻、全通滤波器，而FIR滤波器则要灵活得多，例如频率抽样设计法，可适应各种幅度特性的要求，因而FIR滤波器则要灵活得多，例如频率器可设计出理想正交变换器、理想微分器、线性调频器等各种网络，适应性较广。而且，目前已有许多FIR滤波器的计算机程序可供使用。总结数字滤波器的应用十分广泛,运用MATLAB语言,能很容易地设计出适合各种要求的滤波器。在设计的过程中，我发现，运用MATLAB语言进行数字滤波器设计十分的方便、快捷。这次的程序设计不是很难，根据老师给的提示函数，在网上、书上查找一些例子，认真的编写、调试很快的就得到了老师要求的设计的滤波器，就是在把数字滤波器的性能指标转换成相应的模拟滤波器性能指标上，我出现了一些问题，比如分析混合正弦信号滤波前后频谱幅度变化，FIR滤波参数指标的调整，以及陷波器过渡带及其截止频率的反复计算等等，但是最后通过问老师、查资料还是解决了。通过这次设计，我加强了MATLAB语言的掌握，对于滤波器也有了深入的认识，通过查阅书籍我学到了很多新的知识，锻炼了我的自学与动手能力，为今后的工作打好了坚实的实践基础，使我受益匪浅！由于时间局限，对于MATLAB的知识还仅仅只是了解，并没有完全地掌握。我也明白作为一项工具，MATLAB还有许多需要我去学习去了解，并不是一时的了解就可以的，要学会去深究。MATLAB具有相当强大的矩阵运算和操作功能，其程序编写也比其它语言要简单，是一种对学习非常有帮助的软件，在之后的时间，我依然会对它进行学习。参考文献[1]高西全，丁玉美.数字信号处理（第三版）[M].西安：西安电子科技大学出版社，2009:140～160[2]邢国泉.消除50Hz工频干扰数字滤波器的设计［J］．医疗卫生装备，2008，27（6）:597～599[3]戚仕涛,汤黎明，吴敏，刘铁兵.基于MATLAB的工频干扰陷波器设计[J].医疗卫生装备，2005，26（8）：10-13[4]梁奇，叶明，马文杰.滤除SEMG工频干扰的数字陷波器设计[J].计算机工程与应用，2009,45（17）:61～63[5]蔡建平，黄晓红，孙丽英，朱艺.基于频率采样法的线性相位滤波器设计及Matlab[J].电力自动化设备，2006,26（7）:59～61

附录源码IIR数字陷波器：clc;closeall;wp1=48;wp2=52;ws1=49.9;ws2=50.1;rp=0.1;rs=100;fs=400;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');[hw,w]=freqz(num,den,512);figure(1);subplot(211);plot(w*fs/(2*pi),abs(hw));gridon;xlabel('频率/Hz');ylabel('幅值');title('H(e^j^w)');axis([40,60,0,1.5]);subplot(212);plot(w*fs/(2*pi),20*log10(abs(hw)));gridon;xlabel('频率/Hz');ylabel('幅值');title('损耗函数');axis([40,60,-200,20]);M=4000;n=0:M-1;T=1/fs;xn=sin(96*pi.*n*T)+sin(100*pi.*n*T)+sin(104*pi.*n*T);xk=fft(xn,M);yn=filter(num,den,xn);%复合信号通过Ellipse滤波器滤波后的信号ynk=2*(0:M-1)/M;yk=fft(yn,M);figure(1)subplot(211);plot(w/2/pi*fs,abs(hw));title('幅频特性');xlabel('f/Hz');ylabel('H(e^jw)');gridon;axis([40,60,0,1.2]);gridon;subplot(212);plot(w/2/pi*fs,20*log10(abs(hw))/max(abs(hw)));title('损耗函数');xlabel('f/Hz');ylabel('幅度/dB');axis([40,60,-120,10]);gridon;figure(2)subplot(211);stem(k,abs(xk),'.');xlabel('\omega/\pi');ylabel('幅度');title('滤波前输入信号x(n)的频谱');axis([0.2,0.3,0,2100]);gridon;subplot(212);stem(k,abs(yk),'.');xlabel('\omega/\pi');ylabel('幅度');title('Ellipse滤波后输出信号y(n)的频谱');axis([0.2,0.3,0,2100]);gridon;figure(3)subplot(211);plot(n,xn);xlabel('f/Hz');ylabel('幅度');title('原信号波形图');axis([0,1000,-4,4]);subplot(212);plot(n,yn);xlabel('f/Hz');ylabel('幅度');title('滤波后波形图');axis([0,1000,-4,4]);FIR数字陷波器：窗函数设计法：低通滤波器：functionhd=ideal_lp(wc,M);alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha+eps;%addsmallestnumbertoavoiddividedbyzerohd=sin(wc*m)./(pi*m);clc;clearall;%50Hz工频干扰陷波器%陷波器由一个低通滤波器加上一个高通滤波器组成,而高通滤波器可以用一个全通滤波器减去一个低通滤波器构成Fs=400;T=1/Fs;%采样频率和周期rs=80;%阻带最小衰减wc1=51/(Fs/2)*pi;%wc1为高通滤波器-6dB截止频率，对应51Hzwc2=49/(Fs/2)*pi;%wc2为低通滤波器-6dB截止频率，对应49HzBt=0.0086*pi;N=ceil((rs-7.95)/2.285/Bt)+1;%滤波器阶数%alph=0.1102(rs-8.7);beta=7.865;%衰减系数h=ideal_lp(pi,N)-ideal_lp(wc1,N)+ideal_lp(wc2,N);%h为陷波器冲击响应wn=kaiser(N,beta);notch=h.*wn';%notch为50Hz陷波器冲击响应序列[hw,w]=freqz(notch,1,512);M=4000;n=0:M-1;xn=sin(96*pi.*n*T)+sin(100*pi.*n*T)+sin(104*pi.*n*T);xk=fft(xn,M);yn=fftfilt(notch,xn,M);%复合信号通过凯瑟窗低通滤波后的信号yn%yn=conv(notch,xn);%yn=yn(N/2+1:N/2+M);%延迟校正k=2*(0:M-1)/M;yk=fft(yn,M);figure(1)subplot(211);plot(w/2/pi*Fs,abs(hw));title('幅频特性');xlabel('f/Hz');ylabel('H(e^jw)');gridon;axis([40,60,0,1.2]);subplot(212);plot(w/2/pi*Fs,20*log10(abs(hw))/max(abs(hw)));title('损耗函数');xlabel('f/Hz');ylabel('幅度/dB');gridon;axis([40,60,-100,10]);figure(2)subplot(211);stem(k,abs(xk),'.');xlabel('f/Hz');ylabel('幅度');title('滤波前输入信号x(n)的频谱');axis([0.2,0.3,0,2100]);gridon;subplot(212);stem(k,abs(yk),'.');xlabel('w/\pi');ylabel('幅度');title('通过凯瑟窗滤波后输出信号y(n)的频谱');axis([0.22,0.28,0,2100]);gridon;figure(3)subplot(211);plot(n,xn);xlabel('n');ylabel('幅度');title('原信号波形图');axis([0,1000,-4,4]);gridon;subplot(212);plot(n,yn);xlabel('n');ylabel('幅度');title('滤波后波形图');axis([0,1000,-4,4]);gridon;频率采样法：clc;clearall;%fs=400;%fsl=49.6;fsh=50.4;fpl=48;fph=52;%wsl=fsl*2*pi/fs;wsh=fsh*2*pi/fs;wpl=fpl*2*pi/fs;wph=fph*2*pi/fs;wp1=0.24*pi;ws1=0.248*pi;wp2=0.252*pi;ws2=0.26*pi;delta=(ws1-wp1);M=ceil(2*pi*3/delta)+1;al=(M-1)/2;wl=(2*pi/M);k=0:M-1;T1=0.11;T2=0.59;%过渡带样本点Hrs=[ones(1,ceil((0.24*pi/wl))),T2,T1,zeros(1,ceil((0.008*pi/wl)-1)),T1,T2,ones(1,ceil((1.47*pi/wl)+3)),T2,T1,zeros(1,ceil(0.008*pi/wl)-1),T1,T2,ones(1,ceil((0.24*pi/wl)))];wdl=[00.2350.2450.2550.2651];k1=0:floor((M-1)/2);k2=floor((M-1)/2)+1:M-1;angH=[-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2)];H=Hrs.*exp(j*angH);h=real(ifft(H));figure(1);stem(k,h);title('impulseresponse');xlabel('n');ylabel('h(n)');grid;figure(2);Hf=abs(H);w=k*wl/pi;plot(w,Hf,'*b-');axis([01-0.11.1]);title('amplituderesponse');xlabel('frequencyinpiuni