8. MATLAB在信号与系统中的应用全文-移动互联网-

第六章MATLAB在信号与系统中的应用例6.1连续信号的MATLAB描述单位冲激函数单位阶跃函数：复指数函数例6.4卷积的计算根据卷积公式：因此编程的过程为：（1）写出h(t)的MATLAB表达式；（2）写出u(t)的MATLAB表达式；（3）利用MATLAB的卷积语句y=conv(u，h)求解（4）画曲线plot(t,y)。例6.7方波分解为多次正弦波之和图示的周期性方波，其傅里叶级数为分别计算直到9次谐波，并做图。图6.7-1输入周期性方波例6.13离散信号的MATLAB表述编写MATLAB程序，产生下列基本脉冲序列：（1）单位脉冲序列：起点n0，终点nf，在ns处有一单位脉冲(n0≤ns≤nf)。n1=n0:nf;x1=[zeros(1,ns-n0),1,zeros(1,nf-ns)];用逻辑关系编程：n1

=

n0:nf;

x1=[(n1-ns)==0]例6.13离散信号的MATLAB表述编写MATLAB程序，产生下列基本脉冲序列：（2）单位阶跃序列：起点n0，终点nf，在ns前为0,在ns后为1(n0≤ns≤nf)。n2=n0:nf;x2=[zeros(1,ns-n0),ones(1,nf-ns+1)];用逻辑关系编程：n1

=

n0:nf;

x1=[(n1-ns)>=0]例6.13离散信号的MATLAB表述编写MATLAB程序，产生下列基本脉冲序列：（3）复数指数序列：n3

=

n0:nf;x3=exp((-0.2+0.5j)*n3);例6.11

方波的频谱分析频谱分析用福利也变实现，傅里叶变换公式：

积分可以用离散求和来近似。将积分区间t分成N等份，用求和代替积分。这样，傅立叶变换式可近似为求和可以用f(t)行向量乘以e-jtn列向量来实现。式中Δt是t的增量，在程序中，用dt表示。频域离散化将频域取值区间也等分，在各点上的取值为：利用元素群运算简化表达式：令正变化简写为：同理反变换为：例6.11方波的频谱分析(续)求不同处的F值，都用同一公式，这就可以利用MATLAB中的元素群运算。将

设为一个行数组，代入上式，则可写为（程序中

用w表示）

其中，F是与w等长的行向量，t‘是列向量，w是行向量，t’*w是一矩阵，其行数与t相同，列数与w相同。这样，此式就完成了傅里叶变换。类似地也可得到傅里叶逆变换表示式为例6.11

方波的频谱分析(续)算得的时域信号波形及其频谱图如右。下图为采样周期较低时的情况，有明显的频率混叠。例6.12信号通过滤波器

计算幅度为1，宽度为5s的矩形脉冲（同例6.11）通过下列滤波器的响应。（1）理想低通滤波器，（2）低通滤波器解：滤波器输出的频谱Y(j

)=F(j

)·H(j

)其时间响应y(t)是Y(j)的傅里叶反变换。例6.12信号通过滤波器(续)（1）理想低通滤波器的截止角频率

c=10，故只取F(j)中 =

0～10的部分，用MATLAB语言表述，输出频率分量对应的的下标数组为n2

=

find

((w>=

wc)

&

(w<=

wc)其对应的频率数组为w2=w(n2)频段内的频谱数组为F2=F(n2:)，它就是滤波后的频谱数组Y2，其逆变换即y2

=

F2*exp(j*n2"*t)/pi*dw例6.12信号通过滤波器(续)（2）三阶低通滤波器的频率响应滤波器的输出为例6.14解差分方程的递推程序描述线性时不变离散系统的差分方程为编写解上述方程的通用程序。解：

建模可用递推法解差分方程，移项如下：于是得y(n)=(b*us"-a(2::na)*ys")/a(1)这里需要n=0之前的y和u，而MATLAB变量下标不能取负数。需要作一些技巧性的处理。例6.14

递推解差分方程(续)本例中的处理方法是另设两个变量ym和um，使ym(k)=ys(k-na+1)，um(k)=us(k-na+1)，这相当于把y和u右移na-1个序号，故ym和um的第1到na-1位相当于y和u在起点之前的初值。注意在程序中，随着计算点的右移，要随时更新相应于公式中的向量us和ys。例6.15离散系统对输入的响应描述LTI系统的差分方程为y(n)

-

y(n-1)

+

0.9y(n-2)

-

0.5y(n-3)

=

5u(n)

-

2u(n-1)

+

2u(n-2（1）如已知y(0)=-2，y(-1)=2，y(-2)=-0.5，求零输入的响应，计算20步。求单位脉冲的响应h(n)，计算20步。求单位阶跃的响应g(n)，计算20步。解：利用例6.14的通用程序（1）令us=zeros(1，20);ym=[-1/2，2，-2]；（2）令us=[1，zeros(1，19)]；ym=[0，0，0]（3）令um=ones

(1,20)；ym=[0，0，0]；例6.16二阶数字滤波器的频响。低通滤波器的系统函数（传递函数）为求其频率响应。解：利用多项式求值的函数polyval，分别求出分子分母多项式在z=exp(jw)时的值，求其比值。它是对应于w数组的复数数组，其幅值为幅特性，相角为相特性。然后绘图。也可用信号处理工具箱中的freqz函数快速求解，但为了弄清原理，这里不提倡。6.4.1

模型的典型表达式连续系统状态空间型

设x为状态变量，u为输入，y为输出，系统的状态方程为：

如果系统是n阶的，输入有nu个，输出有ny个，则A为n×n阶，B为n×nu阶，C为ny×n阶，而D为ny×nu阶矩阵，对单输入单输出(SISO)系统，

ny=nu=1。6.4.2

模型转换MATLAB中各种模型转换的函数6.4.2

模型转换（续）传递函数型到极点留数型及反向变换知道传递函数的系数g求其极点p，方法同上，而求其中某极点处留数r，可用专用函数residue，格式为[r，p，h]=residue(f，g)可直接由f，g求出r，p，k。由极点留数型到传递函数型仍可用同一函数。[f，g]

=

residue(r，p，h)residue函数根据输入变元的数目为二或三个，决定变换的方向。例6.17由传递函数模型转换为零极增益和状态空间模型已知描述系统的微分方程为（1）（2）求出它的传递函数模型、零极增益模型、极点留数模型和状态空间模型。解：（1）f=[2，-5，3];g=[2，3，5，9];（2）f

=

[1，3，2];

g

=

[1，5，7，3];本例中用MATLAB的基本函数编程，在熟练后均可

用工具箱函数tf2zp，tf2ss及residue函数求得。例6.18由状态空间转为传递函数设SISO系统的状态方程为：

如果系统是n阶的，则A为n×n阶，B为n×1阶，

C为1×n阶，而D为1×1阶。给定A，B，C，D，就建立了系统模型。对状态方程取拉氏变换，解出H(s)=f(s)/g(s)，得故有f(s)=Cadj(s

I-A)B

+Dg(s)。

g

=poly(eig(A))=det(sI-A)例6.19系统的串联、并联和反馈系统的串联由图所示，YB=WBUB

=WBWAUA

=WU故W(s)

=

WA(s)

WB(s)多项式相乘由卷积函数conv实现，其表示式为：f

=

conv(fA，fB)，

g

=

conv(gA，gB)系统的并联Y=

WAU

+

WBU

=

(WA+WB)

U

=

WU故

W(s)

=

WA(s)

+

WB(s)f

=

polyadd(conv(fA，gB)，conv(fB，gA))g

=

conv(gA，gB)例6.19系统的串并联和反馈系统的反馈系统的连接方法如图6.18-3。复合系统的传递函数故MATLAB表达式为

f=conv(fA，gB)g=polyadd(conv(fA，fB)，conv(gA，gB)例6.20复杂系统的信号流图计算

设信号流图中有ki个输入节点，k个中间和输出节点，它们分别代表输入信号ui（i=1，2，…，ki）和系统状态xj（j=1，2，…，k）。信号流图代表它们之间的联结关系。用拉普拉斯算子表示后，任意状态xj可以表为ui和xj的线性组合用矩阵表示，可写成：从而得到传递函数H

=

(I

-Q)-1P这个简明的公式就等价于梅森公式。例6.20由图列出方程为x1

=

ux2

=

x1-x3-x5x3

=

G1x2x4

=

x1+x3-x5x5

=

G2x4x6

=

x3+x5-x7x7

=

G3x6x8

=

Kx7用矩阵表示例6.20写成X=QX+PU

其中系统传递函数为：H

=(I

-Q)-1P例6.21连续系统状态方程的解线性连续系统在输入信号为零时的状态方程表示式为，其中x为n×1的列向量，A为n×n阶的方阵，求其在初始条件作用下的解。解：建模：按线性方程理论，此齐次方程的解为其中此处要把方阵At作为一个整体求指数函数，就需要调用MATLAB矩阵指数函数expm。注意它与exp函数不同。例6.21连续系统状态方程的解如果A是n×n阶，则expm(At)也是n×n阶，这时如果要算一系列的t值所对应的expm(At)，就不可能像标量指数那样用元素群运算方法了，必须用for循环。不仅如此，如果t的长度为nt，则状态转移矩阵

F(t)=expm(At)将是一个n×n×nt的三维矩阵，需要用MATLAB中高维矩阵的概