吹管音乐滤波去噪使用TUKEYWIN设计的FIR滤波器要点

《吹管音乐滤波去噪使用TUKEYWIN设计的FIR滤波器要点》由会员分享，可在线阅读，更多相关《吹管音乐滤波去噪使用TUKEYWIN设计的FIR滤波器要点（18页珍藏版）》请在装配图网上搜索。

1、吹管音乐滤波去噪一一使用TUKEYWIN计的FIR滤波器第3页共17页吹管音乐滤波去噪使用TUKEYWIN 设计的FIR滤波器学生姓名： 指导老师：摘 要 本课程设计主要内容是设计利用窗口设计法选择TUKEYWIN窗设计一个FIR滤波器，对一段含噪吹管音乐信号进行滤波去噪处理并根据滤波前后的波形和频谱分析 滤波性能。本课程设计仿真平台为 MATLAB7.0 ,开发工具是M语言编程。首先下载一 段吹管音乐信号，并人为加入一单频噪声，然后对信号进行频谱分析以确定所加噪声频 率，并设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱并进行分析。由 分析结果可知，滤波器后的语音信号与原始信号基本一

2、致，即设计的FIR滤波器能够去除信号中所加单频噪声，达到了设计目的。关键词 滤波去噪；FIR滤波器；TUKEYWIN窗;MATLAB1引言本课程设计主要解决在含噪情况下对吹管音乐信号的滤波去噪处理，处理时采用的是利用窗口设计法选择Tukeywin窗设计的FIR滤波器。通过课程设计了解 FIR滤波器 设计的原理和步骤，掌握用 Matlab语言设计滤波器的方法及编程。通过观察音乐信号 滤波前后的时域波形的比较，加深对滤波器作用的理解。1.1 课程设计目的熟悉Matlab语言环境，掌握 Matlab语言的编程规则，利用 TUKEYWIN窗函数设 计法来设计符合要求的FIR滤波器来实现音乐信号的滤波去

3、噪。 并绘制滤波前后的时域 波形和频谱图。根据图形分析判断滤波器设计的正确性。通过本次课程设计熟悉利用 TUKEYWIN窗函数法设计FIR滤波器的过程。增强自己独立解决问题的能力，提高自 己的动手能力。加深对理论知识联系实际问题的理解。为以后的工作奠定坚实的基础。1.2 课程设计要求(1)滤波器指标必须符合工程实际。(2)设计完后应检查其频率响应曲线是否满足指标。(3)处理结果和分析结论应该一致，而且应符合理论(4)独立完成课程设计并按要求编写课程设计报告书1.3 课程设计平台MATLAB软件包括五大通用功能：数值计算功能( Nemeric);符号运算功能 (Symbolic);数据可视化功能

4、(Graphic);数据图形文字统一处理功能(Notebook) 和建模仿真可视化功能(Simulink)。其中，符号运算功能的实现是通过请求MAPLE内核计算并将结果返回到 MATLAB命令窗口。该软件有三大特点：一是功能强大；二 是界面友善、语言自然；三是开放性强。目前，Mathworks公司已推出30多个应用工具箱。MATLAB在线性代数、矩阵分 析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处 理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众 多领域的理论研究和工程设计中得到了广泛应用。MATLAB在信号与系统中的应用主要包括符号

5、运算和数值计算仿真分析。由于信 号与系统课程的许多内容都是基于公式演算，而MATLAB借助符号数学工具箱 提供的符号运算功能能基本满足信号与系统课程的需求。例如，解微分方程、傅里叶正反变换、 拉普拉斯正反变换、z正反变换等。MATLAB在信号与系统中的另一主要应用是数值计 算与仿真分析，主要包括函数波形绘制、函数运算、冲激响应与阶跃响应仿真分析、信 号的时域分析、信号的频谱分析、系统的 S域分析、零极点图绘制等内容。数值计算仿 真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用MATLAB进行信号处理领域的各种分析和实际应用打下基础1 o2设计原理窗函数设计法的基本思想是用 FIR

6、DF逼近希望的滤波特性。设希望逼近的滤波器 的频率响应函数为Hd(ejw),其单位脉冲响应用hd(n)表示。为了设计简单方便，通常选 择Hd(ejw)为具有片段常数特性的理想滤波器。因此 hd(n)是无限长非因果序列，不能直 接作为FIRDF的单位脉冲响应。窗函数设计法就是截取hd(n)为有限长的一段因果序列, 并用合适的窗函数进行加权作为 FIRDF的单位脉冲响应hd(n) 2o2.1 FIR滤波器FIR滤波器是有限长单位脉冲响应数字滤波器的简称，FIR滤波器的最大优点是可以实现线性相位滤波。FIR滤波器满足了在数字通信和图像传输与处理等应用场合对线 性相位的要求，FIR滤波器是全零点滤波器

7、，硬件和软件实现结构简单，不用考虑稳定 性问题。当幅频特性指标相同(不考虑相位特性)时，FIR滤波器的阶数比IIRDF的阶数高得多。所以在要求线性相位滤波的应用场合，一般都用FIR滤波器。FIR数字滤波器设计的基本步骤如下：(1)确定指标：在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器 的技术指标。在很多实际应用中，数字滤波器常常被用来实现选频操作。因此，指标的 形式一般在频域中给出幅度和相位响应。幅度指标主要以两种方式给出。第一种是绝对 指标。它提供对幅度响应函数的要求，一般应用于 FIR滤波器的设计。第二种指标是相 对指标。它以分贝值的形式给出要求。(2)逼近：确定了技术指标后

8、，就可以建立一个目标的数字滤波器模型。通常用理 想的数字滤波器模型。之后，利用数字滤波器的设计方法，设计出一个实际滤波器模型 来逼近给定的目标。(3)性能分析和计算机仿真：上两步的结果是得到以差分或冲激响应描述的滤波器。 根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求； 或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。FIR滤波器的设计问题实际上是确定能满足所要求的转移序列或脉冲响应的常数的 问题，设计方法主要有窗函数、频率响应法和等波纹最佳逼近法等。2.2 窗口设计法吹管音乐滤波去噪一一使用TUKEYWIN计的FIR滤波器第5页共17页窗口设计法是一种通

9、过截断和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的 技术指标。在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形 式一般为在频域中以分贝值给出的相对幅度响应和相位响应。窗口设计法步骤如下：(1)根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度 No窗函数 的类型可根据最小阻带衰减AS独立选择。(2)根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)。(3)由性能指标确定窗函数 W(n)和长度No(4)求得实际滤波器的单位脉冲响应 h(n), h(n)即为所设计FIR滤波器系数向量b(

10、n)0h(n) = hd (n).W(n)(2.1)常见的窗函数性能表如下图2.1所示:图2.1常见窗函数性能表名称滤波器最小阻带 衰减名称滤波器最小阻带衰 减过渡市见过渡市见矩形1.8 兀 /M21dBPARZENWIN6.6 兀 /M56dB巴特利特6.1 兀/M25dBFLATTOPWIN19.6 兀/M108dB汉宁6.2 兀 /M44dBGAUSSWIN5.8 兀 /M60dB汉明6.6 兀 /M51dBBARTHANNWIN3.6 兀 /M40dB11 九/M74dBBLACKMANHARRIS16.1 兀/M109dB布莱克曼BOHMANWIN5.8 兀 /M51.5dBCHEB

11、WIN15.2 兀 /M1113dBNUTTALLWIN15.4 兀 /M108dBTUKEYWIN2.4 兀 /M22dBw =tukeywin(M)在数组w中广生M点的Tukeywin窗函数。由上面查表可知：Tukeywin窗的滤波器的过渡带宽为2.4冗/M最小阻带衰减22db。利用wvtool可绘制出Tukeywin窗函数波形如下：N=128;t=tukeywin(N);wvtool(t);Ojg一Window viewerTim e domain20406080100120SamplesLeakage Factor: 3.56 %.64.2U IO.IO.QFrequency doma

12、in 40 200*120 002 g 0.60.8Normalized Frequency (xnrftditeample)崛inlobe wicfth0.017576Reiativs sidelobe sttenuetion; -15.1 d6图2-2 Tukeywin窗函数波形3设计步骤3.1 设计流程图根据设计的要求，首先从网上下载一段音乐信号，对音乐信号加入噪声干扰，再利 用Tukeywin窗设计合理的FIR滤波器。最后用滤波器对干扰后的音乐信号进行滤波去 噪。具体设计流程图如下图3.1所示：吹管音乐滤波去噪一一使用TUKEYWIN计的FIR滤波器第7页 共17页图3-1音乐信号滤波

13、去噪流程图3.2 录制语音信号在网上下载一段mp3音乐，再利用视频转换器将其转换成.wav格式，且为单声道，再将此.wav格式音乐控制在3秒内，以减少设计中的误差。然后在Matlab软件平台下, 利用函数wavread对语音信号进行采样，记住采样频率和采样点数，采集完成后在信号 中加入一个单频噪声，绘制原音乐信号和加噪后的音乐信号的时域和频域的波形图。具体调用如下:数字信号meilideshenhua.wav); %输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits 是波形文件每样本的编码位数。sound(x,fs,bits); %按指定的

14、采样率和每样本编码位数回放N=length(x);%计算信号x的长度fn=2100;%单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; %计算时间范围，样本数除以采样频率x=x(:,1); y=x+0.1*sin(fn*2*pi*t);sound(y,fs,bits); %应该可以明显听出有尖锐的单频啸叫声X=abs(fft(x);Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换，取幅度谱X=X(1:N/2); Y=Y(1:N/2);% 截取前半部分deltaf=fs/N; %计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; %计算频谱频率范围得到

15、的波形图如图3.2所示：时间(t)频率时间频率图3-2原始音乐信号与加噪后的音乐信号时域、频域波形图吹管音乐滤波去噪一一使用TUKEYWIN计的FIR滤波器第9页共17页3.3 滤波器设计在此用窗函数法设计一个带阻的FIR滤波器fD1 =2000Hz p Ifs2 =2150Hz此处给出各滤波器的性能指标：%=2050Hz、s1, As= 20dB, Ap=1 dB,其截止频率在单频噪声干扰f p2 = 2200Hz附近。再进行加窗设计，利用 Matlab中的函数freqz画出各滤波器的频率响应。 具体调用如下：fpd=2000;fsd=2050;fsu=2150;fpu=2200;Rp=1;

16、As=20;% 带阻滤波器设计指标fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心 频率，和频率间隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将 Hz 为单位的模拟频率换 算为rad为单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(2.4*pi/dw)+1;%计算TUKEYWIN窗设计该滤波器时需要的阶数n=0:M-1;%定义时间范围w_tuk=tukeywin(M); % 产生 M 阶的 TUKEYWIN

17、 窗hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想 带阻滤波器的脉冲响应h_bs=w_tuk*hd_bs; %用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); %调用自编函数计算滤波器的频率特性得到的波形图如图3.3所示：滤波器相位响应图滤波器脉冲响应图图3-3 滤波器波形图注：滤波器的设计指标o(s=20dB,图中可以明显地看到利用TUKEYWIN窗函数所 设计的滤波器基本上可以达到滤波效果。3.4 信号滤波处理用自己设计的各滤波器分别对采集的信号进行滤波，在

18、Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。在同一张大图 里分别绘制原始信号x,加噪信号x,滤波去噪信号y的时域波形和频谱，以便比较和 分析。滤波程序如下：y_fil=fffilt(h_bs,y); % 用设计好的滤波器对y进行滤波Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半得到的波形图如图3.4所示：吹管音乐滤波去噪一一使用TUKEYWI限计的FIR滤波器第18页共17页原始音乐信号幅度诺X原始音乐信号时间X00511 522.6时间t加干扰音乐信号时间对00.511.52

19、25时间t滤波后的音乐信号时间y00511 522.5时间t200010LIO00 600 1000 1600 2000 2600频率:f加干扰后的音乐信号幅度谱X10500 1000 1500 2000 2500频率f滤波后的音乐信号幅度谱Y2000100002000100000 500 1000 1500 2000 2500 频率f图3.4滤波前后音乐信号的时域波形图和频谱图对比图注：从图中可以看出滤波后的信号的幅度有所减小，说明加噪信号通过滤波器之后 不可能毫无变化，也就是说参数的设置以及窗函数的选择都会影响滤波后信号的质量。3.5 结果分析在Matlab中，函数sound可以对经过tu

20、keywin窗设计的FIR滤波器之后的音乐信号 进行回放。其调用格式：sound (y_fil,fs,bits);可以感觉滤波前后的声音有变化，那尖锐 的噪声明显减小。但是所设计的滤波器并不是理想的滤波器，并且由于tukeywin窗所设计的滤波器本身的性能局限，所以在回放滤波后的音乐的时候还可以听到一些噪声。但从滤波器的波 形图中也可以明显地看到滤波器的设计满足所定义的指标，从另一个方面也可以断言滤波器的设计成功。与相同信源不同滤波方法(Nuttallwin窗)的同学的滤波器相对比，可以看出由于 其衰减偏小，所以性能也较差。与不相同信源（拉弦音乐）相同滤波方法的同学的滤波器相对比，可以看出其滤

21、波 效果基本相近。4由现的问题及解决方法1、在录音时，没有将录音的属性改为 wav的格式，当在Matlab软件平台上调用时， 出现无法调用的提示，最后通过与同学讨论发现，不能直接更改文件后缀，应用转换器 将其转换为wav格式。2、在开始录制音乐信号并将其导入 MATLAB中的时候出现过错误，原因是我所下载的 音乐信号是双音频信号，不符合要求，在老师的指导下我将音乐信号变成了单声道信号， 再次导入的时候MATLAB不再报错。3、在最后对加噪音乐信号进行滤波并绘图的时也出现了错误，绘制出最后的滤波信 号的频谱图时纵坐标出现了负值，但是原信号的频谱图中却没有，回头检查了一遍程序 想到幅度值只有正值，

22、想到是不是忘记加绝对值符号了，在对滤波的信号进行fft变换之后未取绝对值，改正之后再运行程序就不会再出错。4、在滤波处理时，在加入单频噪声部分噪声频率设置为 2200,但滤波器参数设置部分 最高带阻为2200,导致滤波处理后不能滤除噪声信号，而后来将其改为 2100,才达到 了效果。5结束语本次的课程设计，我的任务是利用 Tukeywin窗函数设计FIR滤波器对音乐信号滤 波去噪。在本次课程设计之前，我对 Tukeywin窗函数完全没有了解，因此在看到这个 题目时，我是一头雾水。但是通过自己翻阅资料和询问同学，我掌握了用Tukeywin窗函数设计FIR滤波器的方法步骤，了解了窗函数的基本设计流

23、程。经过这两周的课程设计经过为期两周的课程设计，这次课程设计比较的简单，根 据老师给出的模板，再结合以前学过的知识，并在同学和老师的帮助下我顺利的完 成了任务。不同于在教室里上的理论，这次的课程设计需要将我们平时所学习的知识运 用到实践之中，将知识学以致用。因为是以所学理论为基础，所以在课程设计的过程中， 我又重温滤波器、窗函数等知识，更加熟悉了MATLAB的操作。虽说总体上是比较的简单，但是也会出现一些错误，分析并总结这些错误主要是语 句上的一些用法出现错误，以及没有及时地掌握老师所提醒的部分。但是最后还是在老 师的帮助下，解决了这些问题。课程设计是我们运用所学知识，动手实践的一个很好的机会

24、。它既可以帮助我们加 深对所学知识的理解，又能提高我们运用知识，联系实际，动手实践的能力。而且在设 计过程中可能用到我们没学过的知识，需要我们去查阅资料获取相关信息，这又提高了 我们查找信息和学习新知识的能力。在实物的调试与检测过程中，又会遇到许多意想不 到的问题，需要我们去分析原因和解决问题。也体会到真正的去独立地完成一件事情是 很困难的，同学以及老师的帮助和提醒是必须的。通过这次课程设计，我拓宽了知识面，锻炼了实际操作能力，综合素质也得到了提 高，进一步加深了 了我们对专业的认识和激发了我们对专业的兴趣。虽然课程设计结束了，但是我们的学习还没结束，对知识的进一步学习还需要继续，很开心成功地

25、完成了 这次设计。参考文献1维纳.K英格尔，约翰.G普罗克斯.数字信号处理（MATLAB版）.西安：西安交通大 学出版社，20072程佩青.数字信号处理教程.北京：清华大学出版社，20023张圣勤.MATLAB7.0实用教程.北京：机械工业出版社，20084张志涌.精通MATLAB 6.5版M .北京：北京航空航天大学出版社，2003.5谢德芳.数字信号处理M.北京：科学出版社，2005.6张小虹.信号系统与数字信号处理M.第1版.西安：西安电子科技出版社，2002.附录1:吹管音乐滤波去噪一一使用TUKEYWIN 窗设计的FIR滤波器程序一：%放出原音频数字信号 meilideshenhua

26、.wav);%输入参数为文件的全路径和文件名，%输出的第一个参数是每个样本的值，%fs是生成该波形文件时的采样率，%bits是波形文件每样本的编码位数。sound(x,fs,bits); %按指定的采样率和每样本编码位数回放程序二：%原音乐信号和加噪后的音乐信号的时域和频域的波形图绘制N=length(x);%计算信号x的长度fn=2100;%单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; %计算时间范围，样本数除以采样频率x=x(:,1); y=x+0.1*sin(fn*2*pi*t);sound(y,fs,bits); %应该可以明显听出有尖锐的单频啸叫声X=abs(fft(

27、x);Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换，取幅度谱X=X(1:N/2); Y=Y(1:N/2);% 截取前半部分deltaf=fs/N; %计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; %计算频谱频率范围subplot(2,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度);title(原始音乐信号);axis(0,2.5,-1.5,1.5);grid on;subplot(2,2,2);plot(f,X);xlabel(频率(f);ylabel(幅度谱);title(原始音乐信号幅度谱);axis(0,3000,0

28、,3000);grid on;subplot(2,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title(加干扰后的音乐信号);axis(0,2.5,-1.5,1.5);grid on;subplot(2,2,4);plot(f,Y);xlabel(频率(f);ylabel(幅度谱);title(加干扰后的音乐信号幅度谱);axis(0,3000,0,3000);grid on;程序三：%带阻滤波器设计fpd=2000;fsd=2050;fsu=2150;fpu=2200;Rp=1;As=20;% 带阻滤波器设计指标fcd=(fpd+fsd)/2;fcu=(f

29、pu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率, 和频率间隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;%将Hz为单位的模拟频率换算为rad为单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(2.4*pi/dw)+1;%计算Tukeywin窗设计该滤波器时需要的阶数n=0:M-1;%定义时间范围w_tuk=tukeywin(M); % 产生 M 阶的 Tukeywin 窗hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal

30、_lp(wcu,M);%调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_tuk.*hd_bs; %用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); %调用自编函数计算滤波器的频率特性subplot(2,2,1);plot(w/pi,db);title(FIR滤波器的幅度响应图);xlabel(w/pi);ylabel(db);line(0Q5,-1,-1,Color,r,LineWidth,2,LineStyle,-);line(0,0.5,-20,-20,Color,r,LineWidth,2,LineStyle,-);line(wsd/

31、pi,wsd/pi,-30,15,Color,r,LineWidth,2,LineStyle,-);line(wsu/pi,wsu/pi,-30,15,Color,r,LineWidth,2,LineStyle,-);axis(0.05 0.15 -30 15);grid on;subplot(2,2,2);plot(w/pi,mag);title(FIR滤波器的幅度响应图);xlabel(w/pi);ylabel(幅度 mag);axis(0 0.3 0 1.5);grid on;subplot(2,2,3);plot(w/pi,pha);title(滤波器相位响应图);xlabel(w/p

32、i);ylabel(相位 pha);axis(0 1 -4 4);grid on;subplot(2,2,4);stem(n,h_bs);title(滤波器脉冲响应图);xlabel(n);ylabel(h(n);axis(0 1000 0 0.7);grid on;程序四：%对信号进行滤波y_fil=fftfilt(h_bs,y); %用设计好的滤波器对y进行滤波 Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半 subplot(3,2,1);plot(t,x);title(原始音乐信号时间x);grid on;xlabel(时间:t);

33、ylabel(幅度);axis(0,2.5,-1.5,1.5);subplot(3,2,2);plot(f,X);title(原始音乐信号幅度谱X);xlabel(频率:f);ylabel(幅度)；axis(0 2500 0 2500);grid on;subplot(3,2,3);plot(t,y);title(加干扰音乐信号时间x1);grid on;xlabel(时间 t);ylabel(幅度)；axis(0,2.5,-1.5,1.5);subplot(3,2,4);plot(f,Y);title(加干扰后的音乐信号幅度谱X1);xlabel(频率 f);ylabel(幅度)；axis(

34、0 2500 0 2500);grid on;subplot(3,2,5);plot(t,y_fil);title(滤波后的音乐信号时间y);grid on;xlabel(时间 t);ylabel(幅度)；axis(0,2.5,-1.5,1.5);subplot(3,2,6);plot(f,Y_fil);grid on;title(滤波后的音乐信号幅度谱 Y);xlabel(频率 f);ylabel(幅度)；axis(0 2500 0 2500);sound (y_fil,fs,bits);程序五：% freqz_m(b,a)函数设计程序function db,mag,pha,grd,w =

35、freqz_m(b,a);% freqz子程序的改进版本% db,mag,pha,grd,w = freqz_m(b,a);% db = 0到pi弧度区间内的相对振幅(db)% mag = 0到pi弧度区间内的绝对振幅% pha = 0到pi弧度区间内的相位响应% grd = 0到pi弧度区间内的群迟延% w = 0到pi弧度区间内的501个频率样本向量% b = Ha(z)的分子多项式系数(对FIR b=h)% a = Ha(z)的分母多项式系数(对FIR: a=1)H,w = freqz(b,a,1000,whole);H = (H(1:1:501); w = (w(1:1:501);ma

36、g = abs(H);db = 20*log10(mag+eps)/max(mag);pha = angle(H);% pha = unwrap(angle(H);grd = grpdelay(b,a,w);% grd = diff(pha);% grd = grd(1) grd;% grd = 0 grd(1:1:500); grd; grd(2:1:501) 0;% grd = median(grd)*500/pi;程序六：% ideal_lp(wc,M)函数设计程序function hd = ideal_lp(wc,M);%理想低通滤波器计算% hd = ideal_lp(wc,M)% hd = 0 to M-1之间的理想脉冲响应% wc =截止频率(弧度)% M =理想滤波器的长度alpha = (M-1)/2;n = 0:1:(M-1);m = n - alpha + eps;hd = sin(wc*m) ./ (pi*m);