弹拨音乐滤波去噪使用三角窗研究设计FIR滤波器.doc

长沙理工大学城南学院
《数字信号处理》课程设计报告

任健

院 系 城南学院 专 业 通信工程
班 级 通信1104班 学 号 201185250426
学生姓名 任 健 指导教师 熊文杰
课程成绩 完成日期 2014年7月4日

课程设计成绩评定

院 系 城南学院专 业 通信工程
班 级 通信1104班 学 号 201185250426
学生姓名 任 健指导教师 熊文杰
完成日期 2014 年 7月 4 日

指导教师对学生在课程设计中的评价
评分项目
优
良
中
及格
不及格
课程设计中的创造性成果

学生掌握课程内容的程度

课程设计完成情况

课程设计动手能力

文字表达

学习态度

规范要求

课程设计论文的质量

指导教师对课程设计的评定意见

综合成绩 指导教师签字 2014年7月 4日

课程设计任务书
城南 学院通信工程专业
课程名称
数字信号处理课程设计
时间
2013～2014学年第二学期18～19周
学生姓名
任健
指导老师
熊文杰
题 目
弹拨音乐滤波去噪——使用三角窗设计的FIR滤波器
主要内容：
从网站上下载一段弹拨音乐，播放直观感觉声音特点，绘制波形并观察其频谱特点。加入一个带外单频噪声，播放感觉声音变化。用三角窗设计满足指标的FIR滤波器，对该和弦信号进行滤波去噪处理，比较滤波前后波形和频谱并进行分析，根据结果和学过的理论得出合理结论。并播放感觉去噪后声音与原音及加噪后的不同。与不同信源相同滤波方法的同学比较各种信源的特点，与相同信源不同滤波方法的同学比较各种滤波方法性能的优劣。
要求：
（1）滤波器指标必须符合工程实际。
（2）设计完后应检查其频率响应曲线是否满足指标。
（3）处理结果和分析结论应该一致，而且应符合理论。
（4）独立完成课程设计并按要求编写课程设计报告书。
应当提交的文件：
（1）课程设计学年论文。
（2）课程设计附件（主要是源程序）。
语音信号滤波去噪——使用三角窗设计FIR滤波器

学生姓名：任健 指导老师：熊文杰
摘 要本课程设计主要是通过使用三角窗设计一个FIR滤波器以对语音信号进行滤波去噪处理。本设计首先用麦克风采集一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用三角窗设计一个满足指标的FIR滤波器，对该语言信号进行滤波去早处理，比较滤波前后的波形和频谱分析，根据结果和学过的理论的出合理的结论。通过对比滤波前后的波形图，深入了解滤波器的相关技术指标和性能，掌握设计滤波器的基本方法，通过程序调试及完善，该设计基本满足设计要求。

关键词MATLAB；三角窗；FIR滤波器；滤波去噪

1 引 言
数字信号处理是将信号以数字方式表示并处理的理论和技术。数字信号处理与模拟信号处理是信号处理的子集。数字信号处理的目的是对真实世界的连续模拟信号进行测量或滤波。因此在进行数字信号处理之前需要将信号从模拟域转换到数字域，这通常通过模数转换器实现。而数字信号处理的输出经常也要变换到模拟域，这是通过数模转换器实现的[1]。
数字滤波器是一种用来过滤时间离散信号的数字系统,它是通过对抽样数据进行数学处理来达到频域滤波的目的。随着现代通信的数字化,数字滤波器变得更加重要。数字滤波器的种类很多，但总的来说可以分成两大类，一类是经典滤波器，另一类可称为现代滤波器。从滤波特性方面考虑，数字滤波器可分成数字高通、数字低通、数字带通和数字带阻等滤波器。从实现方法上考虑，将滤波器分成两种，一种称为无限脉冲响应滤波器，简称IIR（Infinite Impulse Response）滤波器，另一种称为FIR（Finite Impulse Response）滤波器[2]。设计FIR数字滤波器的方法有窗函数法、频率采样法和等波纹最佳逼近法等。

1.1 课程设计目的
用麦克风采集一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用三角窗设计一个满足指标的FIR滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。与不同信源相同滤波方法的同学比较各种信源的特点，与相同信源不同滤波方法的同学比较各种滤波方法性能的优劣。

1.2 课程设计的要求
（1）滤波器指标必须符合工程设计。
（2）设计完后应检查其频率响应曲线是否满足指标。
（3）处理结果和分析结论应该一致，而且应符合理论。
（4）在老师的指导下，独立完成课程设计并按要求编写课程设计报告。

1.3设计平台
MATLAB是矩阵实验室（Ma trix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发，数学可视化，数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。是一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。它集图示和精确计算于一身，在应用数学、物理、化工、机电工程、医药、金融和其他需要进行复杂数值计算的领域得到广泛应用。它不仅是一个在各类工程设计中便于使用的计算工具，而且也是一个在数学、数值分析和工程计算等课程教学中的优秀的教学工具，在世界各地的高等院校中十分流行，在各类工业应用中更有不俗的表现。MATLAB可以在几乎所有的PC机和大型计算机上运行，适用于Windows、UNIX等各种系统平台。

2 设计原理
2.1数字信号处理
数字信号处理的算法需要利用计算机或专用处理设备如数字信号处理器（DSP）和专用集成电路（ASIC）等。数字信号处理技术及设备具有灵活、精确、抗干扰强、设备尺寸小、造价低、速度快等突出优点，这些都是模拟信号处理技术与设备所无法比拟的。
数字信号处理的核心算法是离散傅立叶变换(DFT)，是DFT使信号在数字域和频域都实现了离散化，从而可以用通用计算机处理离散信号。而使数字信号处理从理论走向实用的是快速傅立叶变换(FFT)，FFT的出现大大减少了DFT的运算量，使实时的数字信号处理成为可能、极大促进了该学科的发展。
数字信号处理的特征和分类。信号(signal)是一种物理体现，或是传递信息的函数。而信息是信号的具体内容。
模拟信号(analog signal)：指时间连续、幅度连续的信号。
数字信号(digital signal)：时间和幅度上都是离散(量化)的信号。
数字信号可用一序列的数表示，而每个数又可表示为二制码的形式，适合计算机处理。
一维(1-D)信号：一个自变量的函数。
二维(2-D)信号：两个自变量的函数。
多维(M-D)信号：多个自变量的函数。
信号处理的内容：滤波、变换、检测、谱分析、估计、压缩、识别等一系列的加工处理。
随着信息时代、数字世界的到来，数字信号处理已成为一门极其重要的学科和技术领域。

2.2窗口设计法
窗口设计法的基本思想是用FIR逼近希望的滤波特性。设希望逼近的滤波器的频率响应为，其单位脉冲响应用表示。为了设计简单方便，通常选择为具有片段常数特性的理想滤波器。因此是无限长非因果序列，不能直接作为FIR的单位脉冲响应。窗口设计法就是截取为有限长的一段因果序列，并用合适的窗口函数进行加权作为FIR的单位脉冲响应 。
常见的窗函数，可以分为以下主要类型：
（1） 幂窗--采用时间变量某种幂次的函数，如矩形、三角形、梯形或其它时间（t）的高次幂；
（2）三角函数窗--应用三角函数，即正弦或余弦函数等组合成复合函数，例如汉宁窗、海明窗等；
（3）指数窗--采用指数时间函数，例如高斯窗等。
常见窗函数性能表如图2.2所示：
名称
滤波器
过渡带宽
最小阻带衰减
名称
滤波器
过渡带宽
最小阻带衰减
矩形
1.8π/M
21dB
PARZENWIN
6.6π/M
56db
巴特利特
6.1π/M
25dB
FLATTOPWIN
19.6π/M
108db
汉宁
6.2π/M
44dB
GAUSSWIN
5.8π/M
60db
汉明
6.6π/M
51dB
BARTHANNWIN
3.6π/M
40db
布莱克曼
11π/M
74dB
BLACKMANHARRIS
16.1π/M
109db
BOHMANWIN
5.8π/M
51.5db
CHEBWIN
15.2π/M
113db
NUTTALLWIN
15.4π/M
108db
TUKEYWIN
2.4π/M
22db
图2.1常见窗函数性能表

2.3 FIR滤波器
FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。
FIR滤波器的种类
目前，FIR滤波器的硬件实现有以下几种方式：
（1）数字集成电路FIR滤波器
一种是使用单片通用数字滤波器集成电路，这种电路使用简单，但是由于字长和阶数的规格较少，不易完全满足实际需要。虽然可采用多片扩展来满足要求，但会增加体积和功耗，因而在实际应用中受到限制。
（2）DSP芯片FIR滤波器
另一种是使用DSP芯片。DSP芯片有专用的数字信号处理函数可调用，实现FIR滤波器相对简单，但是由于程序顺序执行，速度受到限制。而且，就是同一公司的不同系统的DSP芯片，其编程指令也会有所不同，开发周期较长。
（3）可编程FIR滤波器
　　还有一种是使用可编程逻辑器件，FPGA／CPLD。FPGA有着规整的内部逻辑块整列和丰富的连线资源，特别适合用于细粒度和高并行度结构的FIR滤波器的实现，相对于串行运算主导的通用DSP芯片来说，并行性和可扩展性都更好。
FIR的特点
　　有限长单位冲激响应（FIR）滤波器有以下特点：
（1）系统的单位冲激响应h (n)在有限个n值处不为零；
（2）系统函数H(z)在|z|>0处收敛，极点全部在z = 0处（因果系统）；
（3）结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。

3设计步骤
3.1 设计流程图
本课程设计主要是对一段语音信号，加入噪声后，用某种函数法设计出的FIR滤波器对加入噪声后的语音信号进行滤波去噪处理，并且分析对比前后时域和频域波形的程序设计。
程序的设计流程图如下图3.1所示：
开 始
用麦克风采集一段格式为.wav的语音

对语音信号进行采样，并加入一个单频噪声
对语音信号进行频谱分析，画出时域和频域波形图
用三角窗设计FIR滤波器
画出其频率响应
用FIR滤波器对语音信号进行滤波
比较滤波前后语音信号的波形及频谱
结 束

图3.1程序设计流程图

3.2录制语言信号
在左下角点击开始/程序/附件/娱乐/录音机，即可进入如图3.2.1所示的画面。

图 3.2 进入录音机界面
然后点击最右边的红色圆心的标志，即可进行录音。录一段2至3秒的语音，然后左键单击文件/属性，即可得到如下图3.2.2所示。

图 3.3 打开声音属性
在其中点击立即转换，即可得到如下图3.2.3所示界面，将声音格式设置为PCM，属性设置为8.000kHz，8位，单声道，7KB/秒，然后点击确定。然后再点击文件/另存为，即可将该语音文件存在MATLAB文件下的work文件夹里。

图 3.4 设置语音属性

3.3 语音信号的频谱分析
然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。采集完成后在信号中加入一个单频噪声。
程序如下：>>[x,fs,bits]= wavread('D:\renjian.wav');
>>sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放所录语音文件
>>N=length(x); % 计算信号x的长度
>>fn=1700; % 单频噪声频率，此参数可自己设计
>>t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率
>> x=x'; y=x'+0.3sin(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; % 计算频谱频率范围
图3.5 音乐信号加入噪声前后的时域频域图

3.4 设计数字滤波器和画出其频率响应
截止频率也可以任意自选，在单频噪声干扰附近即可。在这里，很重要的是通带截止db值的设置。这个值一定要根据我们使用的设计滤波器的方法来设定。因为我使用的是三角窗法，在三角窗函数中，滤波器的过渡带宽为6.1π/M，最小阻带衰减为25db。所以，一定要将通带截止db值设置的小于25db，所以，我将其设置为15db。
程序如下：
>> fp1=fn-200;fs1=fn-50;fp2=fn+200;fs2=fn+50;
>> Rp=1;As=15;% 带阻滤波器设计指标
>> df=min((fs1-fp1),(fp2-fs2)); % 计算上下边带中心频率，和频率间隔
>> dw=df/fs*2*pi;wc1=(fs1+fp1)/2;wc2=(fs2+fp2)/2; % 将Hz为单位的模拟频率换算为rad为单位的数字频率
>> ws1=fs1/fs*2*pi;ws2=fs2/fs*2*pi;wp1=fp1/fs*2*pi;wp2=fp2/fs*2*pi;
>>M=ceil(6.1*pi/dw)+1; % 计算三角窗设计该滤波器时需要的阶数
>> Rp=-min(db(wp1/dw+1:wp2/dw))
Rp =

0.1784
>> As=max(-round(db(ws1/dw+1:1:501)))
As =

22
>>n=0:M-1; % 定义时间范围
>> w_tri=(triang(M))'; % 产生M阶的三角窗
>> hd=ideal_lp(wc1,M)+ideal_lp(pi,M)-ideal_lp(wc2,M); % 调用自编函数计算理想带阻滤波器的脉冲响应
>> h=hd.*w_tri; % 用窗口法计算实际滤波器脉冲响应
>> [db,mag,pha,grd,w]=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性
通过该程序得到的图如下图3.4所示：

图3.6滤波器的特性波形图

3.5 用滤波器对信号进行滤波
在将滤波器设计好后，我们用自己设计的带阻滤波器对采集的语音信号进行滤波。在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。
程序如下：
>> y_fil=filter(h,1,y);% 用设计好的滤波器对加噪信号y进行滤波
>> Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半
语音信号的时域波形和频谱图比较如下图3.5所示：

图3.7 音乐信号的时域波形和频谱图比较

3.6 结果分析
从图3.5中我们可以看出原信号与滤波去噪信号的时域图基本相似，只有边缘部分有点差异；原信号与滤波去噪信号的频谱图波形也大致相似。通过观察可以看到，加噪信号的时域图中大部分都被加入的噪声给遮盖了，加噪信号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在2500Hz左右处有一个尖脉冲，而滤波去噪信号的频谱图中该尖脉冲已经消失，波形大致与原图相似，可见滤波去噪效果基本不错。
在将三个信号的时域波形和频谱图比较之后，我们还要通过回放去滤波去噪音乐信号，来跟原信号相比，以检验滤波器的效果。在Matlab中，函数sound可以对声音进行回放。其调用格式为：sound (x，fs，bits)。用sound(y_fil,fs,bits)语句回放该滤波去噪信号，便可以感觉到滤波后的语音信号与原信号差不多，该设计基本符合要求。

4出现的问题及解决方法
在课程设计的两周时间里，刚开始的前一周是在机房上机，老师跟我们综述本次课程设计的具体流程，在老师的引导下我们可以请教老师不懂的问题及出现错误等等。刚开始的时候是模仿老师给的模版一步步来，在遇到不懂的地方老师能及时的指导。在整个课程设计的过程中，出现的问题在老师与同学的帮助下，和自己的努力思考，问题得到了圆满解决。
首先，我用录音机录了一段语音文件，然后在Matlab中编写相应程序，但是在编译过程中程序报错，经过摸索，我发现我的语音文件不是单声道的，所以有进行了声道的调节，最终总算解决了这个问题。
语音文件添加好之后就是对语音添加噪声，进行频谱分析了，这些基本上都是以前上Matlab课程中接触到的一些基本知识，只是很多函数是老师给的，也给我们省下了不少时间。在编程过程中遇到和很多不懂得问题，但是在老师和同学们的帮助下，问题基本上都得到了解决。
确实，在本次课程设计过程中遇到了很多的困难，但是同学们的帮助让我本次设计省了不少心，老师在大体上给我们指明方向，同学们则在我实际动手的过程中给予指点和帮助，让我受益匪浅。

5 结束语
为期两周的课程设计终于结束了，在这次课程设计的过程中我又学到了很多，数字信号处理课程设计是数字信号处理课程的重要实践性环节，是我们在校期间一次较全面的能力训练，使得我巩固加深、综合运用所学的专业理论知识，对通信工程领域各种技术的DSP实现的设计有了一定的了解。通过语音信号滤波去噪——使用三角窗设计的FIR滤波器课程设计，我对数字信号方面的知识得到了更深一步的了解。在这次课程设计中，有许多的知识都不是很懂，通过课设中查阅资料等，我拓宽了知识面，增长了见识。最终在和同学们共同进步中圆满完成课程设计任务。
本次数字信号处理课程设计让我深深体会到做事情的时候都必须用心，要有耐心，使我通过自身的实践，对DSP的设计程序、内容和方法有更深入的掌握，提高实际运用的能力。并可综合运用这些知识解决一定的实际问题，使我在所学知识的综合运用能力上以及分析问题、解决问题能力上得到一定的提高。

参考文献
[1] 张圣勤.MATLAB7.0实用教程[M]. 北京：机械工程出版社，2008
[2] 郑国强，彭勃.数字信号处理：理论与实践[M].西安：西安电子科技大学出版社，2009
[3]吴镇扬.数字信号处理[M]. 北京：高等教育出版社，2004
[4] 高西全，丁玉美，阔永红. 数字信号处理——原理、实现及应用[M]. 北京：电子工业出版社，2009
[5] 程佩青. 数字信号处理教程[M]. 北京：清华大学出版社，2002

附录1：弹拨音乐滤波去噪——使用三角窗设计的FIR滤波器源程序清单

%程序名称：语音信号滤波去噪
%程序功能：用三角窗法设计FIR滤波器，并对加了噪声后的语音信号进行滤波去噪。
%程序作者：任健
%最后修改日期：2014-7-4

>> [x,fs,bits]= wavread('D：/renjian.wav');
>>sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放所录语音文件
>>N=length(x); % 计算信号x的长度
>>fn=2500; % 单频噪声频率，此参数可自己设计

>>t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率
>> x=x'; y=x'+0.03*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; % 计算频谱频率范围
>> subplot(2,2,1);plot(t,x);xlabel('时间（单位：s）');ylabel('幅度');title('原始语音信号');axis([0 1.1,-300 100]);
>> subplot(2,2,2);plot(f,X);xlabel('频率（单位：Hz）');ylabel('幅度谱');title('原始语音信号幅度谱') ; axis([0 1.1,0 2]);
>> subplot(2,2,3);plot(t,y);xlabel('时间（单位：s）');ylabel('幅度');title('加入单频干扰后的语音信号'); axis([0 1.1,-4 4]);
>> subplot(2,2,4);plot(f,Y);xlabel('频率（单位：Hz）');ylabel('幅度谱');title('加入干扰后的语音信号幅度谱');axis([0 4000,0 35]);
>> fp1=fn-200;fs1=fn-50;fp2=fn+200;fs2=fn+50;
>> Rp=1;As=15;% 带阻滤波器设计指标
>> df=min((fs1-fp1),(fp2-fs2)); % 计算上下边带中心频率，和频率间隔
>> dw=df/fs*2*pi;wc1=(fs1+fp1)/2;wc2=(fs2+fp2)/2; % 将Hz为单位的模拟频率换算为rad为单位的数字频率
>> ws1=fs1/fs*2*pi;ws2=fs2/fs*2*pi;wp1=fp1/fs*2*pi;wp2=fp2/fs*2*pi;
>>M=ceil(6.1*pi/dw)+1; % 计算三角窗设计该滤波器时需要的阶数
>> Rp=-min(db(wp1/dw+1:wp2/dw))
Rp =

0.1784
>> As=max(-round(db(ws1/dw+1:1:501)))
As =

22
>>n=0:M-1; % 定义时间范围
>> w_tri=(triangM))'; % 产生M阶的三角窗
>> hd=ideal_lp(wc1,M)+ideal_lp(pi,M)-ideal_lp(wc2,M); % 调用自编函数计算理想带阻滤波器的脉冲响应
>> h=hd.*w_tri; % 用窗口法计算实际滤波器脉冲响应
>> [db,mag,pha,grd,w]=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性
>> subplot(2,2,1);plot(w/pi,db);xlabel('w');ylabel('db'); axis([0 1.1,-300 100]);
>> subplot(2,2,2);plot(w/pi,mag);xlabel('w');ylabel('mag'); axis([0 1.1,0 2]);
>> subplot(2,2,3);plot(w/pi,pha);xlabel('w');ylabel('pha'); axis([0 1.1,-4 4]);
>> subplot(2,2,4);stem(n,w_tri);xlabel('n');ylabel('w(n)');
>> y_fil=filter(h,1,y);% 用设计好的滤波器对加噪信号y进行滤波
>> Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半
>> subplot(3,2,1);plot(t,x);xlabel('时间（单位：s）');ylabel('幅度');title('原始语音信号')
>> subplot(3,2,2);plot(f,X);xlabel('频率（单位：Hz）');ylabel('幅度谱');title('原始语音信号幅度谱') ;axis([0 6000,0 100])
>> subplot(3,2,3);plot(t,y);xlabel('时间（单位：s）');ylabel('幅度');title('加入单频干扰后的语音信号')
>> subplot(3,2,4);plot(f,Y);xlabel('频率（单位：Hz）');ylabel('幅度谱');title('加入干扰后的语音信号幅度谱');axis([0 6000,0 100])
>> subplot(3,2,5);plot(t,y_fil);xlabel('时间（单位：s）');ylabel('幅度');title('滤波后的语音信号')
>> subplot(3,2,6); plot(f,Y_fil);xlabel('频率（单位：Hz）');ylabel('幅度谱');title('滤波后的语音信号幅度谱');axis([0 6000,0 100])

附录2：ideal_lp函数源程序清单

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);

附录3： freqz_m函数源程序清单

function [db,mag,pha,grd,w] = 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))';
mag = 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;