数字信号处理实验指导书（审）

电 工 电 子 实 验 中 心

实 验 指 导 书

数字信号处理 实验教程

二○○九 年 三 月

高等学校电工电子实验系列

数字信号处理实验教程

主编 石海霞 周玉荣

攀枝花学院电气信息工程学院

电工电子实验中心

内容简介

内 容 简 介

数字信号处理是一门理论与实践紧密联系的课程，适当的上机实验有助于深入理解和巩固验证基本理论知识，了解并体会数字信号处理的CAD手段和方法，锻炼初学者用计算机和MATLAB语言及其工具箱函数解决数字信号处理算法的仿真和滤波器设计问题的能力。

本实验指导书结合数字信号处理的基本理论和基本内容设计了八个上机实验，每个实验对应一个主题内容，包括常见离散信号的MATLAB产生和图形显示、离散时间系统的时域分析、离散时间信号的DTFT、离散时间信号的Z变换、离散傅立叶变换DFT、快速傅立叶变换FFT及其应用、基于MATLAB的IIR和FIR数字滤波器设计等。此外，在附录中，还简单介绍了MATLAB的基本用法。每个实验中，均给出了实验方法和步骤，还有部分的MATLAB程序，通过实验可以使学生掌握数字信号处理的基本原理和方法。

目录

目 录

绪 论 ............................................................................................................................ 1 实验一 常见离散信号的MATLAB产生和图形显示.................................................... 2 实验二 离散时间系统的时域分析 ............................................................................. 6 实验三 离散时间信号的DTFT .................................................................................... 9 实验四 离散时间信号的Z变换 ............................................................................... 14 实验五 离散傅立叶变换DFT .................................................................................... 18 实验六 快速傅立叶变换FFT及其应用 ................................................................... 24 实验七 基于MATLAB的IIR数字滤波器设计(设计实验) ...................................... 30 实验八 基于MATLAB的FIR数字滤波器设计(设计实验) ...................................... 33 附 录 .......................................................................................................................... 37 参 考 文 献 .............................................................................................................. 40

绪论

绪 论

随着电子技术迅速地向数字化发展，《数字信号处理》越来越成为广大理工科，特别是IT领域的学生和技术人员的必修内容。

数字信号处理是把信号用数字或符号表示成序列，通过计算机或通用（专用）信号处理设备，用数值计算方法进行各种处理，达到提取有用信息便于应用的目的。数字信号处理的理论和技术一出现就受到人们的极大关注，发展非常迅速。而且随着各种电子技术及计算机技术的飞速发展，数字信号处理的理论和技术还在不断丰富和完善，新的理论和技术层出不穷。目前数字信号处理已广泛地应用在语音、雷达、声纳、地震、图象、通信、控制、生物医学、遥感遥测、地质勘探、航空航天、故障检测、自动化仪表等领域。

数字信号处理是一门理论和实践、原理和应用结合紧密的课程，由于信号处理涉及大量的运算，可以说离开了计算机及相应的软件，就不可能解决任何稍微复杂的实际应用问题。Matlab是1984年美国Math Works公司的产品，MATLAB语言具备高效、可视化及推理能力强等特点，它的推出得到了各个领域专家学者的广泛关注，其强大的扩展功能为各个领域的应用提供了基础，是目前工程界流行最广的科学计算语言。早在20世纪90年代中期，MATLAB就己成为国际公认的信号处理的标准软件和开发平台。从1996年后，美国新出版的信号处理教材就没有一本是不用MATLAB的。

本实验指导书结合数字信号处理的基本理论和基本内容，用科学计算语言MATLAB实现数字信号处理的方法和实践， 通过实验用所学理论来分析解释程序的运行结果，进一步验证、理解和巩固学到的理论知识，从而达到掌握数字信号处理的基本原理和方法的目的。

1

实验一 常见离散信号的MATLAB产生和图形显示

实验一 常见离散信号的MATLAB产生和图形显示

一、实验目的

１. 学会用MATLAB在时域中产生一些基本的离散时间信号。 ２. 了解信号的各种运算。

二、实验原理

（一）、序列的产生

由于MATLAB数值计算的特点，用它来分析离散时间信号与系统是很方便的。离散信号是数字信号处理的最基础的内容，由于内存有限，MATLAB无法表示无限序列。在MATLAB中，可以用一个列向量来表示一个有限长度的序列，但是这种表示方法没有包含采样位置的信息，要完全表示x(n)，要用x 和n 两个向量，例如

x(n)={2,1,0,2,3,-1,2,3} ↑ 在MATLAB中表示为

n=[-4,-3,-2,-1,0,1,2,3]; x=[2,1,0,2,3,-1,2,3];

当序列从n=0开始，则不需要采样位置信息，这时可以只用x 来表示

1． 单位抽样序列

?(n)??

?0?1n?0n?0

在MATLAB中可以利用zeros()函数实现。

x?zeros(1,N);x(1)?1;

如果?(n)在时间轴上延迟了k个单位，得到?(n?k)即：

?(n?k)??

?0?1n?kn?0

2．单位阶跃序列

2

实验一 常见离散信号的MATLAB产生和图形显示

n?0?1 u(n)?

n?0?0在MATLAB中可以利用ones()函数实现。

x?ones(1,N)

3．正弦序列

x(n)?Asin(2?fn/Fs??)

在MATLAB中

n?0:N?1x?A*sin(2*pi*f*n/Fs?fai)

4．复正弦序列

x(n)?ej?n

在MATLAB中

n?0:N?1x?exp(j*w*n)

5．指数序列

x(n)?a

n在MATLAB中

n?0:N?1x?a.^n

（二）、简单运算 1．信号加

x(n)?x1(n)?x2(n) 在MATLAB中

x=x1+x2

注意：x1和x2应该具有相同的长度，位置对应，才能相加。否则，需要先通过Zeros函数左右补零后再相加。

2．信号延迟

3

实验一 常见离散信号的MATLAB产生和图形显示

给定离散信号x(n)，若信号y(n)为 y(n)= x(n-k)

那么y(n)就是信号x(n)在时间轴上右移k个采样周期后得到的新的序列。 在MATLAB中

y(n)= x(n-k) 3．信号乘

x(n)?x1(n)?x2(n)

在MATLAB中

x?x1.?x2

这是信号的点乘运算，所以同样需要信号加需要的x1和x2二者的长度要相等这一前提条件。

4．信号变化幅度

y(n)= k×x(n) 在MATLAB中

y?k?x

5．信号翻转

y(n)= x(-n)

在MATLAB中

y=fliplr(x)

三、实验内容与步骤

１. 产生一个单位样本序列x1(n)，起点为ns= -10, 终点为nf=20, 在n0=0时有一单位脉冲并显示它。修改程序，以产生带有延时11个样本的延迟单位样本序列x2(n)= x1(n-11)，并显示它。

２.已知 c= -(1/12)+(pi/6)*i;产生一个复数值的指数序列x2(n)=2*exp(c*n)，起点为ns= 0, 终点为nf=40;并显示它。

3．产生一个正弦序列x3(n) =1.5*cos(2*pi*f*n); 起点为ns= 0, 终点为nf=40;并显示它。

4．复杂信号的产生：复杂的信号可以通过在简单信号上执行基本的运算来产生 试产生一个振幅调制信号

y(n)?(1?m?cos(2?fLn))?cos(2?fHn)?(1?0.4?cos(2??0.01n))?cos(2??0.1n)

4

实验一 常见离散信号的MATLAB产生和图形显示

n=0:100

四、实验仪器设备

计算机，MATLAB软件

五、实验注意事项

预先阅读附录（MATLAB基础介绍）；

六、思考题

１. 讨论复指数序列x2(n)的哪个参数控制该序列的增长或衰减率？哪个参数控制该序列的振幅？

２. 讨论正弦序列x3(n)的哪个参数控制该序列的相位？哪个参数控制该序列的振幅？

3. 讨论算术运算符*和. *之间的区别是什么?

5

实验二 离散时间系统的时域分析

实验二 离散时间系统的时域分析

一、实验目的

１. 运用MATLAB仿真一些简单的离散时间系统，并研究它们的时域特性。 ２. 运用MATLAB中的卷积运算计算系统的输出序列，加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。

二、实验原理

离散时间系统其输入、输出关系可用以下差分方程描述：

N?k?0dky[n?k]?M?k?0pkx[n?k]

当输入信号为冲激信号时，系统的输出记为系统单位冲激响应

?[n]?h[n]，则系统响应为如下的卷积计算式：

?y[n]?x[n]?h[n]??x[m]h[n?m]

m??? 当h[n]是有限长度的（n：[0，M]）时，称系统为FIR系统；反之，称系统为IIR系统。在MATLAB中，可以用函数y=Filter(p,d,x) 求解差分方程，也可以用函数 y=Conv(x,h)计算卷积。 例2.1

clf; n=0:40; a=1;b=2; x1= 0.1*n; x2=sin(2*pi*n); x=a*x1+b*x2; num=[1, 0.5,3];

6

实验三 离散时间信号的DTFT

4．运行下面程序并显示它，验证离散时间傅立叶变换时域卷积性质。

clf;

w=-pi:2*pi/255: pi; x1=[1 3 5 7 9 11 13 15 17]; x2=[1 -2 3 -2 1]; y=conv(x1,x2); h1=freqz(x1,1,w); h2=freqz(x2,1,w); hp=hi.*h2; h3=freqz(y,1,w); subplot(2,2,1)

plot(w/pi,abs (hp));grid title(‘幅度谱的乘积’) subplot(2,2,2)

plot(w/pi,abs (h3));grid title(‘卷积后序列的幅度谱’) subplot(2,2,3)

plot(w/pi,angle (hp));grid title(‘相位谱的和’) subplot(2,2,4)

plot(w/pi,angle (h3));grid title(‘卷积后序列的相位谱’)

四、实验仪器设备

计算机，MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序；

六、思考题

１.讨论实验程序1中的离散时间系统的频率响应是离散的还是连续的，是否是周期的？周期为多少？y与

２.讨论实验程序2中h1和h2的关系是什么？哪个参数控制时移量？

12

实验三 离散时间信号的DTFT

3. 讨论实验程序3中h1和h2的关系是什么？哪个参数控制频移量？

4. 讨论实验程序4中y与x1和x2的关系是什么？h1和h2与x1和x2的关系是什么？h1和hp相等吗？

13

实验四 离散时间信号的Z变换

实验四 离散时间信号的Z变换

一、实验目的

１. 运用MATLAB理解Z变换及其绘制H(z)的零极点图。 ２. 运用MATLAB计算逆Z变换。

二、实验原理

（一）、MATLAB在ZT中的应用。

线性时不变离散时间系统的冲激响应h(n)的z变换是其系统函数H(z)， 在MATLAB中可以利用性质求解Z变换，例如可以利用线性卷积求的Z变换。若H(z)的收敛域包含单位圆，即系统为稳定系统，即系统在单位圆上z?ej?处计算的是系统的频率响应。

（二）、逆Z变换

Z变换对于分析和表示离散线性时不变系统具有重要作用。但是在MATLAB中不能直接计算Z变换，但是对于一些序列可以进行逆Z变换。

已知序列的Z变换及其收敛域， 求序列称为逆Z变换。 序列的Z变换及共逆Z变换表示如下：

?X(z)?x(n)??n???x(n)z?n,Rx??z?Rx?n?1

12?

通常，直接计算逆Z变换的方法有三种：围线积分法、长除法和部分分式展开法。在实际中，直接计算围线积分比较困难，往往不直接计算围线积分。由于序列的Z变换常为有理函数，因此采用部分分式展开法比较切合实际，它是将留数定律和常用序列的Z变换相结合的一种方法。

设x(n)的Z变换X(z)是有理函数，分母多项式是N阶，分子多项式是M阶，将X(z)展成一些简单的常用的部分分式之和，通过常用序列的Z变换求得各部分的逆变换，再相加即得到原序列x(n)。在MATLAB中提供了函数residuez来实现上述过程，调用格式如下：

14

?j?cX(z)zdz,c?(Rx?,Rx?)

实验四 离散时间信号的Z变换

[R，P，K]= residuez（B，A）

其中B、A分别是有理函数分子多项式的系数和分母多项式的系数，输出R是留数列向量，P是极点列向量。如果分子多项式的阶数大于分母多项式的阶数，则K返回为常数项的系数。 例4.1 计算 X(z)?的Z反变换。

由于分母多项式为：1?1.5z?1?0.5z?2，则 MATLAB实现： clear b=1; a=[1,-1.5,0.5]; [R,P,K]=residuez(b,a) 输出R =

2 -1 P =

1.0000 0.5000 K = []

因此：得到X(z)的部分分式展开为X(z)?得x(n)?(2?0.5)u(n)。

三、实验内容与步骤

1. 运行下面程序，利用线性卷积求Z变换。 设X1(z)?z?2?3z?1n1(1?z?1)(1?0.5z?1),z?1

21?z?1??11?0.5z?1，根据常用序列的z变换可

，X2(z)?2z?4z?3?5z2?1，求X3(z)?X1(z)?X2(z)

由变换定义可知：x1(n)?{1,2,3},n?{?1,0,1} x2(n)?{2,4,3,5},n?{?2,?1,0,1}

15

实验四 离散时间信号的Z变换

通过求x3(n)?x1(n)?x2(n)，再求其z变换得到X3(z)?X1(z)?X2(z)。 MATLAB程序： x1=[]1,2,3]; x2=[2,4,3,5]; n1=-1:1; n2=-2:1; x3=conv(x1,x2) nb3=n1(1)+n2(1);

nc3=n1(length(x1))+n2(length(x2)); n3=[nb3:nc3]

２.已知两个线性时不变的因果系统，系统函数分别为

H1(z)?1?z?N ，

H2(z)?1?zN?N?N1?az

分别令N=8，a=0.8，运行下面程序计算并显示示这两个系统的零、极点图及幅频特性。 程序：

b=[1,0,0,0,0,0,0,0,-1]; %H1(z)和H2(z)的分子多项式系数向量 a0=1; %H1(z)分母多项式系数向量

a1=[1,0,0,0,0,0,0,0,-(0.8)^8]; % H2(z)的分母多项式系数向量 [H, w]=freqz(b,a0); [H1, w1]=freqz(b,a1); subplot(2,2,1);zplane(b,a0);

xlabel(‘实部’);ylabel(‘虚部’); title(‘H1 (z)系统的零极点图’)； subplot(2,2,2);zplane(b,a1);

16

实验四 离散时间信号的Z变换

xlabel(‘实部’);ylabel(‘虚部’); title(‘H2(z)系统的零极点图’)； subplot(2,2,3)；plot(w/pi,abs(H)); title(‘H1 (z)系统幅频响应曲线’)； subplot(2,2,4)；plot(w/pi,abs(H1)); title(‘H2(z)系统幅频响应曲线’)；

3.编写MATLAB程序 计算系统

X(z)?1(1?0.3z?1?1?1?1)(1?0.4z)(1?0.5z)(1?0.8z),z?0.8的z反变换。

四、实验仪器设备

计算机，MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序；

六、思考题

１.讨论实验程序1中的n3代表什么含义？根据x3和n3的结果写出

X3(z)?X1(z)?X2(z)的结果。

２.讨论实验程序2中极点与系统稳定性的关系？根据程序运行结果判断该系统的稳定性。 3. 根据实验程序3的运行结果写出z反变换x(n)。

17

实验五 离散傅立叶变换DFT

实验五 离散傅立叶变换DFT

一、实验目的

１. 运用MATLAB计算有限长序列的DFT和IDFT。 ２. 运用MATLAB验证离散傅立叶变换的性质。 3 .运用MATLAB计算有限长序列的圆周卷积。

二、实验原理

（一）、离散傅立叶变换DFT的定义

一个有限长度的序列x(n)（0≤n

（0???2?）上对X(ej?)均匀采样得到

?X(k)?X(ej?)

??2?k/N??n???x(n)e?j2?kn/N 0?k?N?1

可以看到X(k)也是频域上的有限长序列，长度为N。序列X(k)称为序列x(n)的N点DFT。N称为DFT变换区间长度。 通常表示

WN?e?j2?/N

可将定义式表示为

?X(k)?

?x(n)Wn???kn 0?k?N?1

X(k)的离散傅里叶逆变换(IDFT)为

x(n)?1N?

1．圆周移位

?X(k)Wn????kn 0?n?N?1

（二）、DFT的性质

定义序列x(n)的m单位的圆周移位y(n)为：

y(n)?~x(n?m)RN(n)?x((n?m))NRN(n)18

实验五 离散傅立叶变换DFT

(x((n?m))N即对x(n)以N为周期进行周期延拓的序列~RN(n)表x(n)的m点移位，示对此延拓移位后再取主值序列)

2． 圆周卷积

??X1(k) 0?k?N?1 设 x1(n)??NDFT??X2(k) 0?k?N?1 x2(n)??NDFT??X1(k)X2(k) 0?k?N?1 则 x1(n) x2(n)??NDFT这里 x1(n) x2(n) 表示x1(n)与 x2(n)的N点循环卷积。

N?1x1(n) x2(n)??xm?02(m)[x1((n?m))NRN(n)],n?0,1,?,N?1

3． 共轭对称性

x(n)?xep(n)?xop(n),0?n?N?1

19

实验五 离散傅立叶变换DFT

1?*x(n)?[x(n)?x(N?n)]?ep2??,0?n?N?1 1*?xop(n)?[x(n)?x(N?n)]2?DFT??X(k) x(n)??Nxep(n)????NDFT12[X(k)?X(k)]?Re[X(k)]?Xr(k)

*实际应用中，利用上述对称性质可以减少DFT的运算量，提高运算效率。

三、实验内容与步骤

１. 构造离散傅立叶正、反变换函数的MATLAB程序，其中dft(xn,N)为离散傅立叶正变换，idft(xn,N)为离散傅立叶反变换。 function[Xk]=dft(xn,N) n=[0:1:N-1]; k=n;

WN=exp(-j*2*pi/N); nk=n’*k;

WNnk=WN.^nk; Xk=xn*WNnk;

function[xn]=idft(xn,N) n=[0:1:N-1]; k=n;

WN=exp(-j*2*pi/N); nk=n’*k;

WNnk=WN.^(-nk); xn =(Xk*WNnk)/N;

如果x(n)?sin(n?/8)?sin(n?/4)是一个N=16的有限长序列，利用离散傅立叶变换函数求其16点DFT,并显示其DFT结果。

２. 利用MATLAB程序求有限长序列x(n)=8(0.4)n, 0≤n<20的圆周移位

xm(n)?x[(n?10)]20R20(n)，并显示其图形。

程序：

20

实验五 离散傅立叶变换DFT

N=20; m=10; n=0:1:N-1; x=8* (0.4).^n;

n1=mod((n+m),N); %求余

xm=x(n1+1); % 求余后加1是因为MATLAB向量下标从1开始 subplot(2,1,1) stem(n,x); title(‘原序列’)； xlabel(‘n’); ylabel(‘x(n)’); subplot(2,1,2) stem(n,xm);

title(‘圆周移位序列’)； xlabel(‘n’); ylabel(‘xm(n)’);

3.构造一个圆周卷积函数，专门用来计算两个序列（假设均从0点开始）的循环卷积。

Function[y]=circonv(x1,x2) xn2=[x2(1),fliplr(x2)]; xn2(length(xn2))=[]; C=xn2;R=x2; M=toeplitz(C,R); y=x1* M; 如果

x1(n)??(n)?2?(n?1)?3?(n?2)?4?(n?3),

x2(n)?4?(n)?3?(n?1)?2?(n?2)??(n?3)，利用上面的圆周卷积函数求

其4点圆周卷积,并显示其结果。再利用时域圆周卷积定理，调用前面的离散傅立叶正、反变换函数计算它们的4点圆周卷积，并显示其结果。

21

实验五 离散傅立叶变换DFT

4． 调用离散傅立叶正变换函数运行下面MATLAB程序，验证DFT的共轭对

称性。

x=[1 2 4 2 6 32 6 4 2 zeros(1,247)];

x1=[x(1) x(256:-1:2)];

xep=0.5*(x+x1); xop=0.5*(x-x1); XF=dft(x,256); XEPF=dft(xep,256); XOPF=dft(xop,256); clf k=0:255; subplot(2,2,1);

plot(k/128,real(XF));grid; ylabel(‘振幅’)；

title(‘原序列DFT的实部’)；

subplot(2,2,2);

plot(k/128,imag(XF));grid; ylabel(‘振幅’)；

title(‘原序列DFT的虚部’)； subplot(2,2,3);

plot(k/128,real(XEPF));grid; ylabel(‘振幅’)；

title(‘其共轭对称序列DFT’)；

subplot(2,2,4);

plot(k/128,imag(XOPF));grid; ylabel(‘振幅’)；

title(‘其共轭反对称序列DFT’)；

四、实验仪器设备

计算机，MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序；

22

实验五 离散傅立叶变换DFT

六、思考题

１.讨论实验程序1中的.^代表什么含义？

２.讨论实验程序2中圆周移位与线性移位的关系。指令mod的含义是什么？ 中极点与系统稳定性的关系？根据程序运行结果判断该系统的稳定性。 3 .计算实验程序3的理论计算结果，对比两种程序运行结果，结果是否一致？ 4 .讨论实验程序4说明实序列的共轭对称部分的离散傅立叶变换与原序列的傅立叶变换实部的关系，共轭反对称部分的离散傅立叶变换与原序列的傅立叶变换虚部的关系。

23

实验六 快速傅立叶变换FFT及其应用

实验六 快速傅立叶变换FFT及其应用

一、实验目的

１. 利用MATLAB的快速傅立叶变换来计算信号的离散傅立叶变换。 ２. 利用MATLAB程序，理解进一步离散傅立叶变换的物理意义。 3. 利用MATLAB程序，理解快速卷积算法。

二、实验原理

在MATLAB中，使用函数fft可以很容易地计算有限长序列x(n)的离散傅立叶变换X[k]。此函数有两种形式，fft(x)计算序列x(n) 的离散傅立叶变换X(k)，这里X(k)的长度与x(n)的长度相等。fft(x，L)计算序列x(n) 的L点离散傅立叶变换，其中L≥N。若L>N，在计算离散傅立叶变换之前，对x(n)尾部的L-N个值进行补零。同样，离散傅立叶变换序列X(k)的离散傅立叶逆变换x(n)用函数ifft计算，它也有两种形式。 （一）、基本序列的离散傅立叶变换计算

N点离散傅立叶变换的一种物理解释就是，X[k]是x(n)以N为周期的周期延拓序列的

~~离散傅立叶级数系数X(k)的主值区间序列，即X(k)?X(k)RN(k)。例如序列

cos(n)RN(n),当N=16时，cos(n)RN(n)正好是cos(n)的一个周期，所以

888cos(????8n)RN(n)的周期延拓序列就是这种单一频率的正弦序列。而当N=8时，

cos(n)RN(n)正好是cos(n)的半个周期，cos(n)RN(n)的周期延拓就不再是单一频888???率的正弦序列，而是含有丰富的谐波成分，其离散傅立叶级数的系数与N=16时的差别很大，因此对信号进行谱分析时，一定要截取整个周期，否则得到错误的频谱。

（二）、验证N点DFT的物理意义

假如x(n)非周期、有限长，则傅立叶变换存在，那么对X(e?k=2πk/N, k=0,1,…,,N-1取样,则可得X(k)。

24

j?)在N个等间隔频率

实验六 快速傅立叶变换FFT及其应用

X(k)?X(?)??2?k/N???x(n)en????j2?kn/N 0?k?N?1

序列x(n)的N点DFT的物理意义是对X(ω)在[0，2π]上进行N点的等间隔采样。

（三）、利用FFT计算序列的线性卷积

直接计算线性卷积计算量大，并且计算机无法判断y(n)的长度，需要计算多少的y(n)值，若输入为无限长，就更无法计算，其运算量随长度成级数增长。由于可以利用FFT对DFT进行有效的计算，我们希望能够利用DFT来计算线性卷积。

设 x(n) 和 h(n) 是长度分别为M和N的有限长序列， 令 L=M+N-1，定义两个长度L的有限长序列： x'(n)??

h'(n)???h(n),?0,0?n?N?1N?n?L?1?x(n),?0,0?n?M?1M?n?L?1 (3.4.8)

(3.4.9)

通过对x(n) 和 h(n)补充零样本值得到上面两个序列。那么：

yl(n)?x(n)?h(n)?yc(n)?x'(n) h'(n) (3.4.10) 上面的过程如下图所示：

计算线性卷积也可以直接调用函数con来计算，因为MATLAB中的计时比较粗糙，所以只有M和N较大的时候，才能比较两种方法的执行时间快慢。

25

实验六 快速傅立叶变换FFT及其应用

三、实验内容与步骤

１. 对复正弦序列x(n)?ej?8nTLAB程序求当N=16和N=8时的离散傅RN(n)，利用MA

立叶变换，并显示其图形。 程序：

N=16；N1=8； n=0:N-1;k=0: N1-1; x=exp(j*pi*n/8); X1=fft(x,N); X2=fft(x,N1); subplot(2,1,1); stem(n,abs(X1)); axis([0,20,0,20]); ylabel(‘∣X1(k)∣’); title(‘16点的DFT’)；

subplt(2,1, 2);

stem(n, abs(X2)); axis([0, 20, 0, 20]); ylabel(‘∣X2(k)∣’); title(‘8点的DFT’)；

２.已知x(n)?R4(n),X(?)?N=16时的DFT。

1?e?j4??j?1?e, 绘制相应的幅频和相频曲线，并计算N=8和

程序：

N1=8;N2=16;

n=0:N1-1;k1=0:N1-1;k2=0:N2-1; w=2*pi*(0:2047)/2048;

26

实验六 快速傅立叶变换FFT及其应用

Xw=(1-exp(-j*4*w))./(1-exp(-j*));%对x(n)的频谱函数采样2048个点可以近似地看

作是连续的频谱

xn=[(n>=0)&(n<4)];%产生x(n) X1k=fft(xn,N1); X2k=fft(xn,N2);

subplot(3,2,1);plot(w/pi,abs(Xw));xlabel(‘w/π’); title(‘x(n)的幅频曲线’)；

subplot(3,2,2);plot(w/pi,angle(Xw));axis([0,2,-pi,pi]); line([0,2],[0,0]);xlabel(‘w/π’); title(‘x(n)的相频曲线’)；

subplot(3,2,3);stem(k1,abs(X1k),’.’);

xlabel(‘k(w=2πk/N1)’);ylabel(‘∣X1(k)∣’);hold on plot(N1/2*w/pi,abs(Xw));%图形上叠加连续频谱的幅度曲线 subplot(3,2,4);stem(k1,angle(X1k)，’.’); axis([0,N1,-pi,pi]);line([0,N1],[0,0]);

xlabel(‘k(w=2πk/N1)’);ylabel(‘Arg[X1(k)]’);hold on

plot(N1/2*w/pi,angle(Xw));%图形上叠加连续频谱的相位曲线 subplot(3,2,5);stem(k2,abs(X2k),’.’);

xlabel(‘k(w=2πk/N2)’);ylabel(‘∣X2(k)∣’);hold on plot(N2/2*w/pi,abs(Xw));%图形上叠加连续频谱的幅度曲线 subplot(3,2,6);stem(k2,angle(X2k)，’.’);

27

实验六 快速傅立叶变换FFT及其应用

axis([0,N2,-pi,pi]);line([0,N2],[0,0]);

xlabel(‘k(w=2πk/N2)’);ylabel(‘Ang[X2(k)]’);hold on

plot(N2/2*w/pi,angle(Xw));%图形上叠加连续频谱的相位曲线

3．分别利用快速卷积法以及conv函数计算下面两个序列的线性卷积。

h=[3 2 1 -2 1 0 -4 0 3];

↑

x=[1 -2 3 -4 3 2 1 ] ↑

程序1：快速卷积 clf;

h=[3 2 1 -2 1 0 -4 0 3];%冲激 x=[1 -2 3 -4 3 2 1 ]; %输入序列 L=pow2(nextpow2(length(x)+length(h)-1)); Xk=fft(x,L); Hk=fft(h,L); Yk=Xk.*Hk; y=ifft(Yk,L);

nh=0:8;nx=0:6;ny=0:L-1;

subplot(3,1,1);stem(nx,x);title(‘x(n)’);

subplot(3,1,2);stem(nh,h);title(‘h(n)’);

subplot(3,1,3);stem(nx,x); xlabel(‘时间序号n’);ylabel(‘振幅’)；title(‘卷积y(n)’); 程序2： clf;

h=[3 2 1 -2 1 0 -4 0 3];%冲激 x=[1 -2 3 -4 3 2 1 ]; %输入序列 y=conv(h,x); n=0:14; stem(n,y);

28

实验六 快速傅立叶变换FFT及其应用

xlabel(‘时间序号n’);ylabel(‘振幅’)； title(‘用卷积函数conv得到的输出’)；grid;

四、实验仪器设备

计算机，MATLAB软件

五、实验注意事项

课前预先阅读并理解实验程序；

六、思考题

１.分析实验程序1的图形中，两种N值下DFT是否有差别及产生差别的原因。 ２.根据实验程序2的结果图分析由N点DFT的物理意义。

3 .实验程序3两个程序的输出是否一致？计算该卷积的理论结果，与两个程序的输出是否一致？若改变L的长度（改为8时），两个程序的输出是否一致？FFT的变换长度L必须满足什么条件，输出y(n)才等于x(n)和h(n)的线性卷积，这种快速卷积的方法是利用什么运算与线性卷积运算的关系得到的？

29

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

一、实验目的

１. 进一步熟悉IIR数字滤波器的理论知识。

２. 熟悉与IIR数字滤波器设计有关的MATLAB函数。

3 . 学会通过MATLAB，利用脉冲响应不变法和双线性变换法设计IIR数字滤波器，加深对数字滤波器的常用指标和设计过程的理解。

二、实验原理

（一）、低通滤波器的常用指标：

1??P?G(eG(ej?j?)?1??P,for???P

)??S,for?S????

通带边缘频率：?p，阻带边缘频率：?s ， 通带起伏：?p，通带峰值起伏： ?最小阻带衰减：?S??20log（二）、IIR数字滤波器设计

目前，设计IIR数字滤波器的通用方法是先设计相应的低通滤波器，然后再通过双线性变换法和频率变换得到所需要的数字滤波器。模拟滤波器从功能上分有低通、高通、带通及带阻四种，从类型上分有巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器以及贝塞尔滤波器等。

1、利用模拟滤波器设计IIR数字低通滤波器的步骤。

(1)确定数字低通滤波器的技术指标：通带截止频率ωp、通带衰减αp、阻带截止频率ωs、阻带衰减αs。

(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。

脉冲响应不变法：

双线性变换法：

??1tan(?)T2210p??20log10(1??p)[dB]，阻带起伏：?s

(?s)[dB]。

???T(3)按照模拟低通滤波器的技术指标设计模拟低通滤波器。

30

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

(4)将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。

2、下面给出与IIR数字滤波器设计有关的MATLAB文件。

（1）buttord.m

用来确定数字低通或模拟低通滤波器的阶次，其调用格式分别是 a. [N,Wn]=buttord(Wp,Ws,Rp,Rs) b. [N,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)

格式a对应数字滤波器，式中Wp,Ws分别是通带和阻带的截止频率，实际上它们是归一化频率，其值在0-1之间，1对应π（即对π的归一化）。Rp,Rs分别是通带和阻带衰减，单位为dB。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率。

格式b对应模拟滤波器，式中各个变量的含义和格式a相同，但Wp,Ws及Wn是模拟角频率，单位为rad/s。

（2）buttap.m

用来设计模拟低通原型（归一化）滤波器Ha(p)，其调用的格式为 [z , p, k]=buttap(N)

N是欲设计的低通原型（归一化）滤波器的阶次，z, p和k分别是设计出Ha(p)的极点、零点及增益。

（3）lp2lp.m

将模拟低通原型（归一化）滤波器Ha(p)转换为实际的低通滤波器Ha(s)。（去归一化），其调用格式为：

[B，A]=lp2lp(b,a,Wn)

b,a分别是模拟低通原型滤波器Ha(p)的分子、分母多项式的系数向量，其中B，A是去归一化后Ha(s) 的分子、分母多项式的系数向量, Wn为截止频率。

（4）bilinear.m

实现双线性变换，即由模拟滤波器Ha(s)得到数字滤波器H(z)。其调用格式是： [Bz,Az]=bilinear(B,A,Fs)

B，A是去归一化后Ha(s) 的分子、分母多项式的系数向量，Bz，Az是H (z) 的分子、分母多项式的系数向量, Fs是抽样频率。

（4）impinvar.m

由脉冲响应不变法将模拟滤波器Ha(s)转换为数字滤波器H(z)。其调用格式

31

实验七 基于MATLAB的IIR数字滤波器设计(设计实验)

是： [Bz,Az]= impinvar(B,A,Fs)

B，A是去归一化后Ha(s) 的分子、分母多项式的系数向量，Bz，Az是H (z) 的分子、分母多项式的系数向量, Fs是抽样频率。

(5) butter.m

用来直接设计巴特沃斯数字滤波器（双线性变换法），实际上它把buttord.m，buttap.m，lp2lp.m及bilinear.m等文件都包含进去，从而使设计过程更简捷，其调用格式为： a. [B，A]=butter(N, Wn)

b. [B，A]=butter(N，Wn，‘s’)

格式a是设计低通数字滤波器，格式b是设计低通模拟滤波器。B，A是H (z) 的分子、分母多项式的系数向量，Wn是截止频率。

三、实验内容与步骤

１. 设计MATLAB程序，采用脉冲响应不变法设计一个巴特沃斯低通数字滤波器，其通带上限临界频率为400Hz，阻带临界频率为600Hz，抽样频率是1000Hz，在通带内的最大衰减为0.3dB, 阻带内的最小衰减为60dB，并绘出幅频特性曲线。 ２. 设计MATLAB程序，采用双线性变换法设计一个巴特沃斯低通数字滤波器，要求在通带[0，0.2π]内衰减不大于3dB, 在阻带[0.6π,π]内衰减不小于40dB，并绘出幅频特性曲线。

四、实验仪器设备

计算机，MATLAB软件

五、实验要求

根据要求独立编程设计，并根据程序运行结果写出滤波器的系统函数。

32

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

一、实验目的

１. 进一步熟悉FIR数字滤波器的理论知识。

２. 熟悉与FIR数字滤波器设计有关的MATLAB函数。 3. 学会通过MATLAB，利用窗函数法设计FIR数字滤波器。

二、实验原理

设计FIR滤波器实际上是要在满足线性相位的条件下，实现幅度响应的逼近。而一个FIR滤波器若是符合线性相位，则必须满足一定的条件，即：

一个FIR滤波器若是线性相位的，则其单位冲激响应必然满足 h(n)??h(N?1?n) n=0,1,…,N-1 h(n)是关于(N-1)/2对称（奇对称或偶对称） 即,

(1) h(n)是偶对称序列

N?1? ? ? ?

?2 ?,0?n?N?1?h?n??h?N?1?n? (1) h(n) 是奇对称（反对称）序列

N?1??? ? 2 ???h?n???h?N?1?n???

设滤波器要求的理想频率响应为Hd(ejw) , 那么FIR滤波器的设计问题在于

N?1——寻找一系统函数H(z)??h(n)zn?0?n ，使其频率响应H(ejw)?H(z)|z?ejw逼近

Hd(ejw)。若要求FIR滤波器具有线性相位特性，则h(n)必须满足上节所述的对称条件。逼近的方法有三种：

窗口设计法（时域逼近）；频率采样法Frequency-sampling（频域逼近）；

33

最优化设计Optimum Equiripple（等波纹逼近）。

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

窗函数法又称傅立叶级数法，是设计FIR数字滤波器的最简单的方法。FIR数字滤波器的设计问题就是要使所设计的FIR数字滤波器的频率响应H(w)去逼近所要求的理想滤波器的响应Hd(w)。从单位采样响应序列来，就是使所设计滤波器的h(n)逼近理想单位采样响应序列hd(n)，这可以用hd(n)和一个窗函数w(n)的乘积来得到。 （一）、设计原理。

窗函数设计FIR数字滤波器的步骤如下： （1）给定要求的频率响应函数Hd(w)； （2）计算hd(n)；

（3）根据过渡带宽及阻带最小衰减的要求，选定窗的性状以及窗的大小N； （4）根据所选择的合适的窗函数w(n)来修正hd(n)，得到所设计的FIR数字滤波器的单位采样响应序列h (n)= hd(n) w(n)，n=0，1，…，N-1 。 （二）、函数的应用

MATLAB中用fir1函数来设计具有标准频率响应的FIR滤波器。其调用方式： b=fir1(n,wn)——设计n阶低通FIR滤波器，返回的向量b为滤波器的系数（即h(n)的值），它的阶数为n+1；截止频率为wn（对π归一化后的值）。 b=fir1(n,wn，’hign’)——设计n阶高通FIR滤波器 b=fir1(n,wn，’low’)——设计n阶低通FIR滤波器 b=fir1(n,wn，’bandpass’)——设计n阶带通FIR滤波器 b=fir1(n,wn，’stop’)——设计n阶带阻FIR滤波器

b=fir1(n,wn，win)——输入参数win用来指定使用的窗函数的类型，其长度为n+1，缺省情况下，默认为汉明窗。

三、实验内容与步骤

１. 用矩形窗、三角窗、汉宁窗、汉明窗分别设计低通数字滤波器。信号采样频率为1000Hz，数字滤波器的截止频率为 100Hz，滤波器的阶数为80。

程序：

passrad=2*pi*100/1000；

w1=boxcar(81); %矩形窗 w2=triang(81); %三角窗 w3=hanning(81); %汉宁窗

34

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

w4=hamming(81); %矩形窗 n=1:1:81;

hd=sin(passrad*(n-41))./(pi*(n-41)); hd(41)= passrad/pi; h1=hd. *rot90(w1); h2=hd. *rot90(w2); h3=hd. *rot90(w3); h4=hd. *rot90(w4); [MAG1,RAD]=freqz(h1); [MAG2,RAD]=freqz(h2); [MAG3,RAD]=freqz(h3); [MAG4,RAD]=freqz(h4); subplot(2,2,1);

plot(RAD,20*log10(abs(MAG1))); grid on; subplot(2,2,2);

plot(RAD,20*log10(abs(MAG2))); grid on; subplot(2,2,3);

plot(RAD,20*log10(abs(MAG3))); grid on; subplot(2,2,4);

plot(RAD,20*log10(abs(MAG4))); grid on;

运行上面程序，并显示其图形。

２.编写MATLAB程序，利用窗函数法设计线性相位FIR低通数字滤波器，实现对模拟信号采样后进行数字低通滤波，对模拟信号的滤波要求如下：

通带截止频率：fp=2kHz;阻带截止频率：fs=3kHz; 阻带最小衰减：αs=40dB;

采样频率：Fs=10kHz

选择合适的窗函数及其长度，求出h(n)，并画出幅频特性衰减曲线。

35

实验八 基于MATLAB的FIR数字滤波器设计(设计实验)

3．编写MATLAB程序，利用窗函数法设计一线性相位FIR数字低通滤波器，截止频率wc=0.2π，过渡带宽 △w<0.4π, 阻带衰减αs>40dB。

四、实验仪器设备

计算机，MATLAB软件

五、实验注意事项

根据要求独立编程设计。

六、思考题

１. 根据实验程序1的幅频特性曲线，分析不同窗函数下的幅频特性曲线有什么差别。

２. 根据实验内容2和3所编写程序的运行结果，讨论它们的幅频曲线是否符合设计要求？它们的h(n)有什么特点？

36

本文来源：https://www.bwwdw.com/article/72qg.html