MATLAB讲义第二章.docx
《MATLAB讲义第二章.docx》由会员分享，可在线阅读，更多相关《MATLAB讲义第二章.docx（19页珍藏版）》请在冰豆网上搜索。
MATLAB讲义第二章
第二章 序列的付氏变换（DTFT）和Z变换
一、序列的DTFT：
X（ejw）
序列x（n）的付氏变换如下：
方法一:
定义法此法只为让大家熟悉DTFT形成的过程
因为MATLAB无法计算连续变量w,只能在
（或考虑对称性，在
）范围内，把w赋值为很密的、长度很长的向量来近似连续变量。
通常最简单的就是赋以K个等间隔的值：
它表示把基本数字频率的范围
均分成K份后，每一份的大小。
则频率向量表示为w=[w1,w2,…,wk]=k
dw,
序列的位置向量为n=[n1,n2,…,nN]
则DTFT的计算式可用一个向量与矩阵相乘的运算来实现：
[X（w1）,X（w2）,…,X（wk）]
=[x（n1）,x（n2）,…,x（nN）]
=[x（n1）,x（n2）,…,x（nN）]
又因为（jnT*w）=jnT*kdw=j*dw*n’*k
所以
Xw=xn*exp（j*dw*n’*k）
例1：
求有限长序列xn=[13531]，n=-1:
3的DTFT,画出它在w=-8~8rad/s范围内的频率特性，讨论其对称性。
再把xn左右移动，讨论时移对DTFT的影响。
解:
根据DTFT定义得
将w在-8~8rad/s之间分为1000份。
>>xn=[13531];nx=-1:
3;
>>w=linspace（-8,8,1000）;%设定频率向量¿
>>Xw=xn*exp（-j*nx'*w）;%用定义计算DTFT
>>subplot（3,5,1）;stem（nx,xn）;axis（[-2,6,0,6]）;title（'Ô Ê¼ÐòÁÐ'）;ylabel（'xn'）;%画序列图
>>subplot（3,5,2）;plot（w,abs（Xw））;title（'·ù¶È'）;%画幅频曲线
>>subplot（3,5,3）;plot（w,angle（Xw））;title（'Ïà½Ç'）%画相频曲线
>>subplot（3,5,4）;plot（w,real（Xw））;title（'ʵ²¿'）%画实频曲线
>>subplot（3,5,5）;plot（w,imag（Xw））;title（'Ð鲿'）%画虚频曲线
>>nx_2=nx+2;%使xn右移2位
>>Xw_2=xn*exp（-j*nx_2'*w）;
>>subplot（3,5,6）;stem（nx_2,xn）;axis（[-2,6,0,6]）;ylabel（'xn_2'）;title（'ÓÒÒÆ2λ'）
>>subplot（3,5,7）;plot（w,abs（Xw_2））;
>>subplot（3,5,8）;plot（w,angle（Xw_2））;
>>subplot（3,5,9）;plot（w,real（Xw_2））;
>>subplot（3,5,10）;plot（w,imag（Xw_2））;
>>nx_1=nx-1;%使xn左移1位
>>Xw_1=xn*exp（-j*nx_1'*w）;
>>subplot（3,5,11）;stem（nx_1,xn）;axis（[-2,6,0,6]）;ylabel（'xn_1'）;title（'×óÒÆ1λ'）
>>subplot（3,5,12）;plot（w,abs（Xw_1））;
>>subplot（3,5,13）;plot（w,angle（Xw_1））;
>>subplot（3,5,14）;plot（w,real（Xw_1））;
>>subplot（3,5,15）;plot（w,imag（Xw_1））;
由上述结果得出以下结论：
（1）序列的DTFT是连续函数。
（2）序列的DTFT是周期函数，周期为2pi。
（3）本题给出的是实序列，实序列的DTFT具有对称性，幅频和实频偶对称;相频和虚频奇对称。
（4）信号的时移不影响幅频特性，只影响相频。
方法二:
调用内部函数法（Freqz）
MATLAB工具箱中有计算DTFT的专用函数freqz。
它的调用方式为：
（1）[h,w]=freqz（b,a,N）
h—复频率响应
N—N点复频率响应，是把
分割的份数，计算的是频率部分的特性。
若省略N，默认512
b,a—滤波器系数的分子分母向量,且以z－1的升幂排列的系数向量b为输入序列xn的系数向量,a为输出序列yn的系数向量,且a0要归一化,即a0=1.
w—数字频率向量，它把0-π均分为N份，w=[0:
N-1]π/N
此式可计算出正频率区间特性。
（2）[h,w]=freqz（b,a,N，’whole’）
用于画出全奈氏频率范围内的特性。
此时N把0-2π均分。
（3）[h]=freqz（b,a,w）
在自己选定的频点向量W上计算频率特性。
（4）freqz（b,a）
若没有给出左端，则直接在当前窗口给出频率响应的幅频响应和相频响应。
例：
差分方程：
y（n）-0.9y（n-1）=0.5x（n）+0.8x（n-1）,求它的频率响应H（ejw），并画图。
并求输入为x（n）=cos（0.1πn）u（n）的稳态输出ys.
>>b=[0.5,0.8];
>>a=[1,-0.9];
>>[h,w]=freqz（b,a）;%求出频率响应（0到π分成500点）
>>subplot（2,1,1）,plot（w,abs（h））,gridon;title（'|h（ejw）|'）
>>subplot（2,1,2）,plot（w,angle（h））,gridon;title（'angle[h（ejw）]'）
>>
稳态输入x（n）的频率为wx=0.1π,初始相角为θx=0，系统在该频点处的响应为：
>>wx=0.1*pi;
>>nwx=floor（wx/pi*512）
nwx=
51
>>hwx=h（nwx）
hwx=
1.2085-4.0139i
>>Ax=abs（hwx）
Ax=
4.1919
>>thetax=angle（hwx）
thetax=
-1.2784
>>y=Ax*cos（0.1*pi*n-thetax）;
>>subplot（2,1,1）,stem（n,x）,title（'x'）;
>>subplot（2,1,2）,stem（n,y）,title（'y'）;
>>subplot（2,1,1）,stem（n,x,'*'）,title（'x'）;
>>subplot（2,1,2）,stem（n,y,'.'）,title（'y'）;
由运行结果Ax=4.1919，thetax=-1.2784得
即
则ys=4.2694cos（0.1πn-1.2784）=4.1919cos[0.1π（n-4.0693）]
三、Z变换的计算
（1）两个序列的Z变换
例1：
x1（n）=[2,3,4],ns1=0;x2（n）=[3456],ns2=0;求两序列的Z变换，及两序列的卷积y（n）及Y（z）.
由Z变换的定义
知，
X1（z）=2z0+3z-1+4z-2
X2（z）=3z0+4z-1+5z-2+6z-3
∵y（n）=x1（n）*x2（n）
∴Y（z）=X1（z）X2（z）
∴Y（z）=（2z0+3z-1+4z-2）（3z0+4z-1+5z-2+6z-3）
＝6＋17z-1+34z-2+43z-3+38z-4+24z-5
则，y（n）=[61734433824]
y（n）=x1（n）*x2（n）=conv（x1,x2）
求卷积程序：
>>x1=[234];nx1=0:
2;x2=[3456];nx2=0:
3;
>>y=conv（x1,x2）
y=
61734433824
即y（n）=[61734433824]
从而Y（z）=＝6＋17z-1+34z-2+43z-3+38z-4+24z-5
由以上看出，两序列的卷积和它们的Z变换的相乘结果是一致的。
（时域卷积定理）
因为Z变换的计算还涉及到收敛域的确定,这是较为复杂的,所以目前没有一个计算Z变换的很好的函数.
三、Z反变换
MATLAB的内部函数residuez可计算离散系统的极点留数。
调用格式：
[r,p,C]=residuez（b,a）.
b,a――H（z）中，分子分母多项式以z－1的升幂排列的系数向量。
p――极点向量,即分母的根
r――部分分式分解后，各极点所对应的留数Ak
C――无穷项多项式系数向量，仅在M≥N时存在。
则
则Z反变换为：
y（n）=r
（1）p
（1）nu（n）+…+r（N）p（N）nu（n）+C
（1）δ（n）+C
（2）δ（n-1）+…
若H（z）是系统函数，则其反变换h（n）就是它的单位脉冲响应，可用函数impz来计算，可用该命令来核对计算。
其调用方式为：
h=impz（b,a,L）
其中L——待求脉冲序列h（n）的长度
实例1：
计算的
反变换
先用函数poly求出分母多项式的系数，再用residuez求X（z）的极点和留数
程序：
>>b=1;a=poly（[0.9,0.9,-0.7]）;
>>[r,p,C]=residuez（b,a）
则其z变换为：
则
r=
0.2461
0.5625
0.1914
p=
0.9000
0.9000
-0.7000
C=
[]
四、用Z变换求系统输出
（1）方法1：
采用residuez函数
实例1：
系统函数
，输入信号x=[2345],用Z变换求出输出y（n）。
解：
∵Y（z）=H（z）·X（z）=
=
∴ B（z）是X（z）与－3z-1的卷积，用函数conv可求得
求出Y（z）后，用函数residuez可求其反变换，即得y（n）。
>>x=[2345];nsx=0;
>>nfx=length（x）-1;
>>b=-3;a=[2-52];
>>B=conv（-3,x）;A=a;
>>[r,p,k]=residuez（B,A）
r=
-10.2500
32.0000
求得
则
p=
2.0000
0.5000
k=
-24.7500-7.5000
>>nf=input（'nf'）;
nf8
求得Yi（z）=
Yd（z）=-24.75-7.5z-1
>>n=0:
nf;
>>yi=（r
（1）*p
（1）.^n+r
（2）*p
（2）.^n）.*stepseq（0,0,nf）;
>>yd=k
（1）*impseq（0,0,nf）+k
（2）*impseq（1,0,nf）;
>>y=yi+yd
>>subplot（1,2,1）,stem（[nsx:
nfx],x）;title（'x'）,line（[0,10],[0,0]）;
>>subplot（1,2,2）,stem（[0:
nf],y）;title（'y'）,line（[0,10],[0,0]）;
（2）方法2：
采用filter函数
若不需求出表达式，可用函数filter来解。
y=filter（b,a,x）
其中：
b,a——差分方程中x（n）,y（n）的系数数组
b=[b0,b1,…,bM]
a=[a0,a1,…,aN]
x——输入序列
当初始条件为0，输入序列为
时，y即为系统的单位脉冲响应h（n）
例:
>>x=[2,3,4,5,zeros（1,10）];b=-3;a=[2-52];
>>y=filter（b,a,x）;
>>ny=length（y）-1
ny=
13
>>subplot（1,2,1）,stem（[0:
ny],x）;title（'x'）,line（[0,10],[0,0]）;
>>subplot（1,2,2）,stem（[0:
ny],y）;title（'y'）,line（[0,10],[0,0]）;
>>
五、Z变换解差分方程
1、递推法（单位脉冲响应）
法1：
采用filter函数
已知a,b向量，给定输入x（n）,则可求得输出y（n），其长度与输入x（n）相同。
当初始条件为0，输入为单位脉冲函数时，系统的输出为脉冲响应h（n）。
>>b=2*[1,2,3];
>>a=[1,-0.8,0.4];
>>Y=[-2,-3];
>>X=[1,1]
X=
11
>>clear
>>b=1;
>>a=[1,-1,0.9];
>>x=[1,zeros（1,7）];
>>y=filter（b,a,x）
y=
1.00001.00000.1000-0.8000-0.8900-0.17000.63100.7840
法2：
采用impz函数
调用方法：
（h,t）=impz（b,a）:
返回参数h是脉冲响应的数据，t是脉冲响应的采样时间间隔
（h,t）=impz（b,a，N）：
N为指定脉冲信号的长度。
（h,t）=impz（b,a,n,Fs）：
Fs指定脉冲信号的采样频率，其默认值为1。
例：
>>x=[2,3,4,5,zeros（1,10）];b=-3;a=[2-52];
>>N=14;
>>[h,t]=impz（b,a,N）;
h=
1.0e+004*
-0.0001
-0.0004
-0.0008
-0.0016
-0.0032
-0.0064
-0.0128
-0.0256
-0.0512
-0.1024
-0.2048
-0.4096
-0.8192
-1.6384
t=
0
1
2
3
4
5
6
7
8
9
10
11
12
13
>>stem（t,h）,title（'h（n）=impz（b,a,14）'）
六、H（z）的零极点
roots函数求其零极点，
zplane（b,a）可绘出零极点图
例：
，求零极点，并作图。
>>b=[2040-2];
>>a=[10.3000];
>>z=roots（b）
极点
z=
-0.0000+1.5538i
-0.0000-1.5538i
0.6436
-0.6436
>>p=roots（a）
零点
p=
0
得到传输函数
0
0
-0.3000
>>zplane（b,a）
>>
升级会员 