MATLAB讲义第二章Word格式.docx
《MATLAB讲义第二章Word格式.docx》由会员分享，可在线阅读，更多相关《MATLAB讲义第二章Word格式.docx（19页珍藏版）》请在冰豆网上搜索。
%用定义计算DTFT
subplot（3,5,1）;
stem（nx,xn）;
axis（[-2,6,0,6]）;
title（'
Ô
Ê
¼
Ð
ò
Á
'
）;
ylabel（'
xn'
%画序列图
subplot（3,5,2）;
plot（w,abs（Xw））;
·
ù
¶
È
%画幅频曲线
subplot（3,5,3）;
plot（w,angle（Xw））;
Ï
à
½
Ç
）%画相频曲线
subplot（3,5,4）;
plot（w,real（Xw））;
Ê
µ
²
¿
）%画实频曲线
subplot（3,5,5）;
plot（w,imag（Xw））;
é
）%画虚频曲线
nx_2=nx+2;
%使xn右移2位
Xw_2=xn*exp（-j*nx_2'
subplot（3,5,6）;
stem（nx_2,xn）;
xn_2'
Ó
Ò
Æ
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'
subplot（3,5,11）;
stem（nx_1,xn）;
xn_1'
×
ó
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;
|h（ejw）|'
subplot（2,1,2）,plot（w,angle（h））,gridon;
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（'
subplot（2,1,2）,stem（n,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变换
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:
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.7000
C=
[]
四、用Z变换求系统输出
（1）方法1：
采用residuez函数
系统函数
，输入信号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）
-10.2500
32.0000
求得
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）;
）,line（[0,10],[0,0]）;
subplot（1,2,2）,stem（[0:
nf],y）;
（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;
y=filter（b,a,x）;
ny=length（y）-1
ny=
13
subplot（1,2,1）,stem（[0:
ny],x）;
ny],y）;
五、Z变换解差分方程
1、递推法（单位脉冲响应）
法1：
已知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
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。
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
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）
零点
得到传输函数
-0.3000
zplane（b,a）
升级会员 