数字信号处理实验指导书二

实验五 用FFT作谱分析

１、实验目的

（１）进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的性质）

（２）熟悉FFT算法原理及子程序的应用。

（３）掌握用FFT对连续信号和时域离散信号进行频谱分析的基本方法。了解可能出现的分析误差和原因，以便在实际中正确应用FFT。

2、实验仪器：

PC机一台 MATLAB软件

3、实验原理

如果用FFT对模拟信号进行谱分析，首先要把模拟信号转换成数字信号，转换时要求知道模拟信号的最高截止频率，以便选择满足采样定理的采样频率。一般选择采样频率是模拟信号中最高频率的3～4倍。另外要选择对模拟信号的观测时间，如果采样频率和观测时间确定，则采样点数也确定了。这里观测时间和对模拟信号进行谱分析的分辨率有关，最小的观测时间和分辨率成倒数关系。要求选择的采样点数和观测时间大于它的最小值。

用FFT作谱分析时，要求做FFT的点数服从２的整数幂，这一点在上面选择采样点数时可以考虑满足，即使满足不了，可以通过在序列尾部加０完成。

1

如果要进行谱分析的模拟信号是周期信号，最好选择观测时间是信号周期的整数倍。如果不知道信号的周期，要尽量选择观测时间长一些，以减少截断效应的影响。

用FFT对模拟信号作谱分析是一种近似的谱分析。首先一般模拟信号（除周期信号外）的频谱是连续频谱，而用FFT作谱分析得到的是数字谱，因此应该取FFT的点数多一些，用它的包络作为模拟信号的近似谱。另外，如果模拟信号不是严格的带限信号，会因为频谱混叠现象引起谱分析的误差，这种情况下可以预先将模拟信号进行预滤，或者尽量将采样频率取高一些。

一般频率混叠发生在折叠频率附近，分析时要注意因频率混叠引起的误差。最后要注意一般模拟信号是无限长的，分析时要截断，截断的长度和分辨率有关，但也要尽量取长一些，取得太短因截断引起的误差会很大。举一个极端的例子，一个周期性正弦波，如果所取观察时间太短，例如取小于一个周期，它的波形和正弦波相差太大，肯定误差很大，但如果取得长一些，即使不是周期的整倍数，这种截断效应也会小一些。

4、实验步骤及内容

（１）复习DFT的定义、性质和用DFT作谱分析的有关内容。

（２）复习FFT算法原理与编程思想。

（３）编制信号产生程序，产生以下典型信号供谱分析用：

x1(n)R4(n)

n1,0n3x2(n)8n,4n70, 其他n

2

4n,0n3x3(n)n3,4n70, 其他n

x4(n)cosn4

x5(t)cos(8t)cos(16t)cos(20t)

（４）分别以变换区间N8，16，32，对x1(n)进行FFT，画出相应的幅频特性曲线。

（５）分别以变换区间N8，16，对x2(n)，x3(n)进行FFT，画出相应的幅频特性曲线。

（６）分别以变换区间N4，8，16，对x4(n)进行FFT，画出相应的幅频特性曲线。

（７）对模拟信号x5(t)选择采样频率和采样点数。

对x5(t)cos(8t)cos(16t)cos(20t)，选择采样频率fsHz，采样点数N分别为16，32，。将模拟信号x5(t)转换成序列，用x5(n)表示，再分别对它们进行N点FFT，并画出相应的幅频特性曲线。 也可以参考例程得到另一种曲线结果。

5、实验用MATLAB函数介绍

fft(); figure(); plot(); stem(); abs();title(); xlabel(); ylabel(); text();

hold on; axis(); grid on; subplot(); sin(); cos(); 等。

6、思考题

3

（１）在N=8时，x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢?

（２）如果周期信号的周期预先不知道，如何用FFT进行谱分析?

7、实验报告要求

（１）简述实验目的及实验原理。

（２）编程实现各实验内容，列出实验清单及说明。

（３）将实验结果和理论分析结果进行比较，分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。并总结实验所得的主要结论。

（４）简要回答思考题。

例程参考

用DFT对连续信号作谱分析。已知

xa(t)=cos(200*pi*t)+sin(100*pi*t)+cos(50*pi*t);

选取不同的截取长度Tp，观察用DFT进行频谱分析时存在的截取效应（频谱泄漏和谱间干扰）。在计算机上用DFT对模拟信号进行谱分析时，只能以有限大的采样频率fs对模拟信号采样。对有限点样本序列（等价于截取模拟信号一段进行采样）作DFT变换得到模拟信号的近似频谱

clear;close all;

4

fs=400;T=1/fs;

Tp=0.04;N=Tp*fs;

N1=[N,4*N,8*N];%三种长度0.04s 4*0.04s 8*0.04s

%矩形窗截断

for m=1:3

n=1:N1(m);

xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T);

Xk=fft(xn,4096);

fk=fs*[0:4095]/4096;

subplot(3,2,2*m-1);

plot(fk,abs(Xk)/max(abs(Xk)));

if m==1 title('矩形窗截断');

end

end

5

%加海明窗截断

for m=1:3

n=1:N1(m);

wn=hamming(N1(m));

xn=(cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T)).*wn';

Xk=fft(xn,4096);

fk=fs*[0:4095]/4096;

subplot(3,2,2*m)

plot(fk,abs(Xk)/max(abs(Xk)));

if m==1 title('hamming窗截断');

end

6

实验六 IIR滤波器的设计与信号滤波

１、实验目的

（１）熟悉用双线性变换法设计IIR数字滤波器的原理与方法。

（２）掌握数字滤波器的计算机仿真方法。

（３）通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。

2、实验仪器：

PC机一台 MATLAB软件

3、实验原理

利用双线性变换设计IIR滤波器（只介绍巴特沃斯数字低通滤波器的设计），首先要设计出满足指标要求的模拟滤波器的传递函数Ha(s)，然后由Ha(s)通过双线性变换可得所要设计的IIR滤波器的系统函数

H(z)。

如果给定的指标为数字滤波器的指标，则首先要转换成模拟滤波器的技术指标，这里主要是边界频率

wp和ws的转换，对

p和s指标不作变化。边界频率的转换关系为

21tan(w)T2。接着，按照模拟低通

滤波器的技术指标根据相应设计公式求出滤波器的阶数N和3dB截止频率c；根据阶数N查巴特沃斯归一化低通滤波器参数表，得到归一化传输函数Ha(p)；最后，将

psc代入Ha(p)去归一，得到实际的模

21z1s1Ha(s)拟滤波器传输函数。之后，通过双线性变换法转换公式T1z，得到所要设计的IIR滤波器的

7

系统函数H(z)。

利用所设计的数字滤波器对实际的心电图采样信号进行数字滤波器。

4、实验步骤及内容

（１）复习有关巴特沃斯模拟滤波器的设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计一个巴特沃斯IIR低通数字滤波器。设计指标参数为：在通带内频率低于0.2时，最大衰减小于1dB；在阻带内0.3,频率区间上，最小衰减大于15dB。

（２）绘制出数字滤波器的幅频响应特性曲线。

（３）用所设计的滤波器对实际心电图信号采样序列（实验数据在后面给出）进行仿真滤波处理，并分别绘制出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

（４）输入为20Hz正弦和200Hz的正弦的叠加波形，要求用双线性变换法设计一巴特沃斯数字低通滤波器滤除200Hz的正弦，使输出中只保留20Hz的正弦波。并绘制出滤波前和滤波后的波形。

5、实验用MATLAB函数介绍

buttord(); butter(); bilinear(); freqz(); freqs(); filter(); figure(); plot(); stem(); abs();title(); xlabel(); ylabel(); text(); hold on; axis(); grid on; subplot();等

6、思考题

21z1s1（１）用双线性变换法设计数字滤波器过程中，变换公式T1z 中T的取值， 对设计结果有

8

无影响? 为什么?

（２）如果用脉冲响应不变法设计该IIR数字低通滤波器，程序如何改动？

7、实验报告要求

（１）简述实验目的及实验原理。

（２）编程实现各实验内容，列出实验清单及说明。

（３）由绘制的

H(ejw)特性曲线及设计过程简述双线性变换法的特点。

（４）对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。

（５）简要回答思考题。

8、心电图信号采样序列x(n)

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出的数据是一实际心电图信号采样序列样本x(n)，其中存在高频干扰。本实验中，以x(n)作为输入序列，滤除其中的干扰成分。

x(n)[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]

9

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

１、实验目的

（１）掌握用窗函数法设计FIR数字滤波器的原理和方法。

（２）熟悉线性相位FIR数字滤波器特性。

（３）了解各种窗函数对滤波特性的影响。

2、实验仪器：

PC机一台 MATLAB软件

3、实验原理

如果所希望的滤波器的理想频率响应函数为

1hd(n)2Hd(ejw)，则其对应的单位脉冲响应为

Hd(ejw)ejwndw

窗函数设计法的基本原理是用有限长单位脉冲响应h(n)逼近hd(n)。由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数w(n)将hd(n)截断，并进行加权处理，得到：

h(n)hd(n)w(n)

h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数H(ejw)为

10

H(e)h(n)ejwnjwn0N1

式中，Ｎ为所选窗函数w(n)的长度。

这种对理想单位取样响应的加窗处理对滤波器的频率响应会产生以下三点影响：

（1）使理想特性不连续的边沿加宽，形成一过渡带，过渡带的宽度取决于窗函数频谱的主瓣宽度。

（2）在过渡带两旁产生肩峰和余振，它们取决于窗函数频谱的旁瓣；旁瓣越多，余振也越多；旁瓣相对值越大，肩峰则越强。

（3）增加截断长度N，只能缩小窗函数频谱的主瓣宽度而不能改变旁瓣的相对值；旁瓣与主瓣的相对关系只决定于窗函数的形状。因此增加N，只能相对应减小过渡带宽。而不能改变肩峰值。肩峰值的大小直接决定通带内的平稳和阻带的衰减，对滤波器性能有很大关系。例如矩形窗的情况下，肩峰达8.95%，致使阻带最小衰减只有21分贝，这在工程上往往是不够的。怎样才能改善阻带的衰减特性呢？只能从改善窗函数的形状上找出路，所以希望的窗函数频谱中应该减少旁瓣，使能量集中在主瓣，这样可以减少肩峰和余振，提高阻带衰减。而且要求主瓣宽度尽量窄，以获得较陡的过渡带，然而这两个要求总不能同时兼得，往往需要用增加主瓣宽度带换取较大的阻带衰急，于是提出了海明窗、汉宁窗、布莱克曼窗、凯塞窗、切比雪夫窗等窗函数。

所以，用窗函数法设计的滤波器性能取决于窗函数w(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度。设待求滤波器的过渡带用w表示，它近似等于窗函数主瓣宽度。因过渡带w近似与窗口长度成反比，NA/w，A决定于窗口形式。例如，矩形窗A=4π，海明窗A=8π等。按照过渡带及阻带衰减情况，选择窗函数形式。原则是在保证阻带衰减满足要求的情况下， 尽量选择主瓣窄的窗函数。

11

jwjw这样选定窗函数类型和窗口长度N后，求出单位脉冲响应h(n)hd(n)w(n)，再求出H(e)。H(e)是jwH(e)。h(n)否满足要求，要进行验算。一般在的尾部加零使长度满足２的整数次幂，以便用FFT计算jwH(e)不满足要求，则要重新选择窗函数类型和长度N，再如果要观察细节，补零点数增多即可。如果

次验算，直至满足要求。

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

h(n)h(N1n)

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性，可选择h(n)h(N1n)一类， 而不能选

h(n)h(N1n)一类。

4、实验步骤及内容

wp0.2（１）根据下列技术指标，设计一个线性相位的FIR数字低通滤波器。通带截止频率带允许波动

Ap0.25dB，通

；阻带截止频率wS0.4，阻带衰减AS50dB。

（2）调用fir1()函数得到所设计的低通滤波器的单位脉冲响应，调用fft()函数进行频响验证。打印输出各部分结果。

（3）编程验证窗长和窗形状对实际滤波器性能的影响。如要求用窗函数法设计一个线性相位FIR数字低通滤波器，用理想低通滤波器作为逼近滤波器，截止频率

wc4rad，用四种窗函数（矩形窗，汉宁

窗（升余弦窗），哈明窗（改进的升余弦窗），布莱克曼窗）设计该滤波器，选择窗函数的长度N15,33两种情况。

12

（4）输入为20Hz正弦和200Hz的正弦的叠加波形，要求用窗函数法设计一数字低通滤波器滤除200Hz的正弦，使输出中只保留20Hz的正弦波。并绘制出滤波前和滤波后的波形。

5、实验用MATLAB函数介绍

fir1(); fft(); freqz(); boxcar(); hamming(); hanning(); blackman(); sin();figure(); plot(); stem(); abs();angle(); title(); xlabel(); ylabel(); text(); hold on; axis(); grid on; subplot();等

6、思考题

如果要求用窗函数法设计带通滤波器， 且给定上、 下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应hd(n)。

7、实验报告要求

（１）简述实验目的及实验原理。

（２）编程实现各实验内容，列出实验清单及说明。

（３）总结窗函数法设计FIR数字滤波器的步骤。

（４）简要回答思考题。

13

附 录

一、常用的MATLAB函数介绍

１． plot

功能：线型绘图函数。

格式：plot（v）；plot（x,y）

说明：

plot（v），v是长度为N的数值向量。作用是在坐标系中顺序地用直线连接顶点，生成一条折（曲）线。当向量元素充分多时，即可生成一条光滑的曲线。在实验中，若FFT点数足够多时，用plot打印的

jw|X(e)|连续曲线。 幅频特性就很接近

plot（x,y），x和y都是长度为n的向量。作用在坐标系中生成顺序连接顶点的折（曲）线。这种调用可被用来生成参数方程的图形。

２．stem

功能：绘制离散序列图。

格式：stem(y)； stem(x,y) ；stem(…,’线端符号’) ；stem(…,’线型’)；

stem(…,‘线型’，‘线端符号’)

14

说明：

stem(y)和stem(x,y)分别与plot（v）和plot（x,y）的绘图规则相同，只是stem绘制的是离散序列图。实验中用于绘制时域序列x(n)的波形图和序列的离散傅里叶变换X(k)的幅度图。

３．subplot

功能：多坐标设置与定位当前坐标系。

格式：subplot(m,n,k)

说明：

subplot(m,n,k) 将图形窗口分成m行n列的m×n块子区域，按行从上到下，从左到右的顺序，在第i块子区定义一个坐标系，使其成为当前坐标系，随后的绘图函数将在该坐标系输出图形。另外，同一个图形窗口的坐标系可以重叠，这样可以产生前面的坐标系遮住后面坐标系的各种图形效果。

４．figure

功能：创建新的图形窗口（用于输出图形的窗口）。

格式：figure；figure(h)

说明：

figure 函数创建一个新的图形窗口，并成为当前图形窗口，所创建的图形窗口的序号（句柄值）是按同一MATLAB程序中创建的顺序号。使用figure(h)函数，该方法常用在程序设计中，用于控制将

15

各种波形图输出到相应的图形窗口中。打印输出或存储时，一个图形窗口打印一张图纸或存储一个图形文件。

５．axis

功能：设置图形的坐标范围。

格式：axis([xmin xmax ymin ymax])６．grid on

功能：画出图形的分格线。

7．title

功能： 书写图名。

格式：title（‘s’）

8．xlabel

功能： 横坐标名称。

格式：xlabel（‘s’）

9．ylabel

16

功能： 纵坐标名称。

格式：ylabel（‘s’）

10．text

功能：在图面（x,y）位置处书写字符注释。

格式：text（x,y,‘s’）

11．hold on

功能：通过该命令，保持当前图形不变化，准备在当前的图形窗口上绘制新的曲线。

12．abs

功能：求绝对值（模值）。

格式：y=abs(x)

说明：

y=abs(x)用于计算x的绝对值，当x为复数时，得到的是复数的模值。当x为字符串时，abs(x)得到字符串的各个字符的ASCII码，例如，x=’123’,则，abs(x)得到：49 50 51。

13．angle

17

功能：求相角。

格式：ψ＝angle(h)

说明：

ψ＝angle(h)用于求复矢量或复矩阵的相角（以弧度为单位），相角介于和之间。

14．conv

功能：求卷积。

格式：c=conv(a,b)

说明：conv(a,b)用于求矢量a和b的卷积，即

c(n)a(k1)b(nk),n1,2,k0N1

式中Ｎ为矢量a和b的最大长度。例如，当a=[1 2 3],b=[4 5 6]时，则

c=conv(a,b)

c=

4 13 28 27 18

此函数可直接用于求两个有限长序列的卷积。设x(n)和h(n)的长度分别为M和N，则计算二者卷

18

积的MATLAB语句如下：

y=conv(x,h)

y的长度为M+N-1。

15．filter

功能：利用IIR滤波器或FIR滤波器对数据进行滤波。

格式：y=filter(b,a,x); [y,zf]=filter(b,a,x); y=filter(b,a,x,zi)

说明：

filter利用数字滤波器对数据进行滤波，其实现采用直接Ⅱ型结构，因而适用于IIR和FIR两种滤波器。滤波器的系统函数为

b0b1z1bMzMH(z)1a1z1aNzN

即滤波器系数a=[a0 a1 a2 …aN],b=[b0 b1 … bM]。

这里的标准形式为a0=1，如果输入矢量a时，a0≠1，则MATLAB将自动进行归一化系数的操作；如果a0=０，则给出出错信息。

y=filter(b,a,x)利用给定系数矢量a和b对x中的数据进行滤波，结果放入y矢量中。

y=filter(b,a,x,zi)可在zi中指定x的初始状态。

19

[y,zf]=filter(b,a,x)除得到矢量y外，还得到x的最终状态矢量zf。

16．fftfilt

功能：利用重叠相加法实现短序列b(n)和长序列x(n)的块之间的卷积。

格式：y=fftfilt(b,x)

17．freqz

功能：数字滤波器的频率响应。

格式：［H,w］=freqz(b,a,Ｎ)；［Ｈ，f］=freqz(b,a,N,Fs);

H=freqz(b,a,w); H=freqz(b,a,f,Fs); freqz(b,a)

说明：

jwH(e)。 freqz用于计算数字滤波器H(z)的频率响应函数

［H,w］=freqz(b,a,Ｎ)可得到数字滤波器的Ｎ点频率响应值，这Ｎ个点均匀地分布在［０，］上，并将这Ｎ个频点的频率记录在w中，相应的频响值记录在Ｈ中。要求Ｎ为大于零的整数，最好为２的整数次幂，以便采用FFT计算，以提高速度。

jwH(e)在［０，Fs/2］上等间隔采样Ｎ点，采样点频率及相应频［Ｈ，f］=freqz(b,a,N,Fs)用于对

响值分别记录在f和H中。由用户指定Fs (以Hz为单位)值。

20

jwH=freqz(b,a,w)用于对H(e)在［０， 2］上进行采样，采样频率点由矢量w指定。

jwH(e)在［０， Fs］上采样，采样频率点由矢量f给定。 H=freqz(b,a,f,Fs) 用于对

freqz(b,a)用于在当前窗口中绘制出幅频和相频特性曲线。

18．impz

功能：计算H(z)相应的单位脉冲响应h(n)。

格式：［h,t］=impz(b,a); [h,t]=impz(b,a,z);

19．fft

功能：一维快速傅里叶变换(FFT)

格式：y=fft(x); y=fft(x,N)

说明：

fft函数用于计算矢量或矩阵的离散傅里叶变换。

y=fft(x)利用FFT算法计算矢量x的离散傅里叶变换，当x为矩阵时，y为矩阵x每一列的FFT。当x的长度为2的整数次幂时，fft采用基2FFT算法，否则采用稍慢的混合基算法。

y=fft(x,N)采用Ｎ点FFT。当x长度小于Ｎ时，fft函数自动在x尾部补零，以构成Ｎ点数据；当x的长度大于Ｎ时，fft截取x的前面Ｎ点数据进行FFT。

21

20．ifft

功能：一维逆快速傅里叶变换（IFFT）

格式：y=ifft(x); y=ifft(x,N)

二、MATLAB在信号处理中的应用举例

1、线性卷积与循环卷积的计算

对于无限长序列不能用MATLAB直接计算线性卷积，在MATLAB内部只提供了一个conv函数计算两个有限长序列的线性卷积。对于循环卷积MATLAB内部没有提供现成的函数，我们可以按照定义式直接编程计算。

例１： 已知两序列：

0.8nx(n)01h(n)00n11其它0n5其它

求它们的线性卷积yl(n)h(n)x(n)和N点的循环卷积yc(n)h(n)x(n)，并研究两者之间的关系。

MATLAB实现程序：

循环卷积的函数

function yc=circonv(x1,x2,N)

22

%realize circular convolution use direct method

%y=circonv(x1,x2,N)

%y:output sequences

%x1,x2:input sequences

%N:circulation length

if length(x1)>N

error('N must not be less than length of x1');

end

if length(x2)>N

error('N must not be less than length of x2');

end

%以上语句判断两个序列的长度是否小于N

x1=[x1,zeros(1,N-length(x1))];

%填充序列x1(n)使其长度为N1+N2-1(序列%h(n)的长度为N1，序列x(n)的长度为N2)

23

x2=[x2,zeros(1,N-length(x2))];

%填充序列x2(n)使其长度为N1+N2-1

n=[0:1:N-1];

x2=x2(mod(-n,N)+1); %生成序列x2((-n))N

H=zeros(N,N);

for n=1:1:N

H(n,:)=cirshiftd(x2,n-1,N); %该矩阵的k行为x2((k-1-n))N

end

yc=x1*H'; %计算循环卷积

function y=cirshiftd(x,m,N)

%directly realize circular shift for sequence x

%y=cirshiftd(x,m,N);

%x:input sequence whose length is less than N

%m:how much to shift

24

%N:circular length

%y:output shifted sequence

if length(x)>N

error('length of x must be less than N');

end

x=[x,zeros(1,N-length(x))];

n=[0:1:N-1];

y=x(mod(n-m,N)+1);

研究两者之间的关系

clear all;

n=[0:1:11];

m=[0:1:5];

N1=length(n);

N2=length(m);

25

xn=0.8.^n; %生成x(n)

hn=ones(1,N2); %生成h(n)

yln=conv(xn,hn); %直接用函数conv计算线性卷积

ycn=circonv(xn,hn,N1); ny1=[0:1:length(yln)-1];

ny2=[0:1:length(ycn)-1];

subplot(2,1,1); %画图stem(ny1,yln);

subplot(2,1,2);

stem(ny2,ycn);

axis([0,16,0,4]);

运行结果：

%用函数circonv计算N1点循环卷积

26

43线性卷积2100246810121413循环卷积2100246810121416

线性卷积和循环卷积的比较图

2、利用离散傅里叶变换（DFT）分析信号的频谱

MATLAB中计算序列的离散傅里叶变换和逆变换是采用快速算法，利用fft和ifft函数实现。函数fft用来求序列的DFT，函数ifft用来求IDFT。

例２：已知序列x(n)2sin(0.48n)cos(0.52n) 0n100，试绘制x(n)及它的离散傅里叶变换|X(k)|图。

MATLAB实现程序：

clear all

N=100;

n=0:N-1;

27

xn=2*sin(0.48*pi*n)+cos(0.52*pi*n);

XK=fft(xn,N);

magXK=abs(XK);

phaXK=angle(XK);

subplot(1,2,1)

plot(n,xn)

xlabel('n');ylabel('x(n)');

title('x(n) N=100');

subplot(1,2,2)

k=0:length(magXK)-1;

stem(k,magXK,'.');

xlabel('k');ylabel('|X(k)|');

title('X(k) N=100');

运行结果：

28

x(n) N=100321600-1-2-320100X(k) N=10080|X(k)|400x(n)050n100050k100

信号及离散傅里叶变换|X(k)|图

3、利用FFT实现线性卷积

若序列x1(n)、x2(n)为长度分别为N1、N2的有限长序列，两序列的Ｎ点的循环卷积yc(n)x1(n)x2(n)，

yl(n)x1(n)x2(n)。由DFT的性质可知：当NN1N21时有yl(n)yc(n)IDFT{DFT[x1(n)]DFT[x2(n)]}。

序列较长时DFT运算通常用快速算法FFT实现。在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅里叶变换和逆变换。

例３：用FFT实现例1中两序列的线性卷积。

MATLAB实现程序：

n=[0:1:11];

m=[0:1:5];

N1=length(n);

N2=length(m);

29

xn=0.8.^n; %生成x(n)

hn=ones(1,N2); %生成h(n)

N=N1+N2-1;

XK=fft(xn,N);

HK=fft(hn,N);

YK=XK.*HK;

yn=ifft(YK,N);

if all(imag(xn)==0)&(all(imag(hn)==0)) %实序列的循环卷积仍然为实序列

yn=real(yn)

end

x=0:N-1;

stem(x,yn,'.')运行结果：

30

4321005101520

利用傅里叶变换实现线性卷积的结果图

4、IIR滤波器的设计与实现

基于模拟滤波器变换原理IIR滤波器的经典设计：首先，根据模拟滤波器的指标设计出相应的模拟滤波器，然后，将设计好的模拟滤波器转换成满足给定技术指标的数字滤波器。通常算法有脉冲响应不变法和双线性变换法。在MATLAB的数字信号处理工具箱中提供了相应的设计函数。常用的有：

（１）Butterworth滤波器阶数选择函数

[N,Wn]=buttord(Wp,Ws,Rp,Rs)

输入参数：Wp通带截止频率，Ws阻带截止频率，Rp通带最大衰减，Rs阻带最小衰减；

输出参数：N符合要求的滤波器最小阶数，Wn为Butterworth滤波器固有频率（3dB）。

（２）零极点增益模型到传递函数模型的转换

[num,den]=zp2tf(Z,P,K);

31

输入参数：Z,P,K分别表示零极点增益模型的零点、极点和增益；

输出参数：num,den分别为传递函数分子和分母的多项式系数。

（３）从低通向低通的转换

[b,a]=lp2lp(Bap,Aap,Wn);

功能：把模拟滤波器原型转换成截至频率为Wn的低通滤波器。

（４）双线性变换函数

[bz,az]=bilinear(b,a,Fs);

功能：把模拟滤波器的零极点模型转换为数字滤波器的零极点模型，其中Fs是采样频率。例4：用双线性变换法设计一个Butterworth数字低通滤波器，要求其通带截至频率0.2　带截至频率0.4　rad，通带衰减

Rp小于0.75dB，阻带衰减Rs大于20dB。

MATLAB实现程序：

wp=0.2*pi; ws=0.4*pi;

ap=0.75; as=20;

T=1;

32

rad，阻

Wp=2/T*tan(wp/2); %求模拟原型滤波器参数

Ws=2/T*tan(ws/2);

[N,Wc]=buttord(Wp,Ws,ap,as,'s'); %求滤波器的阶数

[BS,AS]=butter(N,Wc,'s');%得到模拟原型滤波器的系统函数

[BZ,AZ]=bilinear(BS,AS,1/T);%双线性变换法转换到数字滤波器的系统函数

[H,W]=freqs(BS,AS);

plot(W/pi,20*log10(abs(H)));%画模拟原型滤波器的频响

hold on;grid;

[H,W]=freqz(BZ,AZ);

plot(W/pi,20*log10(abs(H)),'r'); %画数字滤波器的频响

运行结果：

33

Butterworth低通滤波器的频率响应图

5、FIR滤波器的设计与实现

FIR滤波器的设计方法常用的有窗函数法和频率采样法两种，在MATLAB的数字信号处理工具箱中提供了函数fir1。fir1是采用经典窗函数法设计线性相位FIR数字滤波器，且具有标准低通、带通、高通和带阻等类型。

语法格式：

B=fir1(n,Wn)

B=fir1(n,Wn,’ftype’)

B=fir1(n,Wn,window)

B=fir1(n,Wn,’ftype’,window)

其中，n为FIR滤波器的阶数，对于高通、带阻滤波器n取偶数。Wn为滤波器截止频率，取值范

34

围0～1；对于带通、带阻滤波器，Wn=[W1,W2],且W1Matlab提供了几个函数来实现这些窗函数。

W=boxcar(N); Ｎ点的矩形窗

W=triang(N); Ｎ点的三角窗

W=hanning(N); Ｎ点的汉宁窗

W=hamming(N); Ｎ点的海明窗

W=blackman(N); Ｎ点的布莱克曼窗

W=kaiser(N,beta); Ｎ点的凯泽窗(beta=0.7865)

wp0.2例5： 用窗函数法设计一个线性相位FIR低通滤波器，性能指标：通带截止频率止频率ws0.3，阻带衰减不小于40dB，通带衰减不大于3dB。

，阻带截

MATLAB实现程序：

%程序中采用汉宁窗(hanning)

wp=0.2*pi;ws=0.3*pi;

35

wdelta=ws-wp;

N=ceil(8*pi/wdelta);

Wn=(0.2+0.3)*pi/2;

b=fir1(N,Wn/pi,hanning(N+1));

freqz(b,1,512)

运行结果为：

500Magnitude (dB)-50-100-150-20000.10.20.30.40.50.60.7Normalized Frequency ( rad/sample)0.80.910-500Phase (degrees)-1000-1500-2000-250000.10.20.30.40.50.60.7Normalized Frequency ( rad/sample)0.80.91

FIR滤波器的幅频和相频特性图

36