医学信号处理 实验指导书

医学信号处理
实验1 离散时间信号的表示
一、实验目的
1．加深对离散时间信号特点的理解
2．掌握MATLAB 中一些常用离散时间信号的表示方法 3．学习MATLAB 中M 文件的使用和调用函数方法
二、实验方法与示例
1．正弦信号
离散正弦序列的MATLAB 表示与连续信号类似，只不过是用stem 函数而不是用plot 函数来画出序列的波形。

下面是正弦序列k ⎪⎭
⎫
⎝⎛6sin π的MATLAB 源程序。

%正弦序列实现程序 k=0:39;
fk=sin(pi/6*k); stem(k,fk)
title('正弦序列sin(pi/6*k)的波形') 程序运行结果如图1–1所示。

5
10
15
20
25
30
35
40
-1-0.8-0.6-0.4-0.200.20.40.60.81正弦序列sin(pi/6*k)的波形
图1–1 正弦序列波形
2．实指数序列
离散实指数序列的一般形式为k
ca ，可用MA TLAB 中的数组幂运算（即点幂运算）c*k a .^来实现。

可直接实现 也可调用用MATLAB 编写绘制离散时间实指数序列波形的函数。

下面为实指数序列的波形函数。

function dszsu(c,a,k1,k2) %c ：指数序列的幅度 %a ：指数序列的底数
%k1：绘制序列的起始序号 %k2：绘制序列的终止序号 k=k1:k2; x=c*(a.^k); stem(k,x,'filled')
3．复指数序列
对于复指数序列，其一般形式为[]k
j k
f k r e
ω=，可以通过调用下面绘制复指数序列时
域波形的MATLAB 函数。

function dfzsu(n1,n2,r,w)
%n1:绘制波形的虚指数序列的起始时间序号 %n2：绘制波形的虚指数序列的终止时间序号 %w ：虚指数序列的角频率 %r: 指数序列的底数 k=n1:n2;
f=(r*exp(i*w)).^k; Xr=real(f); Xi=imag(f); Xa=abs(f); Xn=angle(f);
subplot(2,2,1), stem(k,Xr,'filled'),title('实部'); subplot(2,2,3), stem(k,Xi,'filled'),title('虚部'); subplot(2,2,2), stem(k,Xa,'filled'),title('模'); subplot(2,2,4), stem(k,Xn,'filled'),title('相角');
4．单位抽样序列
可以通过借助MA TLAB 中的零矩阵函数zeros 表示。

全零矩阵zeros(1,N)产生一个由N 个零组成的列向量，对于有限区间的[]k δ可以通过以下MATLAB 程序表示。

% 单位抽样序列实现程序 k=-30:30;
delta=[zeros(1,30),1,zeros(1,30)]; stem(k,delta)
程序运行结果如图1–2所示。

-30
-20
-10
10
20
30
00.10.20.30.40.50.60.70.80.91单位抽样序列的波形
-30
-20
-10
10
20
30
00.10.20.30.40.50.60.70.80.91单位阶跃序列的波形
图1–2 单位抽样序列波形 图1–3 单位阶跃序列波形
5．单位阶跃序列
可以通过借助MATLAB 中的单位矩阵函数ones 表示。

单位矩阵ones(1,N)产生一个由N 个
1组成的列向量，对于有限区间的[]k u
可以通过以下MA TLAB程序表示。

% 单位阶跃序列实现程序
k=-30:30;
uk=[zeros(1,30),ones(1,31)];
stem(k,uk)
程序运行结果如图1–3所示。

6．Matlab函数、及其调用方法
在MA TLAB语言中，M文件有两种形式：脚本和函数。

脚本没有输入/输出参数，只是一些函数和命令的组合。

它可以在MA TLAB环境下直接执行，也可以访问存在于整个工作空间内的数据。

由脚本建立的变量在脚本执行完后仍将保留在工作空间中可以继续对其进行操作，直到使用clear命令对其清除为止。

函数是MA TLAB语言的重要组成部分。

MATLAB提供的各种工具箱中的M文件几乎都是以函数的形式给出的。

函数接收输入参数，返回输出参数，且只能访问该函数本身工作空间中的变量，从命令窗或其他函数中不能对其工作空间的变量进行访问。

读者可以在信号处理的基本理论和MATLAB信号处理工具箱函数的基础上，自己编写一些子程序以便调用。

(1) 函数结构
MATLAB语言中提供的函数通常由以下五个部分组成：
函数定义行：以function开头，函数名(必须与文件名相同)及函数输入输出参数在此定义；
H1行：第一注释行，供lookfor和help在线帮助使用；
函数帮助文件；通常包括函数输入输出参数的含义，调用格式说明；
函数体：它包括进行运算和赋值的所有MATLAB程序代码。

函数体中可以包括流程控制、输入/输出、计算、赋值、注释以及函数调用和脚本文件调用等。

在函数体中完成对输出参数的计算；
注释。

这五个部分中最重要的是函数定义行和函数体。

函数定义行是一个MATLAB函数所必需的，其他各部分的内容可以没有，这种函数称为空函数。

例：function dszsu(c,a,k1,k2)
%c：指数序列的幅度
%a：指数序列的底数
%k1：绘制序列的起始序号
%k2：绘制序列的终止序号
k=k1:k2; %定义序列的起止范围
x=c*(a.^k); %实指数序列的运算
stem(k,x,'filled') %画图
此函数为自定义的实现实数序列的函数，此函数包括有函数定义行，函数帮助文件，函数体，注释四部分。

2．函数调用
函数调用的过程实际上就是参数传递的过程。

例如，在一个脚本文件里调用函数“max”可采用如下面的方式进行。

n=1:20;
a=sin(2*pi*n/20);
[Y,I]=max(a);
该调用过程把变量“a ”传给了函数中的输入参数“x ”,然后把函数运算的返回值传给输出参数“Y ”和“I ”。

其中，Y 是a 序列的最大值，I 是最大值Y 对应的坐标值。

在数字信号处理的基本理论和MATLAB 信号处理工具箱函数的基础上，可以自己编写一些子程序以便调用。

（1）单位抽样序列0()n n δ-的生成函数impseq.m function [x,n]=impseq(n0,ns,nf)
n=[ns:nf];x=[(n-n0)==0]; %序列起点为ns ，终点为nf ，在n=n0处生成一个单位脉冲。

（2）单位阶跃序列0()u n n -的生成函数stepseq.m
function [x,n]=stepseq(n0,ns,nf)
n=[ns:nf];x=[(n-n0)>=0]; %序列起点为ns ，终点为nf ，在n=n0处生成单位阶跃。

三、实验内容
1．画出正弦序列2sin()7
k π，计算其周期性，并根据画出的序列图进行验证。

2．在一幅图柄中画出四个实指数波形n
ca ，其中2c =，a 值分别为
4
3
,45,43,45--，并分析a 取值不同时实指数序列的收敛性。

3．画出虚指数波形k j k j e e
24
,π，分析这两个序列的周期性。

4．画出复指数序列波形[]k
j k
e
r k f ω=，其中（1） 1.7,4r ωπ==。

（2）0.7,6r ωπ==。

（3）1,3r ωπ==。

分析当r>1时，当0<r<1时，当r=1时，复指数序列的实部和虚部的变化趋势。

5．编写程序来产生下列基本脉冲序列。

（1）单位抽样序列，起点n s =0，终点n f =10，在n 0=3处有一单位脉冲。

（2）单位阶跃序列，起点n s =0，终点n f =10，在n 0=3前为0，在n 0=3后为1。

四、实验要求
1．每个实验内容都要求建立一个脚本文件来实现。

2．在实验内容实现时，学习函数文本的建立和调用。

3．在实验报告中写出完整的自编程序，并给出实验结果和对应的分析。

4．在绘制图形时要求图形的完整性，每个图形都要有图标标示。

五、思考题
1．无限长信号和周期信号在计算机中是如何表示的？ 2．离散时间信号和连续时间信号区别是什么？
3．如何实现连续时间信号的离散化，离散化后是否会有误差存在？
实验2 序列的基本运算
一、实验目的
1．加深对序列基本运算的理解
2．掌握MATLAB 中常用序列基本运算的实现方法
3．继续熟悉MATLAB 中M 文件的使用和调用函数方法
二、实验方法与示例
1．加法
对于离散序列来说，序列相加是将信号对应时间序号的值逐项相加，在这里不能像连续时间信号那样用符号运算来实现，而必须用向量表示的方法，即在MATLAB 中离散序列的相加需表示成两个向量的相加，因而参加运算的两序列向量必须具有相同的维数。

下面为实现离散序列相加的MATLAB 实用子程序。

function [f,k]=lsxj(f1,f2,k1,k2)
%实现f(k)=f1(k)+f2(k),f1,f2,k1,k2是参加运算的二离散序列及其对应的时间序列向量，f 和k 为返回%的和序列及其对应的时间序列向量
k=min(min(k1),min(k2)):max(max(k1),max(k2));%构造和序列长度 s1=zeros(1,length(k));s2=s1; %初始化新向量
s1(find((k>=min(k1))&(k<=max(k1))==1))=f1;%将f1中在和序列范围内但又无定义的点赋值为零
s2(find((k>=min(k2))&(k<=max(k2))==1))=f2;%将f2中在和序列范围内但又无定义的点赋值为零
f=s1+s2; %两长度相等序列求和 stem(k,f,'filled')
axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)]) 2．乘法
与离散序列加法相似，这里参加运算的两序列向量必须具有相同的维数。

下面为实现离散时间信号相乘的MATLAB 实用子程序。

function [f,k]=lsxc(f1,f2,k1,k2)
%实现f(k)=f1(k)+f2(k),f1,f2,k1,k2是参加运算的二离散序列及其对应的时间序列向%量,f 和k 为返回的和序列及其对应的时间序列向量
k=min(min(k1),min(k2)):max(max(k1),max(k2));%构造和序列长度 s1=zeros(1,length(k));s2=s1; %初始化新向量
s1(find((k>=min(k1))&(k<=max(k1))==1))=f1;%将f1中在和序列范围内但又无定义的点赋值为零
s2(find((k>=min(k2))&(k<=max(k2))==1))=f2;%将f2中在和序列范围内但又无定义的%点赋值为零
f=s1.*s2; %两长度相等序列求积 stem(k,f,'filled')
axis([(min(min(k1),min(k2))-1),(max(max(k1),max(k2))+1),(min(f)-0.5),(max(f)+0.5)]) 3．移位
序列移位0()()y n x n n =-的生成函数sigshift.m 可以用下面的MATLAB 程序实现。

function [y,n]=sigshift(x,m,n0)
n=m+n0;y=x 4．翻褶
序列翻褶()()y n x n =-的生成函数sigfold.m 可以用下面的MA TLAB 程序实现。

function [y,n]=sigfold(x,n) y=fliplr(x);n=-fliplr(n) 5．卷积和
求卷积和直接采用MATLAB 中的函数conv ，即y=conv(x,h)；它默认序列从n=0开始。

但是如果序列是从一负值开始，即
{}{}
21:)(21:)(nh n nh n h nx n nx n x ≤≤≤≤ 其中nx1<0或nh1<0，或两者同时为负，这样就不能直接采用conv 函数。

其卷积结果序列为{}():1122y n nx nh n nx nh +≤≤+，这样就可构成一个新的卷积函数conv_m ，求出带下标的序列卷积。

function [y,ny]=conv_m(x,nx,h,nh)
ny1=nx(1)+nh(1);ny2=nx(length(x))+nh(length(h)); ny=[ny1:ny2]; y=conv(x,h)
三、实验内容
[]{}1232,3,1,2,3,4,3,1,k -34,()(2);()();()(2)
f k k f k f k f k f k f k f k =≤≤=-=-=-+1.已知序列对应的值范围绘出以下信号的图形：2．已知两离散序列分别[]{}
[]{}122,1,0,1,21,1,1f k f k =--=，求出两序列的和
[][]k f k f 21+,两个序列的乘积[][]k f k f 21⨯，并绘出[][]12,f k f k 及[][]k f k f 21+，[][]k f k f 21⨯的波形。

3．已知两离散序列分别[]{}[]{}1,1,12,1,0,1,221=--=k f k f ，求出两序列的线性卷积
和[][]12f k f k *，并绘出波形。

4．已知序列x()(0.4)(),()()(4)n n u n y n u n u n ==--，求两序列的线性卷积和h (n )，并绘出x (n )，y (n )以及两序列的线性卷积和h (n )的波形。

四、实验要求
1．要求熟练掌握建立脚本文件来实现实验内容。

2．在实验内容是实现时，要熟练掌握函数文本的建立和调用。

3．在实验报告中写出完整的自编程序，并给出实验结果。

4．在绘制图形时要求图形的完整性，每个图形都要有图标标示。

5．在绘制图形时，坐标最好要一致，以方便查看对应点的运算及卷积和长度的变化。

五、思考题
1．长度不等的序列能否实现加法和乘法？ 2．如何实现序列的减法和除法，试编制程序。

3．实现序列翻褶的方法还有那些？
实验3 离散系统的MATLAB 实现
一、实验目的
1．加深对离散系统单位抽样响应和差分方程的理解 2．掌握MATLAB 对离散系统的实现方法
二、实验方法与示例
1．已知系统单位抽样响应的MATLAB 实现
例3–1： 设线性移不变系统的抽样响应()(0.9)()n
h n u n =，输入序列
()()(1
x n u n u n =--，求系统的输出y (n )。

解：系统的输出y (n )为输入x (n )与系统的单位抽样响应h (n )的卷积，即y (n )=x (n )* h (n )，可直接采用conv_m 函数求输出序列。

程序清单如下。

n=-5:50;x=stepseq(0,-5,50)-stepseq(10,-5,50); h=((0.9).^n).*stepseq(0,-5,50);
subplot(3,1,1);stem(n,x);axis([-5,50,0,2]);ylabel('x(n)'); subplot(3,1,2);stem(n,h);axis([-5,50,0,2]);ylabel('h(n)');
[y,ny]=conv_m(x,n,h,n);subplot(3,1,3);stem(ny,y);axis([-5,50,0,8]); xlabel('n');ylabel('y(n)');
程序运行结果如图3–1所示。

-505101520253035404550
012
x (n )
-505101520253035404550
012
h (n )
-5
5
10
15
20
25
30
35
40
45
50
0246
8n y (n )
图3–1 线性移不变系统的输出图
2．线性定常系统的系统函数模型的表示与转换 （1）LSI 系统的系统函数常用模型
1）多项式模型：()n
n n m
m m a s a s a b s b s b s G ++++++=-- 1
10110 用[]
[]n m a a a den b b b num ,,,,,,1010 ==；来表示系统函数 2）零极点模型：()()()()
()()()
n m p s p s p s z s z s z s K
s G ------= 2121
用[][]K k p p p P z z z Z n m ===;,,,;,,,2121 来表示系统函数 3）部分分式模型：()()s k p s r s G n
i i
+-=
∑=11
用[][][]n m n n k k k k p p p P r r r r -===,,,;,,,;,,,102121 来表示系统函数
（2）线性定常系统传递函数模型转换格式
1）系统函数多项式形式转换为零极点形式：
[z,p,k]=tf2zp(num,den)
2）系统函数多项式形式和部分分式之间的相互转换：
[num,den]=residue(r,p,k) [r,p,k]=residue(num,den)
3）系统函数零极点形式转换为多项式形式：
[num,den]=ze2tf(z,p,k)
3．已知系统线性常系数差分方程或系统函数的MATLAB 实现
LSI 系统的差分方程或系统函数可以写成如下形式：
在 MA TLAB 中，filter 函数专用来在给定差分方程或系统函数系数和输入时求差分方程的数值解，函数最基本的使用方法为：
y=filter(b, a, x)或者 [y,zf] = filter(b,a,X,zi)
其中b= [b 0, b 1, …, b M ]； a= [a 0, a 1, …, a N ]；x 为输入序列组，zi 为初始条件。

例3–2：设线性时不变系统的差分方程为
)()(0
k n x b k n y a M
k k
N k k
-=-∑∑==其中
b=1，
a=[1,-1,0.9]，零初始状态。

输入序列x (n )=u (n ),求系统的输出y (n )。

解：在MATLAB 中用差分方程的系数b ，a 来表示一个系统。

程序清单如下。

N=100;n=0:N-1; b=1;
a=[1,-1,0.9];
x=[n>=0];y=filter(b,a,x); subplot 211
stem(n,x, '.');axis([0,N,-1,3]);ylabel('x(n)'); subplot 212
stem(n,y,'.');axis([0,N,-1,3]);ylabel('y(n)'); 程序运行结果如图3–2所示。

∑
∑==-=-M m m N k k m n x b k n y a 0
0)
()(00()M
k
k k N
k k k b z H z a z -=-==∑
∑
0102030405060708090100
-1
012
3x (n )
0102030405060708090100
-1
012
3y (n )
图3–2 线性移不变系统的输出
三、实验内容
1．设线性移不变系统为两个子系统级联构成，两个抽样子系统单位抽样响应分别为
12()(0.8)(),()()(4)n h n u n h n n n δδ==--，输入序列x (n )=u (n )，求系统的输出y (n )。

2． 设线性移不变系统的差分方程为：
()0.7(1)0.45(2)0.6(3)0.8()0.44(1)0.36(2)0.02(3)
y n y n y n y n x n x n x n x n +-----=--+-+-
输入序列分别为x 1(n)=δ(n), x2(n)=R 30(n)时求系统的输出y1(n)和y2(n)。

3．设线性移不变系统的系统函数为：12
12
2.24 2.49 2.24()10.40.75z z H z z z
----++=-+，判定此系统的线性特性。

（叠加定理：两信号加权和的输出是否等于输出的加权和？）
四、实验要求
1．在实验报告中写出完整的自编程序，并给出实验结果和分析。

2．在绘制图形时要求图形的完整性，每个图形都要有图标标示。

五、思考题
1．对于一个离散系统，在MA TLAB 中有几种表示方法，分别是什么？ 2．实验内容1中若系统为两个子系统并联构成，求系统的输出？
实验4 离散系统的Z 域分析
一、实验目的
1．加深对离散系统的分析方法的掌握
2．掌握使用MATLAB 绘制离散系统的零极图的方法
3．掌握使用MATLAB 绘制单位冲激响应和幅频响应的方法 4．掌握使用MATLAB 进行Z 域的部分分式展开的方法 5．掌握使用MATLAB 进行Z 变换和Z 反变换的方法
二、实验方法与示例
1．MATLAB 几个信号处理工具箱函数 （1）residuez ：求极点留数分解
离散系统的输出 Y(z)=X(z)H(z)=B(z)/A(z)，它必然是z 的有理分式：
+++-++-+-=+++++++++=
=------------1111)1(1)1(1)2()1()(1)()2(1)2()1(1)1()1()()2()1()1()()2()1()()()(z k k z
N p N r z p r z p r z N A z N A z A A z N B z N B z B B z A z B z Y N
N M
M
从而得出其时域信号为：
+-+++++=)1()2()()1()()()()()2()2()()1()1()(n k n k n u N p N r n u p r n u p r n y n n n δδ
其中k 是当M≥N 时的直接项，即有限序列，而其余的都是无限序列。

其中的r,p,k 向量可以由下列语句求得
[r,p,k]=residuez(b,a) 式中，b 和b 分别为()z F 的分子多项式和分母多项式的系数向量，r 为部分分式的系数向量，p 为极点向量，k 为多项式的系数向量。

（2）．impz ：求H(z)的Z 反变换h(n)
[h,T]=impz(B,A,N)
h 为存放h(n)的列向量，时间变量N 存放在列向量T 中，当N 为标量时，表示T=[0：N-1]‟，计算h(n)，n=0，1，2，…，N-1；当N 为向量时，T=N ，仅计算N 指定的整数点上的h(n)。

（3）freqz ：求数字滤波器H(z)的频率响应函数
H=freqz[B,A,w] ① [H,w]=freqz[B,A,M] ②
①计算由向量w 指定的数字频率点上数字滤波器H(z)的频率响应，结果存在H 向量中。

②计算出M 个频率点上的频率响应，存放在H 向量中，M 个频率存放在w 向量中，freqz 函数自动将这M 个频率点均匀设置在频率范围[0，π]之间。

若缺省w 和M 时，函数自动选取512个频率点计算。

不带输出向量的freqz 函数将自动绘制幅频和相频曲线。

（4）zplane ：绘制H(z)的零极点图
zplane(b,a)
绘制出向量b(分子多项式系数)中的零点（以符号“о”表示）和向量a （分母多项式系数）
中的极点（以符号“×”表示）以及参考单位圆，并在多阶零点和极点的右上角标出其阶数。

如果b 和a 为矩阵，则会以不同颜色绘出b 和a 各列中的零点和极点。

注意：()H z 是以Z 的降幂次序排列，则系数向量一定要由多项式的最高幂次开始，一直到常数项，缺项要用0补齐。

（5）roots ：求多项式的极点
p=roots(A)
其中A 为待求根的多项式的系数构成的行向量，返回向量p 则包含该多项式所有的根位置列向量。

（6）ztrans 和iztrans ：z 变换及反变换
Z 变换的函数：F=ztrans(f) Z 反变换的函数： f=iztrans(F)
上面两式中，右端的f 和F 分别为时域表示式和Z 域表示式的符号表示，可应用函数sym 来实现，其调用格式为：S=sym(A) ，式中，A 为待分析的表示式的字符串，S 为符号化的数字或变量。

2．利用MATLAB 绘制离散系统的零极图 例4–1： 已知某离散系统的系统函数为
()1
31
4
5+-+=
z z z z H 试用MATLAB 求出该系统的零极点，并画出零极点分布图，判断系统是否稳定。

解：% 绘制零极点分布图的实现程序
a=[3 -1 0 0 0 1]; b=[1 1];
zplane(b,a);xlabel('实部');ylabel('虚部');title('零极点分布图') p=roots(a) %求出极点 q=roots(b) %求出零点 程序运行如果为： p =
0.7255 + 0.4633i 0.7255 - 0.4633i -0.1861 + 0.7541i -0.1861 - 0.7541i -0.7455 q = -1
绘制的系统零极点如图4–1所示。

由图可知，该系统的所有极点均位于Z 平面的单位圆内，故该系统为稳定系统。

-1
-0.5
00.51
-1
-0.8-0.6-0.4-0.200.2
0.40.60.814
实部
虚部
零极点分布图
图4–1 系统零极点图
3．利用MATLAB 绘制离散系统的单位冲激响应和幅频响应 例4–2： 已知某离散系统的系统函数为
()1
31
4
5+-+=
z z z z H 求系统的单位脉冲响应和幅频响应，并判断系统的是否稳定。

解： % 由系统函数求解系统脉冲响应,频率响应实现程序
a=[3 -1 0 0 0 1]; b=[1 1]; h=impz(b,a); subplot 211
stem(h);xlabel('k');ylabel('h(n)');title('系统的单位脉冲响应'); [H,w]=freqz(b,a); subplot 212
plot(w/pi,abs(H));
xlabel('频率\omega');ylabel('H(e jw )');title('系统的频率响应（幅频）')
程序运行结果如图4–2所示。

由系统的单位脉冲响应可判断系统是稳定的。

01020
3040506070
-0.2
00.20.4
0.6k
h (n )
系统的单位脉冲响应
0.1
0.2
0.3
0.4
0.50.6
0.7
0.8
0.9
1
00.511.5
2频率ω
H (e j w )
系统的频率响应（幅频）
图4–2 系统的单位脉冲相应和幅频特性
4．利用MATLAB 实现z 域的部分分式展开式 例4–3： 利用MA TLAB 计算()3
21431818
-----+=
z z z z F 的部分分式展开式。

解：% 部分分式展开式的实现程序
num=[18];
den=[18 3 -4 -1];
[r,p,k]=residuez(num,den) 程序运行结果为： r =
0.3600 0.2400 0.4000 p =
0.5000 -0.3333 -0.3333 k = []
也就是可得()1112
0.360.240.4
10.510.33310.333F z z z z ---=
++-++（）
5．利用MATLAB 实现Z 变换和Z 反变换
例4–4： 求（1）()()()k u ak n f cos =的Z 变换；(2)()()
2
a z az
z F -=的Z 反变换。

解：（1）Z 变换的MATLAB 程序：
f=sym('cos(a*k)'); F=ztrans(f)
程序运行结果为： F =
(z-cos(a))*z/(z^2-2*z*cos(a)+1) （2）Z 反变换的MATLAB 程序 F=sym('a*z/(z-a)^2'); f=iztrans(F)
程序运行结果为： f = a^n*n
三、实验内容
1．已知某离散系统的系统函数为
()3
.0005.05.01
22
32+--++=z z z z z z H 试用MATLAB 求出该系统的零极点，并画出零极点分布图，求系统的单位冲激响应和幅频响应，并判断系统的是否稳定。

2．利用MATLAB 极点()1
322+-=z z z
z H 的部分分式展开式。

3．计算9.0||,
))
9.01()9.01(1
)(121>+-=
--z z z z X 的Z 反变换。

4．求常用信号的Z 变换。

（1）指数序列()n u a n
（2）阶跃序列()n u
四、实验要求
1．在实验报告中写出完整的自编程序，并给出实验结果和分析。

2．在绘制图形时要求图形的完整性，每个图形都要有图标标示。

五、思考题
1．系统的零极点与系统的时域和Z 域之间的对应关系？ 2．系统Z 域分析和复数域分析有何不同？
3．实验中实现Z 变换的函数ztrans 为单边Z 变换，并且未给出收敛域，试编写双边Z 变换函数，并给出收敛域。

（可借助symsum 函数）
实验5 序列的离散傅里叶变换（DFT ）
一、实验目的
1．加深对离散傅里叶变换概念和性质的理解
2．掌握使用MATLAB 对序列进行离散傅里叶变换 3．掌握使用MATLAB 对序列的离散频谱进行分析
二、实验方法和示例
1． M ATLAB 信号处理工具箱提供一些函数
（1）fft 和ifft ：一维快速正傅里叶变换和反傅里叶变换
一维快速正傅里叶正变换：X=fft(x,N) 一维快速正傅里叶反变换：x=fft(X,N)
采用FFT 算法计算序列向量x 的N 点DFT ，缺省N 时，fft 函数自动按x 的长度计算DFT ，若N 小于x 长度时，自动取x 前N 点进行DFT 运算。

当N 为2的整数次幂时，fft 按基2算法计算，否则用混合基算法。

Ifft 的调用格式类似。

（2）fftshift ：移动坐标轴，使频谱对称
Y=fftshift(X)
用来重新排列X=fft(x)的输出，当X 为向量时，把X 的左右两半进行交换，从而将零频分量移至频谱的中心；如果X 为二维傅里叶变换的结果，它同时将X 的左右和上下部分进行交换。

2．利用MATLAB 实现的快速傅里的变换
例5–1：已知有限长序列()n x 的长度4=N ，且()⎪⎪⎩⎪
⎪⎨
⎧==-===3
3
211201n n n n n x ，用FFT 求()k X ，再
用IFFT 求()n x 。

解：利用快速傅里叶变换函数求解的MA TLAB 实现程序清单如下：
xn=[1,2,-1,3]; X=fft(xn) x=ifft(X)
程序运行结果为：
X =
5.0000 2.0000 + 1.0000i -5.0000 2.0000 - 1.0000i x =
1 2 -1 3
例5–2：已知复正弦序列()()n R e n x N n
j 8
1
π
=，余弦序列)()8
cos()(2n R n n x N π
=，分别对序列求当16=N 和8=N 时的DFT ，并绘出幅频特性曲线，并分析两种N 值下DFT 是否有
差别，及产生原因。

解：程序清单如下：
N=16;N1=8;
n=0:N-1;k=0:N1-1;
x1n=exp(j*pi*n/8); %产生x1(n)
X1k=fft(x1n,N); %计算N 点DFT[x1(n)] X2k=fft(x1n,N1); %计算N1点DFT[x1(n)] x2n=cos(pi*n/8); %产生x2(n)
X3k=fft(x2n,N); %计算N 点DFT[x2(n)] X4k=fft(x2n,N1); %计算N1点DFT[x2(n)]
Subplot(2,2,1);stem(n,abs(X1k),'.');axis([0,20,0,20]);ylabel('|X1(k)|') title('16点的DFT[x1(n)]')
Subplot(2,2,2);stem(n,abs(X3k),'.');axis([0,20,0,20]);ylabel('|X2(k)|') title('16点的DFT[x2(n)]')
Subplot(2,2,3);stem(k,abs(X2k),'.');axis([0,20,0,20]);ylabel('|X1(k)|') title('8点的DFT[x1(n)]')
Subplot(2,2,4);stem(k,abs(X4k),'.');axis([0,20,0,20]);ylabel('|X2(k)|') title('8点的DFT[x2(n)]') 程序运行结果如图5–1所示。

N 点离散傅里叶变换的一种物理解释就是()k X 是()n x 以N 为周期的周期延拓序列的
离散傅里叶级数系数)(~
k X 的主值区间序列，即)()(~)(k R k X k X N
=。

当16=N 时， )(1n x 和)(2n x 正好分别是n j
e
8
π
、)8
cos(
n π
的一个周期，所以)(1n x 和)(2n x 的周期延拓序列就是这两个单一频率的正弦序列，其离散傅里叶级数的系数分别如图5–1中的16点的DFT[x1(n)]和16点的DFT[x2(n)]所示。

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

5
10
15
20
51015
20|X 1(k )|
16点的DFT[x1(n)]
05101520
51015
20|
X 2(k )|
16点的DFT[x2(n)]
5
10
15
20
051015
20
|X 1(k )|
8点的DFT[x1(n)]
5
10
15
20
051015
20|X 2(k )|
8点的DFT[x2(n)]
图5–1 基本序列的离散傅里叶变换
3．MATLAB 验证N 点DFT 的物理意义
例5–2： 已知()()n R n x 4=，ωωω
j
j j e e
n x FT e X ----==11)]([)(4，绘制相应的幅频和相频曲
线，并计算图示8=N 和16=N 时的DFT 。

解：程序清单如下：
clear all;close all;
N1=8;N2=16; % 两种FFT 的变换长度 n=0:N1-1;k1=0:N1-1; k2=0:N2-1; w=2*pi*(0:2047)/2048;
Xw=(1-exp(-j*4*w))./( 1-exp(-j*w));
%对x(n)的频谱函数采样2048个点可以近似的看作是连续的频谱 xn=[(n>=0)&(n<4)]; %产生x(n)
X1k=fft(xn,N1); %计算N1=8点的X1(k) X2k=fft(xn,N2); %计算N2=16点的X2(k) subplot(3,2,1);plot(w/pi,abs(Xw));xlabel('w/')
subplot(3,2,2);plot(w/pi,angle(Xw));axis([0,2,-pi,pi]);line([0,2],[0,0]);xlabel('w/π') subplot(3,2,3);stem(k1,abs(X1k),'.');xlabel('k(ω=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(ω=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(ω=2πk/N2)');ylabel('|X2(k)|');hold on plot(N2/2*w/pi,abs(Xw))
subplot(3,2,6);stem(k2,angle(X2k),'.');axis([0,N2,-pi,pi]);line([0,N2],[0,0]); xlabel('k(ω=2πk/N2)') ;ylabel('Arg[X2(k)]');hold on plot(N2/2*w/pi,angle(Xw))
程序运行结果如图5–2所示。

0.5
1 1.5
2
2
4w/π
0.5
1 1.5
2
-2
2w/π
2
468
24k(ω=2πk/N1)
|X 1(k )|
2
468
-202k(ω=2πk/N1)
A r g [X 1(k )]
5
101520
24k(ω=2πk/N2)
|X 2(k )|
51015
-202k(ω=2πk/N2)
A r g [X 2(k )]
图5–2 离散傅里叶变换和傅里叶变换的采样关系
序列()n x 的N 点DFT 的物理意义是对)(ωj e X 在[0，2π]上进行N 点的等间隔采样。

图5–2可以直观的看到()k X 与)(ω
j e X 之间的采样关系。

注意：用DFT （或FFT ）对模拟信号分析频谱时，最好将()k X 的自变量k 换算成对应的模拟频率k f ，作为横坐标绘图，便于观察频谱。

这样，不管变换区间N 取信号周期的几倍，画出的频谱图中有效离散谐波谱线所在的频率值不变。

变换公式如下式：
11
, 0,1,2,,1s k p
F f k k k k N N NT T =
===- 三、实验内容
1．设()6()x n R n =，用MA TLAB 程序计算N 取不同长度(8,32,64)时()n x 的离散傅里叶变换，画图并分析三种情况下异同，并分析原因。

2．对以下序列进行谱分析⎪⎩
⎪
⎨⎧≤≤-≤≤-=⎪⎩⎪
⎨⎧≤≤-≤≤+=其它n n n n n n x 其它n n n n n n x ,074,330,4)(,07
4,83
0,1)(21
选择FFT 的变换区间N 为8和16 两种情况进行频谱分析。

分别画出幅频特性曲线。

并进行对比、分析和讨论。

3．对以下周期序列进行谱分析4()cos 4
x n n π
=，5()cos(/4)cos(/8)x n n n ππ=+ ，选择FFT 的变换区间N 为8和16 两种情况分别对以上序列进行频谱分析。

分别画出其幅频特性曲线。

并进行对比、分析和讨论。

4．设)52.0cos()48.0cos()(n n n x ππ+= （1）取90≤≤n ，求()k X 1；
（2）将（1）中的()n x 补零加长到990≤≤n ，求()k X 2； （3）增加采样值的个数，取990≤≤n ，求()k X 3；
（4）要求画出（1）（2）（3）中的时域和频域序列图，分析频域序列之间区别，并分析原因。

四、实验要求
1．在实验报告中写出完整的自编程序，并给出实验结果。

要求自编实验程序不能只采用两
种方法中的一种。

2．实验结果用图像打印显示，结果分析要详细。

3．在绘制图形时要求图形的完整性，每个图形都要有图标标示。

4．在绘制图形时，坐标最好要一致，以方便对比。

五、思考题
1．对于周期序列，如果周期不知道，如何用FFT 进行谱分析？ 2．如何选择FFT 的变换区间？（包括非周期信号和周期信号）
实验6 序列的离散傅里叶变换（DFT ）的应用
一、实验目的
1．加深对序列离散傅里叶变换的理解
2．掌握使用MATLAB 实现卷积的两种方法
3．掌握使用MATLAB 对连续时间信号进行频谱分析
二．实验方法和示例
1．MATLAB 实现快速卷积
例6–1：用快速卷积法计算)()4.0sin()(15n R n n x =，)(9.0)(20n R n h n
=两个序列的卷积。

快速卷积计算原理框图如图6–1所示。

L 点的FFT
L 点的FFT
L 点 IDFT
x (n )
h (n )
y (n )
图6–1 快速卷积计算原理框图
解：程序清单如下：
M=15;N=20;nx=1:15;nh=1:20; xn=sin(0.4*nx);hn=0.9.^nh;
L=pow2(nextpow2(M+N-1));%L 变为2的整数次幂 Xk=fft(xn,L); Hk= fft(hn,L); Yk=Xk.*Hk;
yn=ifft(Yk,L);ny=1:L
subplot(3,1,1);stem(nx,xn,'.');title('x(n)'); subplot(3,1,2);stem(nh,hn,'.');title('h(n)');
subplot(3,1,3);stem(ny,real(yn),'.');title('y(n)'); 程序运行结果如图6–2所示。

05
1015
-101x(n)
02468101214161820
00.51h(n)
010203040506070
-5
05y(n)
图6–2 )(),(n h n x 及其线性卷积波形
若N 和M 分别是)(n x 和)(n h 的长度，FFT 的变换长度L 必须满足1-+≥M N L ，输出)(n y 才等于)(n x 和)(n h 的线性卷积。

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

2．MATLAB 实现用DFT 对连续信号作谱分析
例6–2：设)50cos()100sin()200cos()(t t t t x a πππ++=，用DFT 分析)(t x a 的频谱结构，选择不同的截取长度T P ，观察存在的截断效应。

选取的参数：（1）采样频率400Hz s f =，s f T 1=；（2） 采样信号序列()()()n w nT x n x a =，
()n w 是窗函数，选取矩形窗为窗函数：函数()()n R n w N =；（3）对)(n x 作2048点DFT 作
为)(t x a 的近似连续频谱()jf X a 。

其中N 为采样点数，p s T f N ⋅=，p T 为截取时间长度，
取三种长度：0.04s 、4⨯0.04s 、8⨯0.04s 。

解：程序清单如下：
fs=400;T=1/fs; %采样频率和采样间隔 Tp=0.04;N=Tp*fs; %采样点数N
N1=[N,4*N,8*N]; %设定三种截取长度 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=[0:4095]/4096/T;
subplot(3,1,m);plot(fk,abs(Xk)/max(abs(Xk)));
if m==1 title('时域截取长度为0.04s 的信号的幅频特性');end if m==2 title('时域截取长度为4*0.04s 的信号的幅频特性');end if m==3 title('时域截取长度为8*0.04s 的信号的幅频特性');end end
程序运行结果如图6–3所示。

050100150200250300350400
00.5
1时域截取长度为0.04s 的信号的幅频特性
050100150200250300350400
00.5
1时域截取长度为4*0.04s 的信号的幅频特性
050100150200250300350400
0.5
1时域截取长度为8*0.04s 的信号的幅频特性
图6–3 DFT 对连续信号作谱分析
图6–3从上到下截取的长度依次分别是N N N 8,4,，由于截断使原频谱中的单频谱线展宽（也成为“泄露”），截取的长度越长，泄露越少，频率分辨越高。

当截取长度为N (p T 为0.04s)时，25Hz 和50Hz 两根谱线已经分辨不清楚了。

另外，在本来应该为零的频段上出现了一些参差不齐的小谱包，成为谱间干扰，其大小取决于加窗的类型。

四、实验内容
1．用快速卷积法计算)(2.0)(5n R n x n
=，)()2.0sin()(8n R n n h =两个序列的卷积。

2．对模拟周期信号6()
cos8cos16cos 20x t t t t πππ=++进行谱分析，选择采样频
率Hz f s 64=，变换区间64,32,16=N 三种情况进行谱分析。

分别画出其幅频特性，并进行分析和讨论。

五、实验要求
1．在实验报告中写出完整的自编程序，并给出实验结果。

要求自编实验程序不能只采用两种方法中的一种。

2．实验结果用图像打印显示，结果分析要详细。

3．在绘制图形时要求图形的完整性，每个图形都要有图标标示。

4．在绘制图形时，坐标最好要一致，以方便对比。

五、思考题
1．两种求卷积的方法得到的卷积结果都会一致吗？
2．对连续时间信号进行谱分析时是否有误差存在，如果有那造成误差的原因有哪些？ 3．实验内容2中的采样频率s f 变化时，信号的频谱会有什么变化？。