数字信号处理实验指导书

前 言

数字信号处理是一门理论和工程实践密切结合的课程。为了加深对教学内容的理解，应在学习理论的同时，加强上机实验，深入理解和消化基本理论，锻炼初学者独立解决问题的能力。本课程实验要求学生运用MATLAB编程完成一些数字信号处理的基本功能。 MATLAB是一高效的工程计算语言，它将计算、可视化和编程等功能集于一个易于使用的环境。在MATLAB环境中描述问题计编制求解问题的程序时，用户可以按照符合人们科学思维的方式和数学表达习惯的语言形式来书写程序。MATLAB广泛应用于工业，电子，医疗和建筑等众多领域。其典型应用主要包括以下几个方面： 数学计算； 算法开发； 数据采集； 系统建模和仿真； 数据分析和可视化 科学和工程绘图；

应用软件开发（包括用户界面）。

；

实验1 用MATLAB产生时域离散信号

一、.实验目的：

1、了解常用时域离散信号及其特点

2、掌握用MATLAB产生时域离散信号的方法 二、.实验原理： 1、时域离散信号的概念

在时间轴的离散点上取值的信号，称为离散时间信号。通常，离散时间信号用x(n)表示，其幅度可以在某一范围内连续取值。由于信号处理设备或装置（如计算机、专用的信号处理芯片等）均以有限位的二进制数来表示信号的幅度，因此，信号的幅度也必须离散化。我们把时间和幅度均取离散值的信号称为时域离散信号或数字信号。

在MATLAB语言中，时域离散信号可以通过编写程序直接产生。 2、常用时域离散信号的生成 1) 单位抽样序列 单位抽样序列的表示式为

?(n)??

?1?0n?0?1n?k 或 ?(n?k)?? n?0?0n?0以下三段程序分别用不同的方法来产生单位抽样序列。 例1-1 用MATLAB的关系运算式来产生单位抽样序列。

n1= -5;n2=5;n0=0; n=n1:n2; x=[n==n0]; stem(n,x,'filled');

axis([n1,n2,0,1.1*max(x)]); xlabel('时间(n)');ylabel('幅度x(n)');

title('单位脉冲序列');

运行结果如图1-1所示：

单位脉冲序列10.8幅度x(n)0.60.40.20-5-4-3-2-10时间(n)12345

图1-1

例1-2 用zeros函数和抽样点直接赋值来产生单位抽样序列。

n1=-5;n2=5;k=0; n=n1:n2;

nt=length(n);（取参数的长度） nk=abs(k-n1)+1;{取绝对值}

x=zeros(1,nt);{（行、列）}（zeros生成零矩阵） x(nk)=1;

绘图部分的程序及作图结果与例1-1相同。 例1-3 生成移位的单位脉冲序列。 n1=-5;n2=5;n0=2;

n=n1:n2; x=[(n-n0)==0]; stem(n,x,'filled');

axis([n1,n2,0,1.1*max(x)]); xlabel('时间(n)');ylabel('幅度x(n)'); title('单位脉冲序列');

运行结果如图1-2所示

单位脉冲序列10.8幅度x(n)0.60.40.20-5-4-3-2-10时间(n)12345

图1-2

2) 单位阶跃序列

单位阶跃序列表示式为

n?0n?0?1?1 或 u(n)=?u(n-k)=?n?0n?0?0?0以下三段程序分别用不同的方法来产生单位阶跃序列。 例1-4用MATLAB的关系运算式来产生单位阶跃序列。

n1=-2;n2=8;n0=0; n=n1:n2; x=[n>=n0]; stem(n,x,'filled');

axis([n1,n2,0,1.1*max(x)]); xlabel('时间(n)');ylabel('幅度x(n)'); title('单位阶跃序列'); box

运行结果如图1-3所示：

单位阶跃序列10.8幅度x(n)0.60.40.20-2-10123时间(n)45678

图1-3

例1-5 用zeros和ones函数来产生单位阶跃序列。 n1=-2;n2=8;k=0; n=n1:n2; nt=length(n); nk=abs(k-n1)+1;

x=[zeros(1,nk-1),ones(1,nt-nk+1)];

绘图部分的程序及作图结果与例1-4相同。. 例1-6生成移位的单位阶跃序列。 n1=-10;n2=10;n0=4; n=n1:n2; x=[(n-n0)>=0]; stem(n,x,'filled');

axis([n1,n2,0,1.1*max(x)]); xlabel('时间(n)');ylabel('幅度x(n)'); title('移位的单位阶跃序列'); box

程序运行结果如图1-4所示：

移位的单位阶跃序列10.8幅度x(n)0.60.40.20-10-8-6-4-20时间(n)246810

图1-4

3) 实指数序列

实指数序列的表示式为 x(n)=an 其中a为实数

例1-7 编写产生a=1/2和a=2的实指数连续信号和离散序列的程序 n1=-10;n2=10;a1=0.5;a2=2; na1=n1:0;x1=a1.^na1; na2=0:n2;x2=a2.^na2;

subplot(2,2,1);plot(na1,x1);(图形窗口) title('实指数信号（a<1）'); subplot(2,2,3);stem(na1,x1,'filled'); title('实指数序列（a<1）'); subplot(2,2,2);plot(na2,x2); title('实指数信号（a>1）'); subplot(2,2,4);stem(na2,x2,'filled'); title('实指数序列（a<1）'); box

程序运行结果如图1-5所示：

实指数信号（a<1）15001500实指数信号（a>1）100010005005000-10-5实指数序列（a<1）0005实指数序列（a<1）1015001000150010005005000-10-5000510

图1-5

4）复指数序列

复指数序列的表示式为 x(n)=eσ+j

(

ω)n

当ω=0时，x(n)为实指数序列;当σ=0时，x(n)为虚指数序列，即

ejωn=cos(ωn)+jsin(ωn)

其实部为余弦序列，虚部为正弦序列。

例1-8 编写程序产生σ=-0.1，ω=0.6的复指数连续信号与离散序列。 n1=30;a=-0.1;w=0.6; n=0:n1;

x=exp((a+j*w)*n);

subplot(2,2,1);plot(n,real(x)); title('复指数信号的实部');

subplot(2,2,3);stem(n,real(x),'filled');

title('复指数序列的实部'); subplot(2,2,2);plot(n,imag(x)); title('复指数信号的虚部');

subplot(2,2,4);stem(n,imag(x),'filled'); title('复指数序列的虚部'); box

程序运行结果如图1-6所示

复指数信号的实部10.50-0.5-101复指数信号的虚部0.50102030-0.50102030复指数序列的实部10.50-0.5-10102030010.5复指数序列的虚部-0.50102030

图1-6

5）正（余）弦序列

正（余）弦序列的表示式为 x(n)=Umsin(ω0n+Θ)

例1-9 已知一时域周期性正弦信号的频率为1Hz，振幅值为1V。编写程序在图形窗口上显示两个周期的信号波形，并对该信号的一个周期进行32点采样获得离散信号。

f=1;Um=1;nt=2; N=32;T=1/f;

dt=T/N; n=0:nt*N-1; tn=n*dt;

x=Um*sin(2*f*pi*tn); subplot(2,1,1);plot(tn,x);

axis([0,nt*T,1.1*min(x),1.1*max(x)]); ylabel('x(t)');

subplot(2,1,2);stem(tn,x);

axis([0,nt*T,1.1*min(x),1.1*max(x)]); ylabel('x(n)'); box

程序运行结果如图1-7所示

10.5x(t)0-0.5-100.20.40.60.811.21.41.61.8210.5x(n)0-0.5-100.20.40.60.811.21.41.61.82

图1-7

6）矩形波序列

MATLAB提供有专门函数square用于产生矩形波。其调用格式如下： x=square(t) ：类似于sin(t),产生周期为2π，幅值为±1的方波。 x=square(t，duty)：产生指定周期的矩形波，其中duty用于指定占空比。 将square的参数t换成n，且n取整数，则可以获得矩形序列。

例1-10 一个周期性矩形信号频率为5kHz，信号幅度在0～2V之间，占空比为0.25，。编写程序生成该信号，要求在图形窗口上显示2个周期的信号波形；对信号的一个周期进行16点采样获得离散信号。

f=5000;nt=2; N=16;T=1/f; dt=T/N; n=0:nt*N-1; tn=n*dt;

x=square(2*f*pi*tn,25)+1; subplot(2,1,1);plot(tn,x);

axis([0,nt*T,1.1*min(x),1.1*max(x)]); ylabel('x(t)');

subplot(2,1,2);stem(tn,x);

axis([0,nt*T,1.1*min(x),1.1*max(x)]); ylabel('x(n)'); box

程序运行结果如图1-8所示; 21.5)t(x10.5001234x 10-421.5)(nx10.5001234x 10-4

三、.实验内容：

1、阅读并上机验证实验原理部分的例题程序，理解每一条语句的含义。

改变例题中的有关参数（如信号的频率、周期、幅度、显示时间的取值范围、采样点数等），观察对信号波形的影响。

2、编写程序，产生以下离散序列： （1）f(n)=δ(n) (-3

(0.1+j1.6∏)n

(0

（4）f(n)=3sin(nП/4) (0

3、一个连续的周期性方波信号频率为200Hz，信号幅度在-1～+1V之间，要求在图形窗口上显示其两个周期的波形。以4kHz的频率对连续信号进行采样，编写程序生成连续信号和其采样获得的离散信号波形。

四、.实验预习：

1、预先阅读附录部分的MATLAB基础介绍，认真阅读实验原理，明确本次实验任务，读懂例题程序，了解实验方法。

2、根据实验内容预先编写实验程序。

3、预习思考题：产生单位脉冲序列和单位阶跃序列各有几种方法？如何使用？ 五、实验报告

1、列写调试通过的实验程序，打印实验程序产生的曲线图形。

2、思考题：通过例题程序，你发现采样频率Fs、采样点数N、采样时间间隔dt在程序编写中有怎样的联系？使用时需注意什么问题？

实验2 离散LSI系统的时域分析

一、.实验目的：

1、加深对离散系统的差分方程、单位脉冲响应、单位阶跃响应和卷积分析方法的理解。

2、初步了解用MATLAB语言进行离散时间系统时域分析的基本方法。

3、掌握求解离散时间系统的单位脉冲响应、单位阶跃响应、线性卷积以及差分方程的程序的编写方法，了解常用子函数的调用格式

二、实验原理：

1、离散LSI系统的响应与激励

由离散时间系统的时域分析方法可知，一个离散LSI系统的响应与激励可以用如下框图表示：

x[n]y[n]Discrete-timesystmeNM

其输入、输出关系可用以下差分方程描述：

?ak?0ky[n?k]??bkx[n?m]

k?02、用函数impz和dstep求解离散系统的单位脉冲响应和单位阶跃响应。

例2-1 已知描述某因果系统的差分方程为6y(n)+2y(n-2)=x(n)+3x(n-1)+3x(n-2)+x(n-3) 满足初始条件y(-1)=0，x(-1)=0，求系统的单位脉冲响应和单位阶跃响应。

解： 将y(n)项的系数a0进行归一化，得到

y(n)+1/3y(n-2)=1/6x(n)+1/2x(n-1)+1/2x(n-2)+1/6x(n-3)

分析上式可知，这是一个3阶系统，列出其bk和ak系数： a0=1， a,1=0， a,2=1/3， a,3=0 b0=1/6，b,1=1/2， b,2=1/2， b,3=1/6

程序清单如下： a=[1,0,1/3,0]; b=[1/6,1/2,1/2,1/6]; N=32; n=0:N-1; hn=impz(b,a,n); gn=dstep(b,a,n);

subplot(1,2,1);stem(n,hn,'k'); title('系统的单位序列响应'); ylabel('h(n)');xlabel('n');

axis([0,N,1.1*min(hn),1.1*max(hn)]); subplot(1,2,2);stem(n,gn,'k'); title('系统的单位阶跃响应'); ylabel('g(n)');xlabel('n');

axis([0,N,1.1*min(gn),1.1*max(gn)]); 程序运行结果如图2-1所示：

系统的单位序列响应1.20.51.10.410.90.30.8h(n)系统的单位阶跃响应0.2g(n)0.70.60.10.500.40.30.2010n2030010n2030-0.1

图2-1

3、用函数filtic和filter求解离散系统的单位序列响应和单位阶跃响应。

例2-2 已知描述某因果系统的差分方程为6y(n)-2y(n-4)=x(n)-3x(n-2)+3x(n-4)-x(n-6)，满足初始条件y(-1)=0，x(-1)=0，求系统的单位脉冲响应和单位阶跃响应。时间轴上N取32点。

注意：原式非标准形式，必须化为标准形式后再列出系数b,a。 程序清单如下： x01=0;y01=0; a=[1,0,0,0,-1/3,0,0]; b=[1/6,0,-1/2,0,1/2,0,-1/6]; N=32;n=0:N-1; xi=filtic(b,a,0); x1=[n==0];

hn=filter(b,a,x1,xi); x2=[n>=0]; gn=filter(b,a,x2,xi); subplot(1,2,1);stem(n,hn,'k'); title('系统的单位序列响应'); ylabel('h(n)');xlabel('n');

axis([0,N,1.1*min(hn),1.1*max(hn)]); subplot(1,2,2);stem(n,gn,'k'); title('系统的单位阶跃响应'); ylabel('g(n)');xlabel('n');

axis([0,N,1.1*min(gn),1.1*max(gn)]); 程序运行结果如图2-2所示：

系统的单位序列响应0.60.40.2)n(h0-0.2-0.40102030n系统的单位阶跃响应0.20.10-0.1-0.2-0.30102030ng(n)

图2-2

4、用MATLAB实现线性卷积 1）用函数conv进行卷积运算：

求解两个序列的卷积和，关键在于如何确定卷积结果的时宽区间。MATLAB提供的求卷积函数conv默认两个序列的序号均从n=0开始，卷积结果y对应的序列的序号也从n=0开始。

例2-3 已知两个序列f1=0.8n （0

subplot(2,2,1);stem(n1,f1,'filled'); title('f1(n)'); n2=0:10; N2=length(n2); f2=ones(1,N2);

subplot(2,2,2);stem(n2,f2,'filled'); title('f2(n)'); y=conv(f1,f2);

subplot(2,1,2);stem(y,'filled'); 程序运行结果如图2-3所示：

0

图2-3

2）非零起始序列的卷积运算：

当两个序列不是从0开始时，必须对conv函数稍加扩展。

由卷积原理可知，若待卷积的两个序列序号分别为{x(n)；nx=nxs:nxf}，{h(n)；nh=nhs:nhf}，则卷积和y(n)的序号起点和终点分别为：nys=nxs+nhs，nyf=nxf+nhf。据此可定义通用卷积函数convu：

function[y,ny]=convu(h,nh,x,nx) nys=nh(1)+nx(1);nyf=nh(end)+nx(end); y=conv(h,x);ny=nys:nyf;

例2-4 已知序列f1=0.5n (0

[y,ny]=convu(f1,n1,f2,n2); subplot(2,2,1);stem(n1,f1); subplot(2,2,2);stem(n2,f2); subplot(2,1,2);stem(ny,y); 程序运行结果如图2-4所示：

543210051010.80.60.40.20-505103020100-505101520

图2-4

3）卷积积分的动态过程演示：

为了更深入地理解序列卷积的原理，下面提供一段演示卷积和的动态过程的程序。 例2-5 动态演示例2-3两个序列f1=0.8n （0

lf1=length(n1); nf2=0:10; lf2=length(n2); f2=ones(1,lf2); m=max(lf2,lf1);

if lf2>lf1 nf2=0;nf1=lf2-lf1; elseif lf2

u=[zeros(1,lt),f2,zeros(1,nf2),zeros(1,lt)]; t1=(-lt+1:2*lt);

f1=[zeros(1,2*lt),f1,zeros(1,nf1)]; hf1=fliplr(f1); N=length(hf1); y=zeros(1,3*lt); for k=0:2*lt

p=[zeros(1,k),hf1(1:N-k)]; y1=u.*p; yk=sum(y1); y(k+lt+1)=yk;

subplot(4,1,1);stem(t1,u); subplot(4,1,2);stem(t1,p); subplot(4,1,3);stem(t1,y1); subplot(4,1,4);stem(k,yk); axis([-20,50,0,5]);hold on pause(2); end

程序运行结果如图2-5所示：

10.50-1510.50-1510-1-155-10-5051015202530-10-5051015202530-10-50510152025300-20-1001020304050

图2-5

5、离散LSI系统时域响应的求解：

MATLAB提供了多种方法求解离散LSI系统的响应：

1）用conv函数进行卷积积分，求任意输入的系统零状态响应；

例2-6 已知描述某因果系统的差分方程为6y(n)+2y(n-2)=x(n)+3x(n-1)+3x(n-2)+x(n-3)， 满足初始条件y(-1)=0，x(-1)=0。在该系统的输入端加一个矩形脉冲序列，其占空比为0.25，一个周期取16个采样点，求该系统的响应。 程序清单如下： N=16; n=0:N-1;

x=[ones(1,N/4),zeros(1,3*N/4)]; subplot(2,2,1);stem(n,x); a=[1,0,1/3,0,];

title('极点在单位圆内'); subplot(3,2,2);impz(b1,a1,20); z2=[0.3,0]';p2=[0.6+0.8j,0.6-0.8j]'; [b2,a2]=zp2tf(z2,p2,k); subplot(3,2,3);zplane(z2,p2); title('极点在单位圆上'); subplot(3,2,4);impz(b2,a2,20); z3=[0.3,0]';p3=[1+j,1-j]'; [b3,a3]=zp2tf(z3,p3,k); subplot(3,2,5);zplane(z3,p3); title('极点在单位圆外'); subplot(3,2,6);impz(b3,a3,20);

程序运行结果如图3-3所示。由图可见，这三个系统的极点均为复数且处于z平面的右半平面。由图可知，当极点位于单位圆内，系统的单位序列响应随着频率的增大而收敛；当极点位于单位圆上，系统的单位序列响应为等幅振荡；当极点位于单位圆外，系统的单位序列响应随着频率的增大而发散。由此可知系统1、2为稳定系统。

极点在单位圆内Impulse ResponseImaginary Part0-1-202Real Part极点在单位圆上 Amplitude110-10Imaginary Part1015n (samples)Impulse Response50-1-202Real Part极点在单位圆外 Amplitude120-201015n (samples)Impulse Response5Imaginary Part0-1-20Real Part2 Amplitude15000-500051015n (samples)

图3-3

由以上三例可得结论：系统只有在极点处于单位圆内才是稳定的。 例3-6 已知某离散时间系统的系统函数为

0.2?0.1z?1?0.3z?2?0.1z?3?0.2z?4H(z)?

1?1.1z?1?1.5z?2?0.7z?3?0.3z?4求该系统的零极点及零极点分布图，并判断系统的因果稳定性。 程序清单如下： b=[0.2,0.1,0.3,0.1,0.2]; a=[1,-1.1,1.5,-0.7,0.3]; rz=roots(b) rp=roots(a)

subplot(2,1,1);zplane(b,a); title('系统的零极点分布图');

subplot(2,1,2);impz(b,a,20); title('系统的单位序列响应'); xlabel('n');ylabel('h(n)'); 程序运行结果如下： rz =

-0.5000 + 0.8660i -0.5000 - 0.8660i 0.2500 + 0.9682i 0.2500 - 0.9682i rp =

0.2367 + 0.8915i 0.2367 - 0.8915i 0.3133 + 0.5045i 0.3133 - 0.5045i

系统的零极点分布图1Imaginary Part0.50-0.5-1-3-201Real Part系统的单位序列响应-1230.60.4h(n)0.20-0.202468n1012141618

图3-4

由零极点分布图可见，该系统的所有极点均在单位圆内，因此该系统是一个因果稳定系统。

3、离散系统的频率响应

（1）离散系统的频率响应的基本概念 已知稳定系统传递函数的零极点增益模型为

H(z)?K?(z?cm?1Nn?1Mm)

?(z?d则系统的频响函数为

n)H(ej?)?H(z)z?ej??K?(em?1Nn?1Mj??cm)?Kn?Cm?1Nn?1Mmej?m?H(ej?)ej?(?)

jn?(e??dj)?De?n其中，系统的幅频特性为

H(e)?Kj??Cm?1Nn?1Mm

n?D系统的相频特性为

?(?)???m???n??(N?M)

m?1n?1MN由以上各式可见，系统函数与频率响应有着密切的联系。适当地控制系统函数的零极点分布，可以改变离散系统的频响特性：

① 在原点（z=0）处的零点或极点至单位圆的距离始终保持不变，其值|ej|=1，所以，对幅度响应不起作用；

ω

② 单位圆附近的零点对系统幅度响应的谷值位置及深度有明显影响；

③ 单位圆内且靠近单位圆附近的极点对系统幅度的峰值位置及大小有明显的影响。 （2）系统的频响特性分析

例3-7 已知某离散时间系统的系统函数为

0.1321?0.3963z?2?0.3963z?4?0.1321z?6H(z)?

1?0.34319z?2?0.60439z?4?0.20407z?6求该系统在0～П频率范围内的相对幅频响应与相频响应。 程序清单如下：

b=[0.1321,0,-0.3963,0,0.3963,0,-0.1321]; a=[1,0,0.34319,0,0.60439,0,0.20407]; freqz(b,a);

程序运行结果如图3-5所示。该系统是一个IIR数字带通滤波器。其中幅频特性采用归一化的相对幅度值，以分贝（dB）为单位。

Magnitude (dB)-50-100-1500.10.20.30.40.50.60.70.8Normalized Frequency (?? rad/sample)0.910Phase (degrees)-200-400-600-80000.10.20.30.40.50.60.70.8Normalized Frequency (?? rad/sample)0.91

图3-5

例3-8 已知某离散时间系统的系统函数为

0.2?0.1z?1?0.3z?2?0.1z?3?0.2z?4H(z)?

1?1.1z?1?1.5z?2?0.7z?3?0.3z?4

求该系统在0～П频率范围内的绝对幅频响应与相频响应。 程序清单如下：

b=[0.2,0.1,0.3,0.1,0.2];a=[1,-1.1,1.5,-0.7,0.3]; n=(0:500)*pi/500; [h,w]=freqz(b,a,n);

subplot(2,1,1);plot(n/pi,abs(h));grid; axis([0,1,1.1*min(abs(h)),1.1*max(abs(h))]); xlabel('\\omega/\\pi');ylabel('幅度'); subplot(2,1,2);plot(n/pi,angle(h));grid;

axis([0,1,1.1*min(angle(h)),1.1*max(angle(h))]); xlabel('\\omega/\\pi');ylabel('相位');

程序运行结果如图3-6所示。该系统为一低通滤波器。

10.8幅度0.60.40.200.10.20.30.40.5?/?0.60.70.80.912相位0-200.10.20.30.40.5?/?0.60.70.80.91

图3-6

例3-9 已知某离散时间系统的系统函数为

0.1?0.4z?1?0.4z?2?0.1z?3H(z)? ?1?2?31?0.3z?0.55z?0.2z求该系统在0～П频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及零极点分布图。

程序清单如下： b=[0.1,-0.4,0.4,-0.1]; a=[1,0.3,0.55,0.2]; n=(0:500)*pi/500; [h,w]=freqz(b,a,n); db=20*log10(abs(h));

subplot(2,2,1);plot(w/pi,abs(h));grid; axis([0,1,1.1*min(abs(h)),1.1*max(abs(h))]); title('幅频特性（V）');

xlabel('\\omega/\\pi');ylabel('幅度(V)'); subplot(2,2,2);plot(w/pi,angle(h));grid;

axis([0,1,1.1*min(angle(h)),1.1*max(angle(h))]); xlabel('\\omega/\\pi');ylabel('相位'); title('相频特性');

subplot(2,2,3);plot(w/pi,db);grid axis([0,1,-100,5]); title('幅频特性（dB）'); subplot(2,2,4);zplane(b,a); title('零极点分布');

程序运行结果如图3-7所示：

幅频特性（V）10.82相频特性幅度(V)相位00.5?/?幅频特性（dB）10.60.40.20-200.5?/?零极点分布10Imaginary Part00.511-500-1-101Real Part2-100

图3-7

（3） 一个求解频率响应的实用函数。

在实际使用freqz进行离散系统频响特性分析时。通常需要求解幅频响应、相频响应、群时延，幅频响应又分为绝对幅频和相对幅频两种表示方法。下面定义函数freqz_m，利用该函数，可方便求出上述各项。freqz_m函数定义如下：

function[db,mag,pha,grd，w]=freqz_m(b,a); [H,w]=freqz(b,a,1000,'whole'); H=(H(1:501))';w=(w(1:501))'; mag=abs(H);

db=20*log10((mag+eps)/max(mag)); pha=angle(H); grd=grpdelay(b,a,w);

例3-10 已知某离散时间系统的系统函数为

0.1321?0.3963z?2?0.3963z?4?0.1321z?6H(z)?

1?0.34319z?2?0.60439z?4?0.20407z?6求该系统在0～П频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及群时延。

程序清单如下：

b=[0.1321,0,0.3963,0,0.3963,0,0.1321]; a=[1,0,-0.34319,0,0.60439,0,-0.20407]; [db,mag,pha,grd,w]=freqz_m(b,a); subplot(2,2,1);plot(w/pi,mag);grid axis([0,1,1.1*min(mag),1.1*max(mag)]); title('幅频特性（V）');

xlabel('\\omega/\\pi');ylabel('幅度(V)'); subplot(2,2,2);plot(w/pi,pha);grid; axis([0,1,1.1*min(pha),1.1*max(pha)]); xlabel('\\omega/\\pi');ylabel('相位'); title('相频特性');

subplot(2,2,3);plot(w/pi,db);grid axis([0,1,-100,5]); title('幅频特性（dB）'); subplot(2,2,4);plot(w/pi,grd);grid axis([0,1,0,10]) title('群时延');

程序运行结果如图3-8所示：

幅频特性（V）10.82相频特性幅度(V)相位00.5?/?幅频特性（dB）10.60.40.20-200.5?/?群时延1010-505-10000.51000.51

图3-8

三、实验内容：

1、输入并运行例题程序，理解每条语句的含义。 2、求以下各序列的z变换：

x1(n)?nan x2(n)?sin(?0n) x3(n)?e?ansin(?0n)

3、求下列函数的逆z变换

zzz1?z?3 X1(z)? X2(z)? X3(z)? X4(z)?j?02?1z?a(z?a)z?e1?z4、求一下系统函数所描述的离散系统的零极点分布图，并判断系统的稳定性 （1）

H(z)?z(z?0.3)

(z?1?j)(z?1?j)

4?1.6z?1?1.6z?2?4z?3（2）H(z)?

1?0.4z?1?0.35z?2?0.4z?35、已知某离散时间系统的系统函数为

0.187632?0.241242z?2?0.241242z?4?0.187632z?6H(z)?

1?0.602012z?2?0.495684z?4?0.035924z?6求该系统在0～П频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及群时延。

四、实验预习：

1、认真阅读实验原理部分，明确实验目的，复习有关离散LSI系统频率响应的理论知识。

2、读懂实验原理部分的例题程序，熟悉与本实验有关的MATLAB函数。

3、根据实验内容预先编写实验程序，并思考本实验提出的有关MATLAB函数在调用时应注意哪些问题。

4、预习思考题：

① 系统函数零极点的位置与系统单位序列响应有何关系？ ② 离散系统的零极点对系统幅频响应有何影响？ 五、实验报告：

1、列写调试通过的实验程序，打印实验程序产生的曲线图形。 2、列出本实验提出的有关MATLAB函数在调用时应注意的问题。 3、给出预习思考题答案。

实验4 DFS、DFT与FFT

一、实验目的：

1、加深对周期序列DFS、有限长序列DFT和FFT的基本概念及理论的理解。 2、掌握用MATLAB语言求解DFS、DFT、FFT的以及相应反变换的方法。 3、观察周期序列的重复周期数对序列频谱特性的影响。 二、实验原理：

1、周期序列的离散傅里叶级数（DFS） （1）DFS的基本概念：

?表示。其中，N为离散时间序列x(n)满足x(n)= x(n+rN)，称为离散周期序列，用x(n)信号的周期，x(n)称为离散周期序列的主值。

?可以用离散傅里叶级数（DFS）表示： 周期序列x(n)2πjkn1N-1???x(n)=?X(k)eN=IDFS[X(k)] n=0,1,2,?,N-1

Nk=0?是周期序列DFS第k次谐波分量的系数，也称为周期序列的频谱，可表示为 其中，X(k)??X(k)=?x(n)en=0N-1-j2πknN?=DFS[x(n)] k=0,1,2,?,N-1

以上两式也是周期序列的一对傅里叶级数变换对。

令WN=e-j2πN，以上DFS变换对又可以写成：

nk???X(k)=DFS[x(n)] =?x(n)WN k=0,1,2,?,N-1

n=0N-1N-11-nk???x(n)=IDFS[X(k)]=X(k)W?N n=0,1,2,?,N-1 Nk=0与连续周期信号的傅里叶级数相比，周期序列的离散傅里叶级数有以下特点：

① 连续周期信号的傅里叶级数由无穷多个与基波频率成整数倍的谐波分量叠加而成，而周期为N的周期序列的傅里叶级数仅有N个独立的谐波分量。

?也是一个以N为周期的周期序列。 ② 周期序列的频谱X(k)（2）周期序列的DFS和IDFS

例4-1 已知一个周期性矩形序列的脉冲宽度占整个周期的1/4，一个周期的采样点数

为16点，编程显示3个周期的序列波形，并：

① 用傅里叶级数求信号的幅度谱和相位谱。

② 求傅里叶级数逆变换的图形，并与原序列进行比较。 程序清单如下： N=16;

xn=[ones(1,N/4),zeros(1,3*N/4)]; xn=[xn,xn,xn]; n=0:3*N-1; k=0:3*N-1;

Xk=xn*exp(-j*2*pi/N).^(n'*k); x=(Xk*exp(j*2*pi/N).^(n'*k))/N; subplot(2,2,1);stem(n,xn);

title('x(n)');axis([-1,3*N,1.1*min(xn),1.1*max(xn)]); subplot(2,2,2);stem(n,abs(x));

title('IDFS|X(k)|');axis([-1,3*N,1.1*min(x),1.1*max(x)]); subplot(2,2,3),stem(k,abs(Xk));

title('|X(k)|');axis([-1,3*N,1.1*min(abs(Xk)),1.1*max(abs(Xk))]); subplot(2,2,4),stem(k,angle(Xk));

title('arg|X(k)|');axis([-1,3*N,1.1*min(angle(Xk)),1.1*max(angle(Xk))]); 程序运行结果如图4-1所示：

x(n)10.80.60.40.2001020|X(k)|1210864201020304010-1010230408642010IDFS|X(k)|203040arg|X(k)|203040

图4-1

由离散傅里叶级数逆变换图形可见，与原序列相比，幅度扩大了32倍。这是因为周期序列为原主值序列周期的3倍，做逆变换时未做处理。可将逆变换程序改为：

x=(Xk*exp(j*2*pi/N).^(n'*k))/3*3*N;

由上例可见，周期序列的DFS和IDFS是依据变换公式编程的，无论信号序列如何变化，求解的公式总是一样的。因此，可将其编写成通用子程序：

① 离散傅里叶级数变换通用子程序dfs.m function[Xk]=dfs(xn,N) n=0:N-1; k=0:N-1;

WN=exp(-j*2*pi/N); nk=n'*k; Xk=xn*WN.^nk;

②离散傅里叶级数逆变换通用子程序idfs.m

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

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

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

例4-2 利用上述两个子程序，重做例4-1 程序清单如下： N=16;

xn=[ones(1,N/4),zeros(1,3*N/4)]; Xk=dfs(xn,N); x=idfs(Xk,N);

subplot(2,2,1);stem(n,xn);

title('x(n)');axis([-1,3*N,1.1*min(xn),1.1*max(xn)]); subplot(2,2,2);stem(n,abs(x));

title('IDFS|X(k)|');axis([-1,3*N,1.1*min(x),1.1*max(x)]); subplot(2,2,3),stem(k,abs(Xk));

title('|X(k)|');axis([-1,3*N,1.1*min(abs(Xk)),1.1*max(abs(Xk))]); subplot(2,2,4),stem(k,angle(Xk));

title('arg|X(k)|');axis([-1,3*N,1.1*min(angle(Xk)),1.1*max(angle(Xk))]);

程序运行结果如图4-2所示。由于子程序仅适用于对主值区间进行变换，周期次数无法传递给子程序，因此程序执行结果仅显示一个周期的变换情况。

x(n)10.80.60.40.2001020|X(k)|4321010203040210-1010304010.80.60.40.2010IDFS|X(k)|203040arg|X(k)|203040

图4-2

（2）周期重复次数对序列频谱的影响

理论上讲，周期序列不满足绝对可积条件，因此不能用傅里叶级数来表示。实际处理时可先取K个周期进行处理，然后令K趋于无穷大，分析其极限情况。根据这一分析思路，可以观察序列由非周期到周期变化时，频谱由连续谱逐渐向离散谱过渡的过程。 例4-3 已知一矩形序列的脉冲宽度占整个周期的1/2，一个周期的采样点数为10点，用DFS求序列的重复周期数分别为1、4、7、10时的幅度谱。

程序清单如下： xn=[ones(1,5),zeros(1,5)]; Nx=length(xn); Nw=1000;dw=2*pi/Nw;

k=floor((-Nw/2+0.5):(Nw/2+0.5)); for r=0:3 K=3*r+1;

nx=0:(K*Nx-1); x=xn(mod(nx,Nx)+1); Xk=x*(exp(-j*dw*nx'*k))/K; subplot(4,2,2*r+1);stem(nx,x); axis([0,K*Nx-1,0,1.1]);ylabel('x(n)'); subplot(4,2,2*r+2);plot(k*dw,abs(Xk)); axis([-4,4,0,1.1*max(abs(Xk))]);ylabel('X(k)'); end

程序运行结果如图4-3所示。由图可见，序列的重复周期数越多，频谱约是向几个频点集中。当序列的周期数趋于无穷大时，频谱转化为离散谱。

10.5010.5010.5010.5002040608002040600102030024685X(k)0-45x(n)-2024X(k)0-45x(n)-2024X(k)0-45x(n)-2024X(k)0-4x(n)-2024

图4-3

2、离散傅里叶变换（DFT） （1）DFT与IDFT

在实际中常常使用有限长序列。如果有限长序列为x(n)，则该序列的离散傅里叶变换

对可表示为

nkX(k)=DFT[x(n)] =?x(n)WN k=0,1,2,?,N-1

n=0N-11N-1-nkx(n)=IDFT[X(k)]=?X(k)WN n=0,1,2,?,N-1

Nk=0从离散傅里叶变换定义式可以看出，有限长序列在时域上是离散的，在频域上也是离散的，式中WN=e理带来了方便。

由有限长序列的傅里叶变换和逆变换定义可知，DFT和DFS的变换公式非常相似，因此，在程序编写上也基本一致。

例4-4 已知x(n)=[0,1,2,3,4,5,6,7]，求其DFT和IDFT。要求： ① 画出序列傅里叶变换对应的|X(k)|和arg[X(k)]图形。 ② 画出x(n)图形，并与IDFT[X(k)]图形进行比较。 程序清单如下： xn=[0,1,2,3,4,5,6,7]; N=length(xn); n=0:N-1; k=0:N-1;

Xk=xn*exp(-j*2*pi/N).^(n'*k); x=(Xk*exp(j*2*pi/N).^(n'*k))/N; subplot(2,2,1);stem(n,xn);

title('x(n)');axis([-1,N,1.1*min(xn),1.1*max(xn)]); subplot(2,2,2);stem(n,abs(x));

title('IDFT|X(k)|');axis([-1,N,1.1*min(x),1.1*max(x)]); subplot(2,2,3),stem(k,abs(Xk));

title('|X(k)|');axis([-1,N,1.1*min(abs(Xk)),1.1*max(abs(Xk))]); subplot(2,2,4),stem(k,angle(Xk));

-j2πN，即仅在单位圆上N个等间距的点上取值，这为使用计算机进行处

title('arg|X(k)|');axis([-1,N,1.1*min(angle(Xk)),1.1*max(angle(Xk))]);

程序运行结果如图4-4所示。由图可见，与周期序列不同，有限长序列本身是仅有N点的离散序列，相当于周期序列的主值部分。因此，其频谱也对应序列的主值部分，是长度为N的离散序列。

x(n)6420642IDFT|X(k)|024|X(k)|68024arg|X(k)|68302200-2024680246810

图4-4

（2）DFT与DFS的联系

将周期序列的傅里叶级数变换对和有限长序列的离散傅里叶变换对进行比较可见，两

?换成了有限长序列x(n)，同时，由于式中Wnk的周期性，者的区别仅仅是将周期序列x(n)N因而有限长序列的离散傅里叶变换实际上隐含着周期性。

例4-5 已知周期序列的主值x(n)=[0,1,2,3,4,5,6,7]，求x(n)的周期重复次数为4次时的

DFS。要求：

① 画出主值序列周期序列的波形。

??② 画出周期序列傅里叶变换对应的X(k)和[X(k)]的图形。

程序清单如下： xn=[0,1,2,3,4,5,6,7]; N=length(xn); m=0:N-1; n=0:4*N-1; k=0:4*N-1;

xn1=xn(mod(n,N)+1);

Xk=xn1*exp(-j*2*pi/N).^(n'*k); subplot(2,2,1);stem(m,xn); title('x(n)');

subplot(2,2,2);stem(n,xn1); title('周期序列');

subplot(2,2,3),stem(k,abs(Xk));

title('|X(k)|');axis([-1,4*N,1.1*min(abs(Xk)),1.1*max(abs(Xk))]); subplot(2,2,4),stem(k,angle(Xk));

title('arg|X(k)|');axis([-1,4*N,1.1*min(angle(Xk)),1.1*max(angle(Xk))]);

?程序运行结果如图4-5所示。与例4-4相比，有限长序列x(n)可以看成是周期序列x(n)?可以看成是有限长序列以N为周期的周期延拓。频域的一个周期；反之，周期序列x(n)上的情况也一样。从这个意义上说，周期序列只有有限个序列值有意义。 （3）DFT与DTFT的联系

若离散时间非周期序列为x(n)，则它的离散傅里叶变换（DTFT）对定义为

DTFT[x(n)]=X(e)=?x(n)e-jωn

jωn=-??1πIDTFT[X(e)]=x(n)=?X(ejω)ejωndω

2π-πjωjωjω?(?)jωjω其中X(e)称为序列的频谱。可表示为X(e)=X(e)e，X(e)称为序列的幅

本文来源：https://www.bwwdw.com/article/46tp.html