《数字信号处理》课程设计报告-数字滤波器的设计与信号的滤波处理

《数字信号处理》课程设计报告目录TOC\o"1-4"\h\z\u一、设计题目 1二、设计目的 1三、课程设计的内容 1四、课程设计需要完成的工作 1五、设计过程与结果 15.1三组原始数据的时域图像 15.2三组原始数据的频谱图 25.3数字滤波器的设计与信号滤波 25.3.1IIR数字滤波器设计 35.3.1.1.Butterworth低通滤波器 35.3.1.2．切比雪夫I型低通滤波器 35.3.1.3．切比雪夫II型低通滤波器 45.3.1.4．椭圆低通滤波器 65.3.2FIR（窗函数法）滤波器设计 65.3.2.1．Hamming窗函数滤波器 65.3.2.2．频率采样设计法、最优等波纹设计法 75.4结果分析 7六、滤波器分析 8七、个人思考 107.1FIR滤波器的窗口设计法 107.2IIR数字滤波器设计 12八、拓展 148.1FDATool界面说明 158.2FDATool性能分析介绍 158.3FDATool滤波器系数导出 15九、收获与感想 16十、MATLAB程序代码 16《数字信号处理》课程设计报告——数字滤波器的设计与信号的滤波处理一、设计题目应用Matlab软件对在实验室采集的信号进行频谱分析并设计数字滤波器对信号进行滤波处理。二、设计目的为了巩固所学的数字信号处理理论知识，使学生对信号的采集、处理、读取存储等有一个系统的掌握和理解，学会数字滤波器的设计与应用。三、课程设计的内容采集个人的呼吸流量电压信号、左右肺电压信号，共三组信号，时长5分钟，采样率为7.812Hz，以获得处理所需信号。画出采样信号的时域波形和频谱分析图，并设计合适的数字滤波器对采集的信号数据进行滤波，画出滤波后信号的时域波形，并与滤波前信号图像作对比，分析信号的变化。并通过对比几种数字滤波器的滤波结果，分析它们各自的特点。四、课程设计需要完成的工作电压信号的采集；信号的频谱分析；数字滤波器的设计；对信号进行滤波处理；对滤波前后信号进行对比分析；对比几种数字滤波器的滤波结果，分析它们各自的特点。五、设计过程与结果5.1三组原始数据的时域图像由图像可以看出，三组数据原始数据具有较多的噪声（噪声出现的原因将在后面的信号分析结果中进行说明），从图像上难以对数据特点进行分析，因此需要进行合适的滤波后才能进行数据的分析。为了得出数//据的主频率一遍进行滤波设计，下面将数据进行FFT，画出频谱图进行分析。5.2三组原始数据的频谱图从三组数据进行FFT后的图像中可以看出，三组数据的有效频率都在0到0.5Hz之间，0.5Hz以上的可视为噪声，应选择低通滤波器对噪声加以滤除。这个结果从生理学上来说是正确的，因为根据资料显示，人呼吸的频率约在0.2Hz到0.6Hz之间，因此对信号进行频谱分析后会得出这样的图像结果。下面根据这个结果进行低通滤波器的设计。5.3数字滤波器的设计与信号滤波数字滤波器主要分为两种：FIR(FiniteImpulseResponse)滤波器（又名：“有限长单位冲激响应滤波器”）以及IIR(InfiniteImpulseResponse)数字滤波器（又名“无限脉冲响应数字滤波器”，或“递归滤波器”）。IIR数字滤波器与FIR数字滤波器的区别主要有：1、单位响应IIR数字滤波器单位响应为无限脉冲序列，而FIR数字滤波器单位响应为有限的；FIR滤波器，也就是“非递归滤波器”，没有引入反馈。这种滤波器的脉冲响应是有限的。2、幅频特性IIR数字滤波器幅频特性精度很高，不是线性相位的，可以应用于对相位信息不敏感的音频信号上；FIR数字滤波器的幅频特性精度较之于IIR数字滤波器低，但是线性相位，就是不同频率分量的信号经过FIR滤波器后他们的时间差不变，这是很好的性质。3、实时信号处理FIR数字滤波器是有限的单位响应也有利于对数字信号的处理，便于编程，用于计算的时延也小，这对实时的信号处理很重要。5.3.1IIR数字滤波器设计 5.3.1.1.Butterworth低通滤波器5.3.1.2．切比雪夫I型低通滤波器%通带波纹Rp=0.025的图像：%通带波纹Rp较小效果较好，以下给出Rp=1的图像：%从结果图可以看出，Rp越大，滤波器在通带处的波动越大。5.3.1.3．切比雪夫II型低通滤波器%阻带衰减As=60的图像： %阻带衰减As较大效果较好，以下给出As=15的图像：%从结果图可以看出，As越小，滤波器在阻带处的波动越大。5.3.1.4．椭圆低通滤波器5.3.2FIR（窗函数法）滤波器设计5.3.2.1．Hamming窗函数滤波器利用Matlab可以进行的窗函数滤波器类型还有很多，例如：（1）w=boxcar(M)在数组w中产生M点的矩形窗函数。（2）w=triang(M)在数组w中产生M点的Bartlett（三角形）窗函数。（3）w=hann(M)在数组w中产生M点的Hanning窗函数。（4）w=hamming(M)在数组w中产生M点的Hamming窗函数。（5）w=blackman(M)在数组w中产生M点的Blackman窗函数。（6）w=kaiser(M，beta)在数组w中产生M点的kaiser窗函数。（7）……这些窗函数需要根据不同的实际条件进行选择，这里就不再一一赘述了。5.3.2.2．频率采样设计法、最优等波纹设计法FIR滤波器的设计方法还有很多，例如频率采样设计法、最优等波纹设计法等，这些方法在这里就不去设计了，查阅相关书籍资料便可得知详情。5.4结果分析从信号原始图上看，在310-820信号点之间出现了比较强烈噪声干扰，特别是左肺电极电压信号，结合下面的FFT分析，噪声信号集中在2Hz左右，具体原因难以说明，据了解，这次试验中采集的信号大多都存在着这样一种现象，估计原因为以下几点：心电信号的干扰：左肺信号测试电极位于心脏处，因此更容易受到心电信号的影响，产生噪声干扰。呼吸速率的影响：由于第一次进行这个实验的测试，全程用嘴巴呼吸，有很大程度上的紧张、呼吸不习惯等原因存在，特别卡碍事一段时间后，呼吸频率会变得急促，后来经过调整便能恢复，且呼吸频率也会直接影响心跳速率，因此，信号在中段处的噪声会特别明显。测试过程中身体的移动可能导致噪声的产生。从对信号进行FFT处理后得出的图像上看，有效信号集中在0.1Hz-0.6Hz之间，这个结果从生理学上来说是正确的，因为根据资料显示，人呼吸的频率约在0.2Hz到0.6Hz之间，因此对信号进行频谱分析后会得出这样的图像结果。在滤波器的设计过程中，滤波器设计的参数选择很重要，不同的参数对滤波的效果产生很大的影响，因此在进行滤波器设计的过程中，必须花心思去考虑好设计参数，以便得到更好的滤波器。方才采用了4种不同的滤波器对信号进行了滤波处理，以下选取Butterworth滤波前后图像进行对比放大：左肺信号：右肺信号：呼吸信号：从滤波前后的图像对比可以看出，信号的噪声已经得到较好的滤除，图像清晰，便于进行以后所需要的信号数据分析的工作。六、滤波器分析方才使用了IIR数字滤波器中三种最常用的滤波器类型，共4种滤波器，它们在滤波阶数设置一样的情况下对采集的信号做出滤波，现对比4种滤波器的滤波效果，分析它们各自的特点：四种滤波器滤波后的图像：四种滤波器的幅频特性曲线：由以上结果图分析可知，滤波器取相同阶数时：Butterworth低通滤波器具有单调下降的幅频特性，但过渡带最宽。两种类型的Chebyshev滤波器的过渡带宽度相等，比Butterworth滤波器的过渡带窄，但比椭圆滤波器的过渡带宽。ChebyshevI型滤波器在通带具有等波纹幅频特性，阻带内单调；ChebyshevII型滤波器通带内单调，二阻带内等波纹。椭圆滤波器阻带和通带都是等波纹的，但过渡带最窄。七、个人思考Matlab软件是一款功能很强大的软件，利用Matlab软件可以方便地进行滤波器的设计。课程设计的主要目的在于巩固所学的数字信号处理理论知识，因此如果仅仅是利用Matlab简单地进行滤波器设计并处理数据，这难以达到课程设计的主要目的，因此，为了更好地去理解各种滤波器，也锻炼自己的编写程序的能力，下面通过书本上的理论知识加上Matlab函数对FIR、IIR中的集中常见滤波器进行设计，以加深对滤波器的理解。7.1FIR滤波器的窗口设计法窗函数名过渡带宽最小阻带衰减近似值准确值矩形4π1.8π21dBBartlett/巴特利特86.1π25dBHanning/汉宁86.2π44dBHamming/哈明86.6π53dBBlackman/布莱克曼1211π74dBMatlab提供了命令用于产生这些窗函数：（1）w=boxcar(M)在数组w中产生M点的矩形窗函数。（2）w=triang(M)在数组w中产生M点的Bartlett（三角形）窗函数。（3）w=hann(M)在数组w中产生M点的Hanning窗函数。（4）w=hamming(M)在数组w中产生M点的Hamming窗函数。（5）w=blackman(M)在数组w中产生M点的Blackman窗函数。（6）w=kaiser(M，beta)在数组w中产生M点的kaiser窗函数。利用这些指令便可以用Matlab设计基于窗函数法的FIR滤波器。为了方便使用，在这里给出一个理想低通脉冲响应hd(n)，编写程序如下：functionhd=ideal_lp(wc,M);%wc为截止频率，M为窗函数长度alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha;fc=wc/pi;hd=fc*sinc(fc*m);end例子：设计一个数字FIR低通滤波器，器技术指标如下：wp=0.2π，ws=0.3π，Rp=0.25dB，As=50dB题解: 从表格可知，用Hamming和Blackman窗函数都能提供50dB的阻带衰减，现选用Hamming窗函数，因为它能给出比较小的过渡带，因此有较低的阶。尽管在Hamming窗函数滤波器的设计中没有用到通带波纹值Rp，但是还是必须从设计中校核真正的波纹值，以证实波纹确实符合题目要求。设计结果如下：Matlab代码运行结果：>>M=67>>Rp=0.0394>>As=52结果显示真实的阻带衰减是52dB，通带波纹是0.0394dB，显然满足题目设计要求。我们知道，Kaiser窗函数对同一个M值可以给出不同的过渡带宽，它的这种性能使其可以和Hamming窗相匹敌，下面根据上述指标，选择Kaiser窗设计这个低通滤波器，设计结果如下：运行结果：>>M=61>>beta=4.5513>>Rp=0.0442>>As=527.2IIR数字滤波器设计IIR数字滤波器具有无限长脉冲响应，因此能偶与模拟滤波器相匹敌；一般来说，所有的模拟滤波器都有无限长脉冲响应。因此，IIR数字滤波器设计的基本方法是利用复值映射将模拟滤波器变换为数字滤波器。复值映射的变换方法有很多种，如脉冲响应不变变换、双线性变换、匹配Z变换等，在实际中采用的最为普遍的方法是双线性变换法。在Matlab中，有buttord、cheb2ord、cheb2ord和ellipord函数可以在已知设计指标时给出最小滤波阶数N和截止频率wn；以下程序尝试采用定义中给出的公式编程确定这些参数，以加深对双线性变换的理解。下面利用双线性变换法对三种最广泛使用的滤波器进行设计。例子：设计一个数字FIR低通滤波器，器技术指标如下：wp=0.2π，ws=0.3π，Rp=1dB，As=15dB（1）巴特沃斯低通滤波器：设计结果：>>ButterworthFilterOrder=6>>直接型>>b=0.00060.00350.00870.01160.00870.00350.0006>>a=1.0000-3.31434.9501-4.14332.0275-0.54580.0628>>级联型>>BA= 1.00002.00001.00001.0000-0.94590.2342 1.00002.00481.00481.0000-1.05410.3753 1.00001.99520.99521.0000-1.31430.7149>>b0= 5.7969e-04则这个滤波器系统函数是：H（2）数字切比雪夫I型、II型低通滤波器：运行结果： %切比雪夫I型>>ChebyshevIFilterOrder=4>>直接型>>b=0.00180.00730.01100.00730.0018>>a=1.0000-3.05433.8290-2.29250.5507>>级联型>>BA= 1.00002.00001.00001.0000-1.55480.6493 1.00002.00001.00001.0000-1.49960.8482>>b0=0.0018则这个滤波器系统函数是：H%切比雪夫II型>>ChebyshevIIFilterOrder=4>>直接型>>b=0.1797-0.09160.2525-0.09160.1797>>a=1.0000-1.55081.3423-0.47070.1079>>级联型>>BA= 1.00000.55741.00001.0000-0.41830.1503 1.0000-1.06711.00001.0000-1.13250.7183>>b0=0.1797则这个滤波器系统函数是：H（3）数字椭圆低通滤波器：运行结果：>>EllipticFilterOrder=3>>直接型>>b=0.1214-0.0511-0.05110.1214>>a=1.0000-2.11121.7843-0.5325>>级联型>>BA= 1.00001.000001.0000-0.618301.0000-1.42111.00001.0000-1.49280.8612>>b0= 0.1214则这个滤波器系统函数是：H三类滤波器的比较：原型阶次N最小阻带衰减/dB巴特沃斯615切比雪夫425椭圆327在给定了的相同的滤波器设计参数后，从阶数上看，椭圆滤波器给出了最低阶数的设计而巴特沃斯滤波器给出了三者中最高的阶数设计，但是我们也知道，椭圆滤波器在通带内非线性相位最为严重，因此需要根据实际情况取舍，选择合适的滤波器设计。八、拓展除了采用编程的方法实现滤波器之外，Matlab中自带工具箱FDATool也能很方便快速的实现滤波器的设计。另外FPGA、DSP等实现数字滤波器算法时，经常要用到滤波器系数，采用FDATool工具箱也能快速的得到滤波器系数。方法1：在Matlab主界面输入FDATool，回车即可调用FDATool工具箱。方法2：单击Matlab主窗口下方的“Start”按钮，如图B.1所示，选择菜单“ToolBox”

→“Filter

Design”

→“Filter

Design

&

Analysis

Tool（FDATool）”命令，打开FDATool。其界面如下所示：8.1FDATool界面说明FDATool（FilterDesign&AnalysisTool）是Matlab信号处理工具箱里专用的滤波器设计分析工具，Matlab6.0以上的版本还专门增加了滤波器设计工具箱（FilterDesignToolbox）。FDATool可以设计几乎所有的基本的常规滤波器，包括fir和iir的各种设计方法。它操作简单，方便灵活。FDATool界面总共分两大部分，一部分是DesignFilter，在界面的下半部，用来设置滤波器的设计参数，另一部分则是特性区，在界面的上半部分，用来显示滤波器的各种特性。DesignFilter部分主要分为：FilterType（滤波器类型）选项，包括lowpass（低通）、highpass（高通）、bandpass（带通）、bandstop（带阻）和特殊的fir滤波器。DesignMethod（设计方法）选项，包括IIR数字滤波器的butterworth（巴特沃思）法、chebyshevtypei（切比雪夫i型）法、chebyshevtypeii（切比雪夫ii型）法、elliptic（椭圆滤波器）法和fir滤波器的equiripple法、least-squares（最小乘方）法、window（窗函数）法。FilterOrder（滤波器阶数）选项，定义滤波器的阶数，包括specifyorder（指定阶数）和minimumorder（最小阶数）。在specifyorder中填入所要设计的滤波器的阶数（n阶滤波器，specifyorder＝n-1），如果选择minimumorder则Matlab根据所选择的滤波器类型自动使用最小阶数。FrenquencySpecifications选项，可以详细定义频带的各参数，包括采样频率fs和频带的截止频率。它的具体选项由filtertype选项和designmethod选项决定，例如bandpass（带通）滤波器需要定义fstop1（下阻带截止频率）、fpass1（通带下限截止频率）、fpass2（通带上限截止频率）、fstop2（上阻带截止频率），而lowpass（低通）滤波器只需要定义fstop1、fpass1。采用窗函数设计滤波器时，由于过渡带是由窗函数的类型和阶数所决定的，所以只需要定义通带截止频率，而不必定义阻带参数。MagnitudeSpecifications选项，可以定义幅值衰减的情况。例如设计带通滤波器时，可以定义wstop1（频率fstop1处的幅值衰减）、wpass（通带范围内的幅值衰减）、wstop2（频率fstop2处的幅值衰减）。当采用窗函数设计时，通带截止频率处的幅值衰减固定为6db，所以不必定义。WindowSpecifications选项，当选取采用窗函数设计时，该选项可定义，它包含了各种窗函数。8.2FDATool性能分析介绍在设计完所需的滤波器后，往往需要对设计好的FIR滤波器进行相关的性能分析，以便了解该滤波器是否满足设计要求。分析操作步骤如下：（1）选择FDATool的菜单“Analysis”→“Magnitude

Response”，启动幅频响应，x轴为频率，y轴为幅度值（单位为dB）。（2）选择菜单“Analysis”→“Phase

Response”，启动相频响应分析；（3）选择菜单“Analysis”→“Magnitude

and

Phase

Response”来启动滤波器幅频特性与相频特性的比较分析；（4）选择菜单“Analysis”→“Group

Delay

Response”，启动群时延分析。FDATool还提供了以下几种分析工具：群时延响应分析（GroupDelayResponse）。冲激响应分析（Impulse

Response）；阶跃响应分析（Step

Response）；零极点图分析（Pole/Zero

Plot）。8.3FDATool滤波器系数导出为导出设计好的滤波器系数，选择FDATool菜单的“File”→“Export”命令，打开Export（导出）对话框，如下图所示。在该窗口中，选择导出到工作区（Workplace）。这时滤波器系数就存入到一个一维变量Num中了。这时Num中的元素是以小数形式出现的。九、收获与感想滤波器的设计在实际生活中必不可少，当然也比课程设计的题目要复杂得多，但是它们设计的成功与否，归结到底还是参数的设置以及滤波器类型的选择。通过前面的滤波器设计，我们已经知道每一种滤波器都有各自的优缺点，在实际情况下到底该选取哪一种滤波器，主要取决于每种滤波器的优点在设计中的重要性。因此，滤波器的设计必须做到以下几点：确定所需的滤波器的用途，并对其进行性能需求分析，以明确该系统所需的各种性能指标；必须有系统的设计思维，把每一个细节都放到整个系统中考虑，考虑整个系统的可行性、完整性、稳定性和功能的实现；经过这次的课程设计，我对《数字信号处理》这门课程有了更进一步的认识，对所学的知识也有了更深一层的理解和掌握，不再仅仅是停留在“学过”这个层面上，还学会了如何在实际应用中运用自己所学到的知识。或许在平时的课堂上，学习理论知识的时候，我们会觉得“这种东西”有什么用啊，因而没有过多地去在意它，当进行这个课程设计的时候，次贸然发现其实那些理论是多么的重要，而且在实际运用中并不像理论课上那样的枯燥乏味，更多的是需要灵活地去运用所学的理论解决实际问题，不再像写作业那样是几乎的过程那样几乎是定势的一种计算。在设计的过程中遇到过不少的困难，然后需要查阅相关的资料并学习相关的知识一解决问题，这是个考验耐心的过程，在这个过程中必须保持良好的心态，切记心浮气躁，不然问题只会无法有效地解决，造成恶性循环。面对一个问题，或许思索良久都没有好的解决办法，这是不妨放开思维，大胆地尝试，说不定通过另一种思路问题就迎刃而解了，有时候只要稍微转一下弯，不要钻牛角尖，结果就会不一样了。总而言之，通过这次的课程设计，不仅让我巩固了课堂上学习到的专业知识，也锻炼了我各方面的能力，让我收获了很多。十、MATLAB程序代码5.1画出三组原始数据的时域图像：Matlab代码：clearA=xlsread('data.xls');a=A(:,1)';b=A(:,2)';c=A(:,3)';subplot(3,1,1);plot(a);xlim([100,1000]);title('OriginalData_a');subplot(3,1,2);plot(b);xlim([100,1000]);title('OriginalData_b');subplot(3,1,3);plot(c);xlim([100,1000]);title('OriginalData_c');5.2画出三组原始数据的频谱图：Matlab代码：fs=7.812;N=1000;%采样频率和数据点数n=0:N-1;t=n/fs;%时间序列ya=fft(a,N);%对信号进行快速Fourier变换mag1=abs(ya);%求得快速Fourier变换后的振幅yb=fft(b,N);mag2=abs(yb);yc=fft(c,N);mag3=abs(yc);f=n*fs/N;%频率序列subplot(3,1,1);plot(f,mag1);%绘出随频率变化的振幅xlim([0,7.812/2]);ylim([0,10]);title('FFTofdata_a');subplot(3,1,2);plot(f,mag2);xlim([0,7.812/2]);ylim([0,10]);title('FFTofdata_b');subplot(3,1,3);plot(f,mag3);xlim([0,7.812/2]);ylim([0,10]);title('FFTofdata_c');xlabel('frequency/Hz');滤波器的设计与信号滤波：5.3.1IIR数字滤波器设计： 5.3.1.1．巴特沃斯低通滤波器：Matlab代码：fs=7.812;Wn=2*0.5/fs;%归一化N=20;%设置阶数[B,A]=butter(N,Wn,'low');X=filter(B,A,a);%a可以换成b、c以对b、c组数据滤波[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用butterworth滤波后图像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅频特性曲线');xlabel('单位为π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相频特性曲线');xlabel('单位为π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅频特性曲线(dB)');xlabel('单位为π');5.3.1.2．切比雪夫I型低通滤波器：Matlab代码：fs=7.812;wn=2*0.5/fs;N=10;%归一化，并设定阶数Rp=0.025;%通带波纹[B,A]=cheby1(N,Rp,wn);X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用切比雪夫I型滤波后图像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅频特性曲线');xlabel('单位为π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相频特性曲线');xlabel('单位为π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅频特性曲线(dB)');xlabel('单位为π');5.3.1.3．切比雪夫II型低通滤波器：Matlab代码：fs=7.812;wn=2*0.5/fs;N=10;As=60;%阻带衰减[B,A]=cheby2(N,As,wn);X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用切比雪夫II型滤波后图像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅频特性曲线');xlabel('单位为π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相频特性曲线');xlabel('单位为π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅频特性曲线(dB)');xlabel('单位为π');5.3.1.4．椭圆低通滤波器：Matlab代码：fs=7.812;wn=2*0.5/fs;Rp=0.025;As=60;N=5;%设定设计参数[B,A]=ellip(N,Rp,As,wn);X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用ellip滤波后图像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅频特性曲线');xlabel('单位为π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相频特性曲线');xlabel('单位为π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅频特性曲线(dB)');xlabel('单位为π');5.3.2．FIR（窗函数法）滤波器设计：5.3.2.1．Hamming窗函数滤波器：Matlab代码：wp=0.4;ws=0.6;Rp=0.2;As=60;fs=7.812;delta=ws-wp;N=ceil(6.6*pi/delta);wn=(wp+ws)/2;[B,A]=fir1(N,wn/pi,hamming(N+1));%选择窗函数X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用hamming窗函数滤波后图像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅频特性曲线');xlabel('单位为π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相频特性曲线');xlabel('单位为π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅频特性曲线(dB)');xlabel('单位为π');7.1Hamming窗函数滤波器设计：Matlab代码：clearwp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)+1n=[0:1:(M-1)];wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_ham=(hamming(M))';h=hd.*w_ham;[H,w]=freqz(h,[1],1000,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));delta_w=2*pi/1000;Rp=-(min(db(1:1:(wp/delta_w+1))))As=-round(max(db(ws/delta_w+1:1:501)))%%%%%%%%%%plots%%%%%%%%%%subplot(2,2,1);stem(n,hd);title('IdealImpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h_d(n)');subplot(2,2,2);stem(n,w_ham);title('HammingWindow');axis([0M-101.1]);xlabel('n');ylabel('w(n)');subplot(2,2,3);stem(n,h);title('ActualimpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h(n)');subplot(2,2,4);plot(w/pi,db,w/pi,-50,'-k');title('MagnitudeResponseindB');axis([01-12010]);xlabel('frequencyinpiunits');ylabel('dB');7.1Kaiser窗函数滤波器设计：Matlab代码：clearwp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;As=50;M=ceil((As-7.95)/(2.285*tr_width)+1)+1n=[0:1:(M-1)];beta=0.1102*(As-8.7)wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_kai=(kaiser(M,beta))';h=hd.*w_kai;[H,w]=freqz(h,[1],1000,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));delta_w=2*pi/1000;Rp=-(min(db(1:1:(wp/delta_w+1))))As=-round(max(db(ws/delta_w+1:1:501)))%%%%%%%%%%plots%%%%%%%%%%subplot(2,2,1);stem(n,hd);title('IdealImpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h_d(n)');subplot(2,2,2);stem(n,w_kai);title('kaiserWindow');axis([0M-101.1]);xlabel('n');ylabel('w(n)');subplot(2,2,3);stem(n,h);title('ActualimpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h(n)');subplot(2,2,4);plot(w/pi,db,w/pi,-50,'-k');title('MagnitudeResponseindB');axis([01-12010]);xlabel('frequency