音乐信号滤波去噪-基于Gaussian的频率采样型FIR滤波器设计

音乐信号滤波去噪

—基于Gausswin设计的FIR滤波器

学生姓名：王杰指导老师：高明

摘要本次课程设计主要内容是设计利用Gausssian窗口函数设计法设计一个FIR滤波器，将一段加入噪声干扰的音乐信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先在因特网上下载工具一段音乐信号，并加入预设好频率的单频噪声，然后对信号进行频谱分析得到加噪前后的信号对比图，然后设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱并进行分析。由分析结果得到，滤波器后的语音信号与原始信号基本一致，即设计的FIR滤波器能够去除信号中所加单频噪声，达到了设计目的。

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

1引言

本次课程设计主要是将一段音乐信号加入噪声，然后用某种函数法设计出的FIR滤波器对加入噪声后的音乐信号进行滤波去噪处理，处理时采用的是利用窗口设计法选择GAUSSWIN 设计的FIR 滤波器，通过课程设计了解FIR 滤波器设计的原理和步骤，掌握用Matlab语言设计滤波器的方法，观察音乐信号滤波前后的时域波形的比较，加深对滤波器作用的理解。

1.1 课程设计的目的

数字信号处理（Digital Signal Proccessing，简称DSP）是一门涉及许多学科而广泛应用于许多领域的新型学科。20世纪60年代后，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并并得到迅速的发展。在过去的二十多年时间里，数字信号处理已经在通信领域得到极为广泛的应用。数字信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、滤波、估值、增强、压缩、识别等处理，以得到符

合人们的需要的信号形式。

在本次课程设计中，最主要的设计是设计FIR 滤波器，FIR 滤波器的设计方法主要分为两类，第一类是基于逼近理想滤波器特性的方法，包括窗函数法、频域采样法和等波纹最佳逼近法；第二类是最优设计法。

本次的课程设计主要采用的是第一类设计方法，是利用GAUSSWIN法设计FIR滤波器对一段音乐信号进行滤波去噪，通过这一过程，对滤波前后波形进行对比分析得到结论。此课程设计比较简单，主要是将书本中的知识应用到现实中，并且根据自己对设计题目的理解，运用软件编写出程序实现这一设计，也是我们对数字信号处理的原理进行验证的一个过程。对此，也可以加深我们对所学知识的理解，培养我们的动手能力。

1.2 课程设计的要求

（1）滤波器指标必须符合工程实际。

（2）设计完后应检查其频率响应曲线是否满足指标。

（3）处理结果和分析结论应该一致，而且应符合理论。

（4）独立完成课程设计并按要求编写课程设计报告书。

1.3 设计平台MATLAB

MATLAB是由美国Math Works公司20世纪80年代中期推出的数学软件。MATLAB 是“Matric Laboratory”的缩写，意及“矩阵实验室”，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。Matlab已经发展成为多学科、多种工作平台的功能强大的大型软件。在欧美的高校和研究机构中，MATLAB是一种非常流行的计算机语言，许多重要的学术刊物上发表的论文均是用MATLAB来分析计算以及绘制出各种图形。

MATLAB是一完整的并可扩展的计算机环境，是一种进行科学和工程计算的交互式程序语言。它的基本数据单元是不需要指定维数的矩阵，它可直接用于表达数学的算式和技术概念，而普通的高级语言只能对一个个具体的数据单元进行操作。因此，解决同样的数值计算问题，使用MATLAB 要比使用Basic、Fortran 和C 语言等提高效率许多倍。许多人赞誉它为万能的数学“演算纸”。MATLAB 采用开放式的环境，你可以读到它的算法，并能改变当前的函数或增添你自己编写的函数

MATLAB 包含的内容非常丰富，功能强大，可以概括为以下几个方面：

(1)可以在多种操作系统下运行，如DOS、Windows 95/98/2000/2000/NT、CompaqAlpha、LinuxSun Solaris 等。

(2)有超过500 种的数学、统计、科学及工程方面的函数，使用简单快捷，并且有很强的用户自定义函数的能力。

(3)有强大的图形绘制和可视化功能，可以进行视觉数据处理和分析，进行图形、图像的显示及编辑，能够绘制二维、三维图形，使用户可以制作高质量的图形，从而写出图文并茂的文章。

(4)有从外部文件及外部硬件设备读入数据的能力。

(5)有丰富的工具箱〔toolbox〕。各个领域的专家学者将众多学科领域中常用的算法编写为一个个子程序，即m 文件，这些m 文件包含在一个个工具箱中。其工具箱可以分为两大类，即功能性工具箱和科学性工具箱。功能性工具箱主要用来扩充MATLAB 的符号计算、图形可视化、建模仿真、文字处理等功能以及与硬件实时交互的功能。学科性工具箱是按学科领域来分类的，如信号处理、控制、通信、神经网络图像处理、系统辨识、鲁棒控制、模糊逻辑、小波等工具箱。

MATLAB 中的信号处理工具箱内容丰富，使用简便。在数字信号处理中常用的算法，如FFT，卷积，相关，滤波器设计，参数模型等，几乎都只用一条语句即可以调用。数字信号处理所常用的函数有波形的产生、滤波器的分析和设计、傅里叶变换、Z 变换等。

2 设计原理

2.1 FIR滤波器

滤波器根据其冲激响应函数的时域特性，可分为2 种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。FIR 和IIR 的滤波原理都是进行卷积，就是对输入信号进行某种计算。FIR 用处就在于对数字信号进行必要的处理，得到所需的输出信号。

FIR 系统有自己突出的优点：系统总是稳定的；易实现线性相位；允许设计多通带（或多阻带）滤波器，后两项都IIR系统不易实现的。FIR 数字滤波器的设计方法有多种，如窗函数设计法、频率采样法和Chebyshev逼近法等。随着Matlab软件尤其是Matlab

以使设计达到最优化。

FIR 数字滤波器设计的基本步骤如下：

（1）确定技术指标

在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以2种方式给出。第一种是绝对指标。他提供对幅度响应函数的要求，一般应用于FIR 滤波器的设计。第二种指标是相对指标。他以分贝值的形式给出要求。本文中滤波器的设计就以线性相位FIR 滤波器的设计为例。

（2）逼近

确定了技术指标后，就可以建立一个目标的数字滤波器模型（通常采用理想的数字滤波器模型）。之后，利用数字滤波器的设计方法（窗函数法、频率采样法等），设计出一个实际滤波器模型来逼近给定的目标。

（3）性能分析和计算机仿真

上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。

2.2 窗口设计法

窗口法设计的基本想法是要选取某一种合适的理想频率选择性滤波器（这种滤波器总是有一个非因果，无限长的脉冲响应），然后将它的脉冲响应截断（或加窗）以得到一个线性相位和因果的FIR 滤波器。因此，这种方法的重点在于选择某种恰当的窗函数和一个合适的理想滤波器。现用()j d H e ω代表一理想频率选择性滤波器，它在整个通带内有单位幅度增益和线性相位特性，而阻带内有零响应。

用窗口设计法基本步骤如下[1]：

(1) 给定要求的理想频率响应()j d H e ω,一般给定分段常数的理想频率特性。 (2) 由于是在时域设计故必须求出()d h n

1

()[()]()2j j j d d

h n IDTFT H e H

e e d π

ω

ωωπω

π

-

==

? (2-1)

以截断(相乘)，窗的点数是N 点。 截断后的序列为h(n)：

()()(),01d h n h n n n N ω=≤≤- (2-2)

窗的点数N 及窗的形状是两个极重要的参数。 (4) 求出加窗后实际的频率响应()j H e ω

()

()[()][()()]11

[()*()]()()22j d j j j j d d

H e DTFT h n DTFT h n n H e W e H

e W e

d ωπ

ωωθωθπωθ

ππ

--

====?(2-3)

(5) 检验()j H e ω是否满足()j d H e ω的要求，不满足，则需考虑改变窗形状或改变窗长的点数N ，重复第(3)、(4)两步，到满足要求为止。 常见的窗函数性能表如下表2-1所示。

表2-1常见的窗函数性能表

2.3 GAUSSIAN 窗

高斯窗是一种指数窗。主瓣较宽，故而频率分辨力低；无负的旁瓣，第一旁瓣衰减达一55dB 。常被用来截短一些非周期信号，如指数衰减信号等。对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。 Gausswin 的时域表达式可以表示为：

212()2/2

[1]k N N w k e

α--+=(2-4)

其中k 与α的取值范围必须满足0,2k N α≤≤≥[2]

2.4滤波器结构

本次课程设计的滤波器采用的是频率采样型结构。 (1) 频率采样型结构的导入

若FIR DF 的冲激响应为有限长(N 点)序列h(n)，则有：

图2-1 关于h(n)的推导

所以，对h(n)可以利用DFT 得到H(k),再利用内插公式：

11

01()

()(1)1N N

k k N H k H z z N W z

----==--∑(2-5) 来表示系统函数。 (2) 频率采样型滤波器结构

由式2-5得到FIR 滤波器的另外一种结构：频率采样型结构。它是由两部分级联而成。

1

1()()()N c k k H

z H z H z N -==∑ (2-6)

其中：级联中的第一部分为梳状滤波器

()(1)N

c H z z -=-， 第二部分由N 个谐振器组成的谐振贵。

1

()

()1k k N H k H z W z

--=

-(2-7) (3) 频率采样型结构流图如图2-2所示

图2-2 频率采样型结构流程图

(4) 频率采样型结构特点

它的系数H(k)直接就是在2k w k N

π

=

处的频率响应。因此，控制滤波器的频率响应是很直接的。但所有的相乘系数及H(k)都是复数，应将它们先化成二阶的实数，这样乘起来比较复杂，增加了乘法次数及储存量；且所有的谐振器的极点都是在单

位圆上，由

k

N W -决定。考虑到系数量化的影响，当系数量化时，极点会移动，有些

极点就不能被梳状滤波器的零点锁抵消。(零点由延时单元绝对，不受量化的影响)系统就会变得不稳定[1]。

3设计步骤

3.1设计流程图

滤波器设计流程图如图3-1所示。

图3-1 滤波器设计流程图

3.2下载并截取音乐信号

从网上下载一段音乐，从中截取一段格式为.wav 的语音信号，时间为2—3秒，并加入噪声。在MALAB 平台上，观察原始语音信号与加入噪声后的时域和频谱图。

源程序如下所示：

[x,fs,bits]=wavread('G:\数字信号处理课程设计_rE\funk.wav');%读取音乐信号的数据sound(x,fs,bits);

t=0:1/fs:(N-1)/fs; %计算时间范围，样本除以采样频率

x=x';

fn=2436; %单频噪声功率

y=x+0.1*sin(fn*2*pi*t); %加单频噪声

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-2所示，由时频域分析图可以看出，音乐信号能量主要集中在1000Hz以内。

图3-2 原始信号与含噪信号的时频域分析图

3.3滤波器设计

截取好原始信号，接下来的工作是设计一个Gaussian窗滤波器。

利用公式[1]：

RP=-min(db([1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end])) (3-1)

AS=-max(db(wsd/dww+1:wsu/dww+1)) (3-2) 分别得到通带波纹及通带最小阻带衰减，并输入相应滤波器的频率参数，接着用freqz_m求得滤波器的响应频率特性。

源程序如下所示:

fp1=1836;fs1=2406;fs2=2476;fp2=3036; %带阻滤波器设置指标

Rp=6;As=28;

fcd=(fp1+fs1)/2; %计算上下边带中心频率

fcu=(fp2+fp2)/2;

df=min((fs1-fp1),(fp2-fs2));

wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率

wsd=fs1/fs*2*pi;wsu=fs2/fs*2*pi;

M=ceil(5.8*pi/dw)+1; %计算Gauss窗设计该滤波器时需要的阶数

n=0:M-1; % 定义时间范围gausswf=gausswin(M); %利用gausswin函数产生一个Gaussian窗函数

hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); %计算理想阻带滤波器的脉冲响应

h_bs=gausswf'.*hd_bs; %窗口法计算实际滤波脉冲响应

[db,mag,pha,grd,w]=freqz_m(h_bs,[1]); %计算滤波器的频率特性

wp1=fp1/fs*2*pi;wp2=fp2/fs*2*pi;dww=2*pi/1000;

RP=-min(db([1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end])) %计算通带最大衰减

AS=-max(db(wsd/dww+1:wsu/dww+1)) %计算出阻带最小衰减滤波器幅度响应图、相位响应图、脉冲响应图如图3-3所示。由图可知，滤波器的实际通带最大衰减高于设定的通带最大衰减，而实际的阻带最小衰减低于设定的阻带最小衰减，性能达标。

图3-3滤波器性能分析图

3.4信号滤波处理

设计好滤波器后，我们要对语音信号进行滤波，对比滤波前和滤波后的音乐信号。

源程序如下所示：

y_fil=fftfilt(h_bs,y); %用设计好的滤波器对y进行滤波

Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:length(Y_fil)/2);% %计算频谱取前一半

sound (y_fil,fs,bits); %音乐信号回放wavwrite(y_fil,fs,'d:\funk_filted.wav'); %产生去噪后的音乐文件原始音乐信号、滤波后信号的时域、频域分析图，纵坐标为对数坐标的去噪前后的幅度频谱图对比图如图3-4所。

图3-4 音乐信号时域频域分析图

3.5结果分析：

通过观察上图，滤波后的语音信号发生了衰减，说明滤波器起到了滤波作用，同时通过频谱对比，可以看出滤波器滤掉了一部分频率范围内的信号。分别听原始语音和滤波后的语音信号，发现滤波后的语音信号噪声减小了，同时原始信号强度稍有减弱，基本达到了滤波的效果。

3.6滤波器结构设计

本次课程设计采用的是频率采样型的滤波器结构，将按照之前步骤设计出的直接型滤波器通过tf2fs函数转换成频率采样型[1]。

源程序如下所示：

%直接型到频率抽样型的转换

%----------------------

%[C,B,A]=tf2fs(h)

%C=包含各并联部分增益的行向量%B=包含按行排列的分子系数矩阵%A=包含按行排列的分母系数矩阵%h=FIR滤波器单位冲激响应向量%----------------------

N=length(h);

H=fft(h,N);

r=input('r=');

magH=abs(H);phaH=angle(H)';

%check even or odd N

if (N==2*floor(N/2))

L=N/2-1;

A1=[1,-r,0;1,r,0];

C1=[real(H(1)),real(H(L+2))]; else

L=(N-1)/2;

A1=[1,-r,0];

C1=[real(H(1))];

end

k=[1:L]';

%初始数组B和数组A

B=zeros(L,2);A=zeros(L,3);

%计算分母系数

A(1:L,1)=1;

A(1:L,2)=-2*r*cos(2*pi*k/N);

A(1:L,3)=r.^2

%计算分子系数

B(1:L,1)=cos(phaH(2:L+1));

B(1:L,2)=-cos(phaH(2:L+1)-(2*pi*k/N)); %计算增益系数 C=[2*magH(2:L+1),C1]';

根据计算得出的C(各并联部分增益的行向量),B(分子系数向量),A(分母系数向量)画出频率采样型滤波器结构图，如图3-4所示，C,B,A 参数如图3-5所示。

图3-4 频率采样型滤波器结构图

()

y n

图3-5滤波器的C，B，A参数表

4出现的问题及解决方法

1、下载得到的音乐文件为.mp3文件，用录音机软件转换格式时无法识别，后来用Adobe

Audition 3.0软件将下载得到的音乐文件转换成.wav文件就能顺利转换文件格式了。

2、在画原始信号，加噪后的信号及滤波信号的时频域分析图时，发现画出的图形坐标

系上下限不一致，导致不能产生直观的对比效果。经过参考matlab的help文档之后得知可以采用axis[]函数手动确定坐标系上下限，直观地画出了对比图。

3、在生成加噪之后的音乐文件时，按照老师的提示使用wavwrite()函数生成，但是在括

号里填写好目标存储目录之后发现matlab报错，不能成功的生成。后来经常参考matlab的help文件，得知wavwrite()函数需要填写源数据作为参数，于是将加噪之后的音频源数据y作为参数，填写到wavwrite()成功生成了加了噪声之后的文件。4、在调用自编函数计算理想阻带滤波器的脉冲响应及计算滤波器的频率特性时，发现

ideal_lp函数以及freqz_m函数报错，仔细检查才发现没有将已经编译好的m文件放入matlab的m文件夹内，导致无法调用自编函数。将编译好的函数m文件复制进目标文件夹后，顺利地进行了计算。

5、设计好滤波器后，对滤波器的性能进行验证时，不知道该如何检验。后来按照老师

所给的提示，通过公式计算出滤波器的实际Rp以及As与预设的Rp和As进行比较并且在性能曲线图中画出作为参考的预设的Rp,As与实际的性能曲线相比较判断性能是否达标，经对比发现滤波器符合设计标准，顺利的完成了滤波器性能的检验6、进行的最后的滤波器结构图设计时，发现由tf2fs函数计算得出的C,B,A矩阵长度过

长，达到了几百个，将会使得结构图的绘制非常困难。经过仔细调整通带和阻带截止频率将高斯窗的阶数设定在了43阶，将C,B,A矩阵减少到了22行，简化了滤波器的结构图，使得绘制顺利的完成了。

5结束语

这次课程设计，是运用上学期的所学数字信号处理理论知识的一次实践，加深了对所学知识的理解，所以这次课程设计对我来说是收获颇多的。

本次的课程设计，我选到的课题主要内容是通过数字信号处理技术实际处理信号，这里就是指对音乐信号的处理。作为存储于计算机中的音乐信号，其本身就是离散化了的向量，我们只需将这些离散的量提取出来，就可以对其进行处理了。

本次课程设计的主要设计平台是MATLAB，通过MATLAB 里调用几个命令函数，很轻易的在实际化音乐与数字信号的理论之间搭了一座桥。课题的特色在于它将音乐看作了一个向量，于是音乐数字化了，则可以完全利用数字信号处理的知识来解决。我们可以像给一般信号做频谱分析一样，来给音乐信号做频谱分析，也可比较容易地用数字滤波器来对音乐进行滤波处理。改变参数，理论结合实际，分析各参数对图形的影响，从而加深对各个参数的理解。在完成这次课程设计过程中学到了许多东西，进一步理解了滤波器设计方法和各参数意义，通过分析信号时域和频域的关系等，加深了对滤波性能的理解，而且学会了使用Matlab一些基本函数，增加了进一步学习Matlab软件的兴趣。同时，通过本次课程设计，锻炼了我的动手能力，和提高了我分析问题，解决问题的能力。

同学们在一起思考问题，通过亲自动手设计一个FIR 滤波器，发现了许多潜在的问题，而这些问题是在平常被我们所忽略的，甚至认为是不成问题的。让我再次感受到我们应当把所学知识和实践相结合，才能够提高自己的知识水平，取得更大的发展。

参考文献

[1]程佩青.数字信号处理教程[M].北京：清华大学出版社，2013:307,476,278-283.

[2]yangpan011.几种常见的窗函数以及其Matlab实现[EB/OL].https://www.wendangku.net/doc/461558210.html,/yangpan011/article/details/52247355?locationNum=5&fps= 1,2017-2-24.

[3]林爱英, 谷小青, 郑宝周. 用频率采样法设计FIR滤波器[J]. 现代电子技术, 2010, 33(17):85-87.

附录1：音乐信号滤波去噪设计源程序清单

% 程序名称：funk_filter.m

% 程序功能：采用基于gausswin的窗口设计法，设计FIR 滤波器对含噪音乐信号进行滤波去噪处理。

% 程序作者：王杰

% 最后修改日期：2017-2-28

[x,fs,bits]=wavread('G:\数字信号处理课程设计_rE\funk.wav');%读取音乐信号的数据sound(x,fs,bits);

N=length(x); %计算信号x的长度

t=0:1/fs:(N-1)/fs; %计算时间范围，样本除以采样频率

x=x';

fn=2436; %单频噪声功率

y=x+0.1*sin(fn*2*pi*t); %加单频噪声

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; %计算频谱的频率范围

figure(1);

subplot(2,2,1),plot(t,x);grid on;

axis([0,2.5,-2,2]);title('原始音乐信号');xlabel('时间（单位：s）');ylabel('幅度');

subplot(2,2,2),plot(f,X);grid on;

title('音乐信号幅度谱图');xlabel('频率（单位：Hz）');ylabel('幅度谱');

subplot(2,2,3),plot(t,y);grid on;

axis([0,2.5,-2,2]);title('加入噪声干扰后的音乐信号');xlabel('时间（单位：s）');ylabel('幅度');

subplot(2,2,4),plot(f,Y);grid on;

title('加入噪声干扰后的音乐信号幅度谱图');xlabel('频率（单位：Hz）');ylabel('幅度谱'); wavwrite(y,fs,'d:\funk_converted.wav');

Rp=6;As=28;

fcd=(fp1+fs1)/2; %计算上下边带中心频率

fcu=(fp2+fp2)/2;

df=min((fs1-fp1),(fp2-fs2));

wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; %将Hz为单位的模拟频率换算为rad为单位的数字频率

wsd=fs1/fs*2*pi;wsu=fs2/fs*2*pi;

M=ceil(5.8*pi/dw)+1; %计算Gauss窗设计该滤波器时需要的阶数

n=0:M-1; % 定义时间范围

gausswf=gausswin(M); %利用gausswin函数产生一个Gaussian窗函数

hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); %计算理想阻带滤波器的脉冲响应

h_bs=gausswf'.*hd_bs; %窗口法计算实际滤波脉冲响应

[db,mag,pha,grd,w]=freqz_m(h_bs,[1]); %计算滤波器的频率特性

wp1=fp1/fs*2*pi;wp2=fp2/fs*2*pi;dww=2*pi/1000;

RP=-min(db([1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end])) %计算通带最大衰减

AS=-max(db(wsd/dww+1:wsu/dww+1)) %计算出阻带最小衰减

figure(2);

subplot(2,2,1);plot(w/pi,db);xlabel('w/pi');ylabel('db');title('FIR滤波器的幅度响应图');axis([0.3 0.9 -80 5]);

line([0,1.5],[-RP,-RP],'Color','r','LineWidth',2,'LineStyle','--');

line([0,1.5],[-AS,-AS],'Color','r','LineWidth',2,'LineStyle','--');

line([wp1/pi,wp1/pi],[-100,10],'Color','r','LineWidth',2,'LineStyle','--');

line([wp2/pi,wp2/pi],[-100,10],'Color','r','LineWidth',2,'LineStyle','--');

line([wsd/pi,wsd/pi],[-100,10],'Color','r','LineWidth',2,'LineStyle','--');

line([wsu/pi,wsu/pi],[-100,10],'Color','r','LineWidth',2,'LineStyle','--'); grid on;

subplot(2,2,2);plot(w/pi,mag);xlabel('w/pi');ylabel('幅度mag');title('FIR滤波器的幅度响应图');axis([0 1 -0.2 1.2]);grid on;

subplot(2,2,3);plot(w/pi,pha);xlabel('w/pi');ylabel('相位pha');title('滤波器相位响应图