工學(xué)MATLAB在信號(hào)與系統(tǒng)中的應(yīng)用

第6章在信號(hào)與系統(tǒng)中的應(yīng)用例6.1連續(xù)信號(hào)的MATLAB描述單位沖激函數(shù)單位階躍函數(shù)：復(fù)指數(shù)函數(shù)時(shí)域變換:主要是在時(shí)域范圍內(nèi)對(duì)信號(hào)進(jìn)行信號(hào)的平移、反折、倒相以及信號(hào)的尺度變換。、移位y=subs(f,t,t-t0);ezplot(y)、反折y=subs(f,t,-t);ezplot(y)、尺度變換

y=subs(f,t,a*t) ezplot(y)、倒相y=-fezplot(y)信號(hào)的運(yùn)算與變換syms

tf=sym("(t/2+1)*(heaviside(t+2)-heaviside(t-2))")subplot(2,3,1),ezplot(f,[-3,3])y1=subs(f,t,t+2)subplot(2,3,2),ezplot(y1,[-5,1])y2=subs(f,t,t-2)subplot(2,3,3),ezplot(y2,[-1,5])y3=subs(f,t,-t)subplot(2,3,4),ezplot(y3,[-3,3])y4=subs(f,t,2*t)subplot(2,3,5),ezplot(y4,[-2,2])y5=-fsubplot(2,3,6),ezplot(y5,[-3,3])已知信號(hào);f2(t)=

sin(

2

t

),用MATLAB繪出滿(mǎn)足下列要求的信號(hào)波形。syms

tf1=sym("(-1*t+4)*(heaviside(t)-heaviside(t-4))"subplot(2,3,1),ezplot(f1)f2=sym("sin(2*pi*t)")subplot(2,3,4),ezplot(f2,[-4,4])y1=subs(f1,t,-t)f3=f1+y1subplot(2,3,2),ezplot(f3)f4=-f3subplot(2,3,3),ezplot(f4)f5=f2*f3subplot(2,3,5),ezplot(f5)f6=f1*f2subplot(2,3,6),ezplot(f6)例6.2LTI系統(tǒng)的零輸入響應(yīng)n階線(xiàn)性時(shí)不變連續(xù)系統(tǒng)的微分方程已知y及其各階導(dǎo)數(shù)的初始值為y(0),y(1)(0)，…，y(n-1)(0)，求系統(tǒng)的零輸入響應(yīng)。解：方程的解為p1,

p2…,pn是方程a1

n+a2

n-1+…+

an

+

an+1

=0的根,C1，…，Cn由y及其各階導(dǎo)數(shù)的初始值來(lái)確定。例6.2LTI系統(tǒng)的零輸入響應(yīng)(續(xù))C1+

C2+…+Cn

=

y0

y0

=

y(0)p1C1+

p2C2+…+

pnCn=Dy0例6.2LTI系統(tǒng)的零輸入響應(yīng)(續(xù))即

V·C

=

Y0其解為C

=V

\Y0式中

V為范德蒙矩陣，在MATLAB的特殊矩陣庫(kù)中有

vander。調(diào)用方法:V=vander(p)例6.3

n階LTI系統(tǒng)的沖激響應(yīng)n階微分方程，寫(xiě)成系統(tǒng)函數(shù)為：

沖擊響應(yīng)就是H(s)的拉普拉斯反變換，可以把H(s)展開(kāi)為極點(diǎn)留數(shù)式。其反變換為MATLAB為用戶(hù)提供了專(zhuān)門(mén)用于求連續(xù)系統(tǒng)沖激響應(yīng)及階躍響應(yīng)，并繪制其時(shí)域波形的函數(shù)

impulse和step。在調(diào)用函數(shù)impulse()和step()時(shí)，我們需要用向量來(lái)對(duì)連續(xù)系統(tǒng)進(jìn)行表示。例6.4卷積的計(jì)算根據(jù)卷積公式：因此編程的過(guò)程為：（1）寫(xiě)出h(t)的MATLAB表達(dá)式；（2）寫(xiě)出u(t)的MATLAB表達(dá)式；（3）利用MATLAB的卷積語(yǔ)句y=conv(u，h)求解（4）畫(huà)曲線(xiàn)plot(t,y)。例6.5

LTI系統(tǒng)的零狀態(tài)響應(yīng)設(shè)二階連續(xù)系統(tǒng)求其沖激響應(yīng)。若輸入為u

=3t

+cos

(0.1t)，求其零狀態(tài)響應(yīng)y(t)。解：特征方程

2+2 +

8

=

0

求出其特征根為p1，p2及相應(yīng)的留數(shù)r1，r2，則沖擊響應(yīng)為：輸出y(t)可用輸入u(t)與沖擊響應(yīng)h(t)的卷積求得。例6.6

有重極點(diǎn)時(shí)的計(jì)算

n級(jí)放大器，每級(jí)的傳遞函數(shù)均為

0/(s+

0)，求階躍響應(yīng)，畫(huà)出n不同時(shí)的波形和頻率特性。解：階躍信號(hào)下系統(tǒng)的輸出為求Y(s)的拉普拉斯反變換，即得到階躍響應(yīng)y(t)。遇到的困難是重極點(diǎn)，公式復(fù)雜，且結(jié)果不穩(wěn)定。為了避開(kāi)重極點(diǎn)問(wèn)題，可以有意把極點(diǎn)拉開(kāi)一些，例如設(shè)n個(gè)極點(diǎn)散布在-0.95

0到1.05

0之間，那樣也就可當(dāng)非重極點(diǎn)來(lái)列程序。例6.7方波分解為多次正弦波之和圖示的周期性方波，其傅里葉級(jí)數(shù)為分別計(jì)算直到9次諧波，并做圖。圖6.7-1輸入周期性方波--x例6.8

全波整流信號(hào)的頻譜周期信號(hào)f(t)可展開(kāi)為直流與各次諧波之和，即其中，是基波角頻率，T為周期。

全波整流電壓Us(t)的波形如圖所示。用傅里葉級(jí)數(shù)可求得其偶次諧波幅值k=2,4,6,…奇次諧波為零。例6.8

全波整流信號(hào)功率（續(xù)）信號(hào)每周期的有效值Us1可由數(shù)字積分求得取前n項(xiàng)分量的功率求出的有效值Us2為用ε可以衡量所取諧波是否包含了原波形的主體。例6.9

周期信號(hào)的濾波

圖示濾波電路，如激勵(lì)電壓us(t)為全波整流信號(hào)，求負(fù)載R兩端的直流和各次諧波分量。解：輸入信號(hào)為多頻的網(wǎng)絡(luò)的輸出由分壓確定例6.9

周期信號(hào)的濾波(續(xù))信號(hào)幅度隨頻率而變?cè)拖到y(tǒng)函數(shù)都是頻率的函數(shù)因此輸出電壓為由此式可求得UR的各次諧波。例6.10調(diào)幅信號(hào)通過(guò)帶通濾波器已知帶通濾波器的系統(tǒng)函數(shù)為激勵(lì)電壓u1(t)=(1+cost

)cos

(100t)求（1）帶通濾器的頻率響應(yīng)；（2）輸出的穩(wěn)態(tài)響應(yīng)u2(t)并畫(huà)出波形。解：將激勵(lì)信號(hào)u1(t)展開(kāi)為傅立葉級(jí)數(shù)例6.10調(diào)幅信號(hào)通過(guò)濾波器(續(xù))帶通濾波器的頻率響應(yīng)為

可畫(huà)出其頻率響應(yīng)，并求各個(gè)分量通過(guò)濾波器后的幅度和相位變化，再將各分量疊合，得到按此信號(hào)作出圖形與原有信號(hào)比較。例6.11

方波的頻譜分析

將積分上下限定為0～10s，并將t分成N等份，用求和代替積分。這樣，傅立葉變換式可寫(xiě)為

求和可以用f(t)行向量乘以e-j

tn列向量來(lái)實(shí)現(xiàn)。式中Δt是t的增量，在程序中，用dt表示。例6.11方波的頻譜分析(續(xù))求不同處的F值，都用同一公式，這就可以利用MATLAB中的矩陣運(yùn)算。將

設(shè)為一個(gè)行數(shù)組，代入上式，則可寫(xiě)為（程序中

用w表示）

其中，F(xiàn)是與w等長(zhǎng)的行向量，t‘是列向量，w是行向量，t’*w是一矩陣，其行數(shù)與t相同，列數(shù)與w相同。這樣，此式就完成了傅里葉變換。類(lèi)似地也可得到傅里葉逆變換表示式為例6.11

方波的頻譜分析(續(xù))

算得的時(shí)域信號(hào)波形及其頻譜圖如右。

下圖為采樣周期較低時(shí)的情況，有明顯的頻率泄漏。例6.12信號(hào)通過(guò)濾波器

計(jì)算幅度為1，寬度為5s的矩形脈沖（同例6.11）通過(guò)下列濾波器的響應(yīng)。（1）理想低通濾波器，（2）低通濾波器解：濾波器輸出的頻譜Y(j

)=F(j

)·H(j

)其時(shí)間響應(yīng)y(t)是Y(j)的傅里葉反變換。例6.12信號(hào)通過(guò)濾波器(續(xù))（1）理想低通濾波器的截止角頻率

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

0～10的部分，用MATLAB語(yǔ)言表述，輸出頻率分量對(duì)應(yīng)的的下標(biāo)數(shù)組為n2

=

find

((w>=

wc)

&

(w<=

wc)其對(duì)應(yīng)的頻率數(shù)組為w2=w(n2)

頻段內(nèi)的頻譜數(shù)組為F2=F(n2:)，它就是濾波后的頻譜數(shù)組Y2，其逆變換即y2

=

F2*exp(j*n2"*t)/pi*dw例6.12信號(hào)通過(guò)濾波器(續(xù))（2）三階低通濾波器的頻率響應(yīng)濾波器的輸出為例6.13離散信號(hào)的MATLAB表述編寫(xiě)MATLAB程序，產(chǎn)生下列基本脈沖序列：?jiǎn)挝幻}沖序列，起點(diǎn)n0，終點(diǎn)nf，在ns處有一單位脈沖(n0≤ns≤nf)。n1=n0:nf;

x1=[zeros(1,ns-n0),1,zeros(1,nf-ns)];用邏輯關(guān)系編程：n1=n0:nf;x1=[(n1-ns)==0]單位階躍序列，起點(diǎn)n0，終點(diǎn)nf，在ns前為0,在ns后為1(n0≤ns≤nf)。n2=n0:nf;x2=[zeros(1,ns-n0),ones(1,nf-ns+1)];復(fù)數(shù)指數(shù)序列。n3

=

n0:nf;

x3=exp((-0.2+0.5j)*n3);例6.14解差分方程的遞推程序描述線(xiàn)性時(shí)不變離散系統(tǒng)的差分方程為編寫(xiě)解上述方程的通用程序。解：

建?？捎眠f推法解差分方程，移項(xiàng)如下：于是得y(n)=(b*us"-a(2::na)*ys")/a(1)這里需要n=0之前的y和u，而MATLAB變量下標(biāo)不能取負(fù)數(shù)。需要作一些技巧性的處理。例6.14

遞推解差分方程(續(xù))

本例中的處理方法是另設(shè)兩個(gè)變量ym和um，使ym(k)=ys(k-na+1)，um(k)=us(k-na+1)，這相當(dāng)于把y和u右移na-1個(gè)序號(hào)，故ym和um的第1到na-1位相當(dāng)于y和u在起點(diǎn)之前的初值。

注意在程序中，隨著計(jì)算點(diǎn)的右移，要隨時(shí)更新相應(yīng)于公式中的向量us和ys。例6.15離散系統(tǒng)對(duì)輸入的響應(yīng)描述LTI系統(tǒng)的差分方程為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，求零輸入的響應(yīng)，計(jì)算20步。求單位脈沖的響應(yīng)h(n)，計(jì)算20步。求單位階躍的響應(yīng)g(n)，計(jì)算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二階數(shù)字濾波器的頻響。低通濾波器的系統(tǒng)函數(shù)（傳遞函數(shù)）為求其頻率響應(yīng)。解：利用多項(xiàng)式求值的函數(shù)polyval，分別求出分子分母多項(xiàng)式在z=exp(jw)時(shí)的值，求其比值。它是對(duì)應(yīng)于w數(shù)組的復(fù)數(shù)數(shù)組，其幅值為幅特性，相角為相特性。然后繪圖。也可用信號(hào)處理工具箱中的freqz函數(shù)快速求解，但為了弄清原理，這里不提倡。6.4.1

模型的典型表達(dá)式連續(xù)系統(tǒng)狀態(tài)空間型

設(shè)x為狀態(tài)變量，u為輸入，y為輸出，系統(tǒng)的狀態(tài)方程為：

如果系統(tǒng)是n階的，輸入有nu個(gè)，輸出有ny個(gè)，則A為n×n階，B為n×nu階，C為ny×n階，而D

為ny×nu階矩陣，對(duì)單輸入單輸出(SISO)系統(tǒng)，

ny=nu=1。6.4.1

模型的典型表達(dá)式(續(xù))傳遞函數(shù)型單輸入單輸出(SISO)n階系統(tǒng)的傳遞函數(shù)為知道分子系數(shù)矢量f

=［f

(1)，f

(2)，…，f

(m

+1)］，分母系數(shù)矢量g

=［g(1)，g(2)，…，g

(n+1)］，就惟一地確定了系統(tǒng)的模型（注意系統(tǒng)的階次n）。而對(duì)物理可實(shí)現(xiàn)的系統(tǒng)，必有n≥m。6.4.1

模型的典型表達(dá)式(續(xù))零極增益型對(duì)傳遞函數(shù)分子和分母進(jìn)行因式分解，可得令z

=[z(1)，z(2)，…，z(m)]為系統(tǒng)的零點(diǎn)矢量，p=[p(1)，p(2)，…，p(n)]為系統(tǒng)的極點(diǎn)矢量，k為系統(tǒng)增益，它是一個(gè)標(biāo)量?？梢钥闯觯琀(s)有m個(gè)零點(diǎn)、n個(gè)極點(diǎn)。物理可實(shí)現(xiàn)系統(tǒng)的n≥m，系統(tǒng)的模型將由矢量z，p及增益k惟一確定，故稱(chēng)為零極增益模型。零極增益模型通常用于描述SISO系統(tǒng)。并可以推廣到MISO系統(tǒng)。6.4.1

模型的典型表達(dá)式(續(xù))極點(diǎn)留數(shù)型如果零極增益模型中的極點(diǎn)都是單極點(diǎn)，將它分解為部分分式，可得其中p

=［p(1)，p(2)，…，p(n)］仍為極點(diǎn)矢量，而r

=［r(1)，r(2)，…，r(n)］為對(duì)應(yīng)于各極點(diǎn)的留數(shù)矢量，p，r兩個(gè)矢量及常數(shù)h惟一地決定了系統(tǒng)的模型。6.4.1

典型表達(dá)式的比較

比較一下這四種情況下模型系數(shù)的總個(gè)數(shù)。假定都是SISO系統(tǒng)，階數(shù)為n，則狀態(tài)空間型有

n2+2n+1個(gè)系數(shù)；傳遞函數(shù)型為m+n+1個(gè)(不含

g(1))(注意，由于m≤n，系數(shù)的數(shù)目小于等于

2n+1)；零極增益型的系數(shù)個(gè)數(shù)為n+m+1；而極

點(diǎn)留數(shù)型為2n+1。

因此，傳遞函數(shù)法的待定系數(shù)最少，而狀態(tài)空間法的待定系數(shù)最多。這說(shuō)明了狀態(tài)空間法中有許多冗余的系數(shù)。事實(shí)上，同一個(gè)系統(tǒng)可以有無(wú)數(shù)個(gè)狀態(tài)空間矩陣A，B，C，D的組合來(lái)描述，其

他描述方法則都是惟一的。6.4.1

離散模型的表達(dá)式2．離散系統(tǒng)以上四種表示模型的方法可以全部推廣至離散系統(tǒng)。只是將系數(shù)矩陣后面加小寫(xiě)字母d，便有：狀態(tài)空間型x[n

+

1]

=

Adx[n]

+

Bdu[n]y[n]

=

Cdx[n]

+

Ddu[n]傳遞函數(shù)型零極增益型6.4.1

離散模型的表達(dá)式(續(xù))極點(diǎn)留數(shù)型數(shù)字信號(hào)處理模型二階環(huán)節(jié)型表6.1列出了連續(xù)和離散線(xiàn)性系統(tǒng)的模型式。6.4.2

模型轉(zhuǎn)換MATLAB中各種模型轉(zhuǎn)換的函數(shù)6.4.2

模型轉(zhuǎn)換（續(xù)）傳遞函數(shù)型到零極增益型已知f，g，求z，p，k，即知道多項(xiàng)式求根。可用MATLAB內(nèi)部函數(shù)roots，z=roots(f)，p=roots(g)，k=f(1)/g(1) 零極增益型到傳遞函數(shù)型已知z，p，k，求f，g，即已知根求多項(xiàng)式。可用MATLAB內(nèi)部函數(shù)poly，它是roots的逆運(yùn)算，即有f=poly(z)*k，g=poly(p)6.4.2

模型轉(zhuǎn)換（續(xù)）傳遞函數(shù)型到極點(diǎn)留數(shù)型及反向變換知道傳遞函數(shù)的系數(shù)g求其極點(diǎn)p，方法同上，而求其中某極點(diǎn)處留數(shù)r，可用專(zhuān)用函數(shù)residue，格式為[r，p，h]=residue(f，g)可直接由f，g求出r，p，k。由極點(diǎn)留數(shù)型到傳遞函數(shù)型仍可用同一函數(shù)。[f，g]

=

residue(r，p，h)residue函數(shù)根據(jù)輸入變?cè)臄?shù)目為二或三個(gè)，決定變換的方向。例6.17由傳遞函數(shù)模型轉(zhuǎn)換為零極增益和狀態(tài)空間模型已知描述系統(tǒng)的微分方程為（1）（2）求出它的傳遞函數(shù)模型、零極增益模型、極點(diǎn)留數(shù)模型和狀態(tài)空間模型。解：（1）f=[2，-5，3];g=[2，3，5，9];（2）f

=

[1，3，2];

g

=

[1，5，7，3];本例中用MATLAB的基本函數(shù)編程，在熟練后均可

用工具箱函數(shù)tf2zp，tf2ss及residue函數(shù)求得。例6.18由狀態(tài)空間轉(zhuǎn)為傳遞函數(shù)設(shè)SISO系統(tǒng)的狀態(tài)方程為：

如果系統(tǒng)是n階的，則A為n×n階，B為n×1階，

C為1×n階，而D為1×1階。給定A，B，C，D，就建立了系統(tǒng)模型。對(duì)狀態(tài)方程取拉氏變換，解出H(s)=f(s)/g(s)，得故有f(s)=Cadj(s

I-A)B

+Dg(s)。

g

=poly(eig(A))=det(sI-A)例6.19系統(tǒng)的串聯(lián)、并聯(lián)和反饋系統(tǒng)的串聯(lián)由圖所示，YB=WBUB

=WBWAUA

=WU故W(s)

=

WA(s)

WB(s)多項(xiàng)式相乘由卷積函數(shù)conv實(shí)現(xiàn)，其表示式為：f

=

conv(fA，fB)，

g

=

conv(gA，gB)系統(tǒng)的并聯(lián)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系統(tǒng)的串并聯(lián)和反饋系統(tǒng)的反饋系統(tǒng)的連接方法如圖6.18-3。復(fù)合系統(tǒng)的傳遞函數(shù)故MATLAB表達(dá)式為

f=conv(fA，gB)g=polyadd(conv(fA，fB)，conv(gA，gB)例6.20復(fù)雜系統(tǒng)的信號(hào)流圖計(jì)算

設(shè)信號(hào)流圖中有ki個(gè)輸入節(jié)點(diǎn)，k個(gè)中間和輸出節(jié)點(diǎn)，它們分別代表輸入信號(hào)ui（i=1，2，…，ki）和系統(tǒng)狀態(tài)xj（j=1，2，…，k）。信號(hào)流圖代表它們之間的聯(lián)結(jié)關(guān)系。用拉普拉斯算子表示后，任意狀態(tài)xj可以表為ui和xj的線(xiàn)性組合用矩陣表示，可寫(xiě)成：從而得到傳遞函數(shù)H

=

(I

-