数字信号处理实验报告

前言

《数字信号处理》是信息电子，通信工程等本科专业及其他相近专业的一门专业必修课。通过本课程的学习，学生应掌握以下基本概念、理论和方法：采样定理、离散序列的变换、离散信号的频谱分析；离散系统的传递函数、频率响应、离散系统的基本分析方法；数字滤波器的设计理论、滤波器的软件实现；离散傅立叶变换理论、快速傅立叶变换方法；有限字长效应。

为了使学生更好地理解和深刻地把握这些知识，并在此基础上，训练和培养学生掌握离散系统的基本概念和分析方法，数字滤波器的设计和实现，以及如何利用快速傅立叶变换等DSP技术对数字信号进行分析、滤波等处理，设置了以下三个实验：

（1）离散时间序列卷积和MATLAB实现；

内容：使用任意的编程语言编制一个程序，实现两个任意有限序列的卷积和。

目的：理解线性非移变系统I/O关系和实现

要求：掌握使用计算机实现数字系统的方法

（2）FFT算法的MATLAB实现；

内容：使用MATLAB编程语言编制一个程序，实现任意有限序列的FFT。

目的：理解FFT算法的意义和实现

要求：掌握使用计算机实现FFT算法的方法

（3）数字滤波器的设计；

内容：使用MATLAB编程语言编制一个程序，实现FIR或IIR滤波器的设计目的：理解数字滤波器的设计技术

要求：掌握使用计算机进行数字滤波器设计的方法

（4）窗函数设计FIR滤波器；

内容：使用MATLAB编程语言编制一个程序，实现FIR或IIR滤波器的设计目的：理解数字滤波器的设计技术

要求：掌握使用计算机进行数字滤波器设计的方法

该实验指导书是参照该课程的教学大纲而编制的，适合于信息电子工程、通信工程等本科专业及其他相近专业。

数字信号知识预备：

一．典型的离散信号 1．单位抽样信号

?

?

?≠==000

1)(n n n δ

2．脉冲串序列p(n)

将)(n δ在时间轴上延迟k 个抽样周期，得

?

?

?≠==-k n k

n k n 01)(δ

若k 从∞-变到∞+，则的所有移位可形成一个无限长的脉冲串序列p(n)

∑∞

-∞

=-=

k k n n p )()(δ

3．单位阶跃序列

??

?<≥=000

1)(n n n u 4．正弦序列 5．复正弦序列 6．指数序列

二．离散信号的运算 1．信号的延迟

给定离散信号x(n),若信号)(1n y ,)(2n y 分别定义为

)(1n y ＝x(n-k) )(2n y =x(n+k)

序列下x(n)在某一时刻k 时的值可用)(n δ的延迟来表示，即

δ

(k

)

x)

n-

(k

x=)

(n

2.两个信号的相加和相乘

x(n)=x1(n)+x2(n)

y(n)=x1(n)x2(n)

y(n)=cx(n)

3.信号时间尺度的变化

y(t)=x(t/a),式中a>0

若a>1或a<1呢？

4.信号的分解

5．信号的变换

三．信号的分类

1．连续时间信号和离散时间信号

2．周期信号和非周期信号

3．确定性信号和随机信号

4．能量信号和功率信号

5．一维信号、二维信号及多通道信号

四．噪声

x(n)=s(n)+u(n)

x(n)=s(n) u(n)

五．离散时间系统的基本概念

y(n)=T[x(n)]

例：一个离散时间系统的输入、输出关系是y(n)=ay(n-1)+x(n)

求所给系统的单位抽样响应。例：系统

∑=-

=

2

)

(

)

(

)

(

k

k

n

x

k

b

n

y

式中)0(b,)1(b,)2(b为常数。求所给系统的单位抽样响应

六．有关离散系统的几个重要定义

1．线性

对ax1(n)+bx2(n)的响应是ay1(n)+by2(n)

2．移不变性

T[x(n)]=y(n)

T[x(n-k)]=y(n-k)

3．因果性

4．稳定性

七．LSI系统的输入输出关系

例：令}1,1{

)}

1(

),

0(

{

)

(=

=h

h

n

h，}4,3,2,1{

)}

3(

)0(

{

)

(=

?

=x

x

n

x，，求)

(n

x和)

(n

h的线性卷积。

实验一：离散时间序列卷积和MATLAB 实现

实验学时：2 实验类型：验证

实验要求：必修

（一）实验目的：学会用MATLAB 对信号与系统分析的方法，理解离散序列卷积和的计算对进行离散信号与系统分析的重要性。

（二）实验原理：

1、离散时间序列f1(k)和f2(k)的卷积和定义：

f(k)=f1(k)*f2(k)=

∑∞

-∞

=-?

i i k f i f )(2)(1

2、在离散信号与系统分析中有两个与卷积和相关的重要结论：

a 、f(k)=

∑

∞

-∞

=-?i i k i f )()(δ=f(k)* δ(k)即离散序列可分解为一系列

幅度由f(k)决定的单位序列δ(k)及其平移序列之积。

b 、对线性时不变系统，设其输入序列为f(k)，单位响应为h(k)，其零状态

响应为y(k)，则有：y(k)=

∑

∞

-∞

=-?i i k h i f )()(

（三）实验内容：

conv.m 用来实现两个离散序列的线性卷积。其调用格式是：

y=conv(x,h)

若x 的长度为N ，h 的长度为M ，则y 的长度L=N+M-1。

题一：令x(n)= {

}5,4,3,2,1，h(n)＝{}246326，，，，，，y(n)=x(n)*h(n)，求y(n)。 要求用subplot 和stem 画出x(n),h(n),y(n)与n 的离散序列图形。

题二：已知序列

f1(k)=???≤≤其它0201k f2(k)=????

??

?===其它

33

2211k k k

调用conv()函数求上述两序列的卷积和

题三：编写计算两离散序列卷积和f(k)=f1(k)*f2(k)的实用函数dconv().要求该程序在计算出卷积和f(k)的同时，还绘出序列f1(k),f2(k)和f(k)的时域波形图，并返回f(k)的非零样值点的对应向量。 function[f,k]=dconv(f1,f2,k1,k2)

%f1(k),f2(k)及f(k)的对应序号向量分别为k1，k2和k 。

题四：试用MATLAB 计算如下所示序列f1(k)与f2(k)的卷积和f(k)，绘出它们的时域波形，并说明序列f1(k)与f2(k)的时域宽度与序列f(k)的时域宽度的关系。提示：可用上述dconv()的函数来解决。

f1(k)=??????

?==-=其它

11021

1

k k k

f2(k)=??

?≤≤-其它

221k

题五：已知某LTI 离散系统，其单位响应h(k)=e(k)-e(k-4)，求该系统在激励为f(k)=e(k)-e(k-3)时的零状态响应，并绘出其时域波形图。 提示：可用dconv()的函数来解决。

（四）实验报告：

1．根据实验原理，编写代码，得出实验结果，并画出波形图。 2．归纳、总结实验结果。 3．心得体会及其他。

实验二：FFT 算法的MATLAB 实现

实验学时：3 实验类型：验证

实验要求：必修

（一）实验目的：理解离散傅立叶变换时信号分析与处理的一种重要变换，特别是FFT 在数字信号处理中的高效率应用。

（二）实验原理：

1、有限长序列x(n)的DFT 的概念和公式：

???

????-≤≤=-≤≤=∑∑-=--=1

01

01

0)(1)(1

0)()(N k kn N N n kn N N n W k x N n x N k W n x k x

)/2(N j N e W π-=

2、FFT 算法 调用格式是

X= fft(x)

或 X=fft(x,N)

对前者，若x 的长度是2的整数次幂，则按该长度实现x 的快速变换，否则，实现的是慢速的非2的整数次幂的变换；对后者，N 应为2的整数次幂，若x 的长度小于N ，则补零，若超过N ，则舍弃N 以后的数据。Ifft 的调用格式与之相同。 例1： N=8;

n=0:N-1;

xn=[4 3 2 6 7 8 9 0];

Xk=fft(xn)

→

Xk =

39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.000

0 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i

Xk与xn的维数相同，共有8个元素。

例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。

clf;

fs=100;N=128; %采样频率和数据点数

n=0:N-1;t=n/fs; %时间序列

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号

y=fft(x,N); %对信号进行快速Fourier变换

mag=abs(y); %求得Fourier变换后的振幅

f=n*fs/N; %频率序列

subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=128');grid on;

subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');

ylabel('振幅');title('N=128');grid on;

%对信号采样数据为1024点的处理

fs=100;N=1024;n=0:N-1;t=n/fs;

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号

y=fft(x,N); %对信号进行快速Fourier变换

mag=abs(y); %求取Fourier变换的振幅

f=n*fs/N;

subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=1024');grid on;

subplot(2,2,4)

plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');

ylabel('振幅');title('N=1024');grid on;

运行结果：

https://www.docsj.com/doc/c210040246.html,/photo/kHUiRaoKm75k5hcPSt06ug==/3383892170016551016.j pghttps://www.docsj.com/doc/c210040246.html,/photo/1K5gl1MolBNNUF84nRAnjg==/3383892170016551 031.jpg

fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT 变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。

（三）实验内容：

题一：若x(n)=cos(n*pi/6)是一个N=12的有限序列，利用MATLAB计算它的DFT并画出图形。

题二：一被噪声污染的信号，很难看出它所包含的频率分量，如一个由50Hz 和120Hz正弦信号构成的信号，受均值随机噪声的干扰，数据采样率为1000Hz，通过FFT来分析其信号频率成分，用MATLAB实现。

题三：调用原始语音信号mtlb，对其进行FFT变换后去掉幅值小于1的FFT 变换值，最后重构语音信号。

（要求有四幅语音信号的频谱图在同一图形窗口以便比较：分别是1、原始语音信号；2、FFT变换；3去掉幅值小于1的FFT变换值；4、重构语音信号）

（四）实验报告：

1．根据实验原理，编写代码，得出实验结果，并画出波形图。

2．归纳、总结实验结果。

3．心得体会及其他。

实验三：数字滤波器的设计

实验学时：3 实验类型：设计

实验要求：必修

（一）实验目的：掌握IIR 数字低通滤波器的设计方法。

（二）实验原理： 1、滤波器的分类

滤波器分两大类：经典滤波器和现代滤波器。

经典滤波器是假定输入信号)(n x 中的有用成分和希望取出的成分各自占有不同的频带。这样，当)(n x 通过一个线性系统（即滤波器）后可讲欲去除的成分有效的去除。

现代滤波器理论研究的主要内容是从含有噪声的数据记录（又称时间序列）中估计出信号的某些特征或信号本身。

经典滤波器分为低通、高通、带通、带阻滤波器。每一种又有模拟滤波器（AF ）和数字滤波器（DF ）。对数字滤波器，又有IIR 滤波器和FIR 滤波器。

IIR DF 的转移函数是：

∑∑=-=-+==

N k k

k M

r r

r z a z

b z X z Y z H 1

01)

()

()(

FIR DF 的转移函数是：

∑-=-=1

0)()(N n n z n h z H

FIR 滤波器可以对给定的频率特性直接进行设计，而IIR 滤波器目前最通用的方法是利用已经很成熟的模拟滤波器的设计方法进行设计。

2、滤波器的技术要求

低通滤波器：p ω：通带截止频率（又称通带上限频率） s ω：阻带下限截止频率

p α：通带允许的最大衰减 s α：阻带允许的最小衰减 （p α，s α的

单位dB ）

p Ω：通带上限角频率 s Ω：阻带下限角频率

（s p p T ω=Ω，s s s T ω=Ω）即 C p

p F ωπ

2=Ω C s

s

F ωπ2=Ω

3、IIR 数字滤波器的设计步骤：

1）按一定规则将给出的数字滤波器的技术指标转换魏模拟低通滤波器的技术指标。

2）根据转换后的技术指标设计模拟低通滤波器)(s G ； 3）再按一定的规则将)(s G 转换成)(z H 。

4）若是高通、带通或带阻数字滤波器则将它们的技术指标先转化为低通模拟滤波器的技术指标，然后按上述步骤2）设计出低通)(s G ，再将)(s G 转换为所需的)(z H 。

4．几种不同类型的滤波器的介绍：

因为我们设计的滤波器的冲击响应一般都为实数，所以有

2

*)()()()()(s j s j G s G s G s G s G Ω=-=Ω=

这样，如果我们能由p α，p Ω，s α，s Ω求出2

)(Ωj G ，那么就容易得到所需要的)(s G 。

不同类型的2

)(Ωj G 的表达式，代表了几种不同类型的滤波器。 （1）巴特沃思(Butterworth)滤波器：

n

C j G )

(11

)(222

Ω+=

Ω C 为待定常数，N 为待定的滤波器阶次。 （2）切比雪夫I 型(Chebyshev – I )滤波器：

)

(11

)(2

22

Ω+=

Ωn C j G ε 5．巴特沃思模拟低通滤波器的设计

由于每一个滤波器的频率范围将直接取决于设计者的应用目的，因此必然是千差万别。为了使设计规范化，我们需要将滤波器的频率参数作归一化处理。设所给的实际频率为Ω（或f ），归一化后的频率为λ，对低通模拟滤波器，令

λ＝p ΩΩ/

显然，1=p λ，p s s ΩΩ=/λ。又令归一化复数变量为p ，λj p =，显然

p p s j j p Ω=ΩΩ==//λ

所以巴特沃思模拟低通滤波器的设计可按以下三个步骤来进行。 (1)将实际频率Ω规一化 (2)求C 和N

11010/2-=P C α

s p s N λααlg 1

10

110lg

10

/10/--=

这样C 和N 可求。

若令p α＝3dB ，则C ＝1，这样巴特沃思滤波器的设计就只剩一个参数N ，这时

N

p N

j G 222

)/(11

11)(ΩΩ+=

+=

λλ

(3)确定)(s G

因为λj p =，根据上面公式有

N

N N p

j p p G p G 22)1(11

)/(11)()(-+=+=

- 由

0)1(12=-+N N p 解得

)221

2exp(πN

N k j

p k -+=，k ＝1，2，

(2)

这样可得

1

)21

2cos(21

)

)((1

)(21+-+-=

--=

-+πN

N k p p p p p p p G k N k k

求得)(p G 后，用p s Ω/代替变量p ，即得实际需要得)(s G 。 6．用双线性Z 变换法设计IIR 数字低通滤波器 s 平面到z 平面的映射关系

1

1

2+-=

z z T s s

称为双线性Z 变换，由此关系求出

s

T s

T z s s )2/(1)2/(1-+=

及 )

2/cos()

2/sin(2ωωs T j

j =Ω

即 )2/tan(2

ωs

T =

Ω

)2/arctan(2s T Ω=ω

因为设计滤波器时系数s

T 2

会被约掉，所以又有 s

s

z -+=

11 )2/tan(ω=Ω

Ω=arctan 2ω

（三）实验内容：

题一：试用双线性Z 变换法设计一低通数字滤波器，给定技术指标是

100=p f Hz ，300=s f Hz ，3=p αdB ，20=s α dB ，抽样频率1000=s F Hz 。

提示：首先应该得到角频率ω，然后再

（1） 将数字滤波器的技术要求转换为模拟滤波器的技术要求。 （2） 设计低通滤波器)(s G

由λ＝p ΩΩ/依次求出p λ，s λ，再求出N ，可得)(p G 然后由p

s p p G s G Ω=

=)()(转换成)(s G

（3） 由)(s G 求)(z H 设计步骤： （1）求出角频率

s

p

p F f w π

2=

s s

s F f w π

2=

（2）对角频率做预畸变)2/tan(p p ω=Ω )2/tan(s s ω=Ω

（3）求出模拟低通滤波器的阶次，利用函数[N ，Wn]=buttord(Wp ，Ws ，Rp ，Rs ，’s’) 注意：Wp ，Ws 应该为（2）中的p Ω，s Ω。

（4）设计模拟低通原型滤波器)(p G ，其调用格式是 [z ，p ，k]=buttap(N)。 N 是欲设计的低通原型滤波器的阶次，z ，p ，k 分别是设计出的)(p G 的极点、零点及增益。

（5）求模拟低通原型滤波器)(p G 的分子分母系数，

[b,a]=zp2tf(z,p,k)。

（6）求出)(p G 的分子、分数系数。[B ，A]＝lp2lp （b ，a ，Wo ） （7）求出)(z H 的分子、分母系数，利用bilinear 函数。 （8）求频率响应)(ωj e H ，利用Freqz 函数。 （9）画)(ωj e H 的图 。

（四）实验报告：

1．根据实验原理，编写代码，得出实验结果，并画出波形图。 2．归纳、总结实验结果。 3．心得体会及其他。

实验四 用窗函数法设计FIR 滤波器

实验学时：2 实验类型：验证

实验要求：必修

（一）实验目的

1. 掌握窗函数法设计FIR 滤波器的原理和方法，观察用几种常用窗函数设计的FIR 数字滤波器技术指标；

2. 掌握FIR 滤波器的线性相位特性；

3. 了解各种窗函数对滤波特性的影响。

（二）实验原理

如果所希望的滤波器的理想频率响应函数为Hd(e jω)，则其对应的单位脉冲响应为ωπ

ωπ

π

ωd e e H n h n j j d ?-

=

)(21

)(，用窗函数wN(n)将hd(n)截断，并进行加权处理，

得到实际滤波器的单位脉冲响应h(n)=hd(n)wN(n)，其频率响应函数为

n j N n j e n h e H ωω

--=∑=1

0)()(。如果要求线性相位特性，则h(n)还必须满足

)1()(n N h n h --±=。可根据具体情况选择h(n)的长度及对称性。

（三）实验内容

题一：生成四种窗函数：矩形窗、三角窗、汉宁窗、海明窗，并观察其频率响应。

题二：根据下列技术指标，设计一个FIR 数字低通滤波器：wp=0.2π，ws=0.4π，ap=0.25dB ， as=50dB ，选择一个适当的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。

提示：根据窗函数最小阻带衰减的特性表，可采用海明窗可提供大于50dB的衰减，其过渡带为6.6π/N，因此具有较小的阶次。

（四）实验用到的MATLAB函数

可以调用MATLAB工具箱函数fir1实现本实验所要求的线性相位FIR-DF的设计，调用一维快速傅立叶变换函数fft来计算滤波器的频率响应函数。

fir1是用窗函数法设计线性相位FIRDF的工具箱函数，调用格式如下：

hn=fir1(N, wc, ‘ftype’, window)

fir1实现线性相位FIR滤波器的标准窗函数法设计。

hn=fir1(N,wc)可得到6 dB截止频率为wc的N阶(单位脉冲响应h(n)长度为N+1)FIR 低通滤波器，默认(缺省参数windows)选用hammiing窗。其单位脉冲响应h(n)满足线性相位条件：h(n)=h(N-1-n)

其中wc为对π归一化的数字频率，0≤wc≤1。

当wc=［wc1, wc2］时，得到的是带通滤波器。

hn=fir1(N,wc,’ftype’)可设计高通和带阻滤波器。

当ftype=high时，设计高通FIR滤波器；

当ftype=stop时，设计带阻FIR滤波器。

应当注意，在设计高通和带阻滤波器时，阶数N只能取偶数(h(n)长度N+1为奇数)。不过，当用户将N设置为奇数时，fir1会自动对N加1。

hn=fir1(N,wc,window)可以指定窗函数向量window。如果缺省window参数，则fir1默认为hamming窗。可用的其他窗函数有Boxcar, Hanning, Bartlett, Blackman, Kaiser和Chebwin窗。例如：

hn=fir1(N,wc,bartlett(N+1))使用Bartlett窗设计；

hn=fir1(N,wc,chebwin(N+1,R))使用Chebyshev窗设计。

hn=fir1(N,wc,’ftype’,window)通过选择wc、ftype和window参数(含义同上)，可以设计各种加窗滤波器。

（五）实验报告：

1．根据实验原理，编写代码，得出实验结果，并画出波形图。

2．归纳、总结实验结果。

3．心得体会及其他。

附录一：

1、fft.m

调用格式是

X= fft(x) 或X=fft(x,N)

对前者，若x的长度是2的整数次幂，则按该长度实现x的快速变换，否则，实现的是慢速的非2的整数次幂的变换；对后者，N应为2的整数次幂，若x 的长度小于N，则补零，若超过N，则舍弃N以后的数据。Ifft的调用格式与之相同。

2、randn

本文件可用来产生均值为零、方差为1、服从高斯（正态）分步的白噪声信号u(n)，其调用格式是

u=randn(N) 或randn=(M，N)

前者表示u为N维向量，后者表示u为M*N的矩阵。

3、Buttord.m

本文件用来确定数字低通或模拟低通滤波器的阶次，其调用格式分别是

1）[N，Wn]=buttord(Wp，Ws，Rp，Rs)

2）[N，Wn]=buttord(Wp，Ws，Rp，Rs，’s’)

格式1)对应数字滤波器，式中Wp,Ws分别是通带和阻带的截止频率，实际上它们是归一化频率，其值在0～1之间，1对应抽样频率的一半。对低通和高通滤波器，Wp,Ws都是标量，对带通和带阻滤波器，Wp,Ws都是1×2的向量。Rp,Rs分别是通带和阻带的衰减，单位为dB。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。格式2）对应模拟滤波器，式中各个变量的含义和格式1)相同，但Wp,Ws及Wn的单位为rad/s ，因此，它们实际上式频率 。

4、Buttap.m

本文件用来设计模拟低通原型滤波器)(p G ，其调用格式是

[z ，p ，k]=buttap(N)

N 是欲设计的低通原型滤波器的阶次，z ，p ，k 分别是设计出的)(p G 的极点、零点及增益。 5、Lp2lp.m 6、Lp2hp.m 7、Lp2bp.m 8、Lp2bs.m

从文件名可以看出，上述4个文件的功能分别是将模拟低通原型滤波器)(p G 转换为实际的低通、高通、带通及带阻滤波器，其调用格式分别为： （1）[B ，A]＝lp2lp （b ，a ，Wo ） 或 [B ，A]＝lp2hp （b ，a ，Wo ） （2）[B ，A]＝lp2bp （b ，a ，Wo ）

或 [B ，A]＝lp2bs （b ，a ，Wo ）

式中b ，a 分别是模拟低通原型滤波器)(p G 的分子、分母多项式的系数向量，B ，A 分别是转换后的)(s H 的分子、分母多项式的系数向量；在（1）中，Wo 是低通或高通滤波器的截止频率；在（2）中，Wo 是带通或带阻滤波器的中心频率，Bw 是其带宽。 9、Bilinear.m

本文件实现双线性变换，即由模拟滤波器)(s H 得到数字滤波器)(z H ，而s 和z 的关系由1

1

2+-=

z z T s s 给出。其调用格式是：

[Bz ，Az]＝ bilinear （B ，A ，Fs ）

式中B ，A 分别是)(s H 的分子、分母多项式的系数向量，Bz ，Az 分别是)(z H 的分子、分母多项式的系数向量， Fs 是抽样频率。 10、

Freqz.m

本文件用来在已知)(z B ,)(z A 的情况下求出系统的频率响应)(ωj e H ，格式是

[H ，w]＝freqz （b ，a ，N ，’whole’，Fs ）