西安電子科技大學(xué)MATLAB基礎(chǔ)及應(yīng)用課件-第9章MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例

第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例連續(xù)信號(hào)和系統(tǒng)離散信號(hào)和系統(tǒng)系統(tǒng)函數(shù)頻譜及其幾何意義場(chǎng)的計(jì)算和偏微分方程的數(shù)值解西安電子科技大學(xué)《MATLAB基礎(chǔ)及應(yīng)用》第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例9.1

連續(xù)信號(hào)和系統(tǒng)【例

9-1-1】

列出單位脈沖、單位階躍、復(fù)指數(shù)函數(shù)等連續(xù)信號(hào)的MATLAB表達(dá)式。解：建模嚴(yán)格說來，MATLAB是不能表示連續(xù)信號(hào)的，因?yàn)樗o出的是各個(gè)樣本點(diǎn)的數(shù)據(jù)，只有當(dāng)樣本點(diǎn)取得很密時(shí)才可看成連續(xù)信號(hào)。所謂“密”，是相對(duì)于信號(hào)變化的快慢而言，形象地說，在相鄰樣本點(diǎn)之間的數(shù)據(jù)變化必須非常小才能看成“密”，其嚴(yán)格的數(shù)學(xué)定義此處不予討論。在編程中，先設(shè)定共同的時(shí)間坐標(biāo)，然后分別列出生成三種信號(hào)的程序。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序clear，t0=0;

tf=5;

dt=0.05;

t1=1;

t=［t0∶dt∶tf］;%(1)單位脈沖信號(hào)%在t1(t0≤t1≤tf)處有一持續(xù)時(shí)間為dt，面積為1的脈沖信號(hào)，其余時(shí)間均為零t=［t0∶dt∶tf］;

st=length(t);n1=floor((t1-t0)/dt);x1=zeros(1，st);x1(n1)=1/dt;%求t1對(duì)應(yīng)的樣本序號(hào)%把全部信號(hào)先初始化為零%給出t1處幅度為1/dt的脈沖信號(hào)subplot(2，2，1)，stairs(t，x1)%繪圖，注意為何用stairs而不用plot命令

axis(［0，5，0，1.1/dt］)%為了使脈沖頂部避開圖框，改變圖框坐標(biāo)第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例%(2)單位階躍信號(hào)%信號(hào)從t0到tf，在t1(t0≤t1≤tf)前為0，到t1處有一躍變，為1x2=［zeros(1，n1-1)，ones(1，st-n1+1)］;subplot(2，2，3)，stairs(t，x2)axis(［0，5，0，1.1］)%產(chǎn)生階躍信號(hào)%繪圖%為了使方波頂部避開圖框，改變圖框坐標(biāo)%(3)復(fù)數(shù)指數(shù)信號(hào)u=-0.5;

w=10;

x3=exp((u+j*w)*t);subplot(2，2，2)，plot(t，real(x3))

%繪圖subplot(2，2，4)，plot(t，imag(x3))

%繪圖第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序運(yùn)行結(jié)果x1、x2、x3的波形如圖9-1所示。注意：若要顯示連續(xù)信號(hào)波形中的不連續(xù)點(diǎn)，用stairs命令；若要使波形光滑些，則用plot命令較好。復(fù)數(shù)指數(shù)信號(hào)x3可以分解為余弦和正弦信號(hào)，它們分別是復(fù)數(shù)信號(hào)的實(shí)部和虛部，右圖中的兩個(gè)衰減振蕩信號(hào)就代表了這兩個(gè)相位差90°的分量。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-1例9-1-1中x1、x2、x3對(duì)應(yīng)的四種波形第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例9-1-2】編寫求任意高階連續(xù)常系數(shù)線性系統(tǒng)沖擊響應(yīng)的程序。解:這個(gè)問題在第4章4.3.5節(jié)介紹多項(xiàng)式函數(shù)庫時(shí)已經(jīng)打下基礎(chǔ)，在第7章機(jī)械振動(dòng)的例7-3-1又討論過二階常系數(shù)線性微分方程的解法，讀者可以先看懂那些例題的解法，再看本題。任意階次的連續(xù)線性系統(tǒng)可用下列線性常微分方程表述:第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例寫成傳遞函數(shù)形式為因此，其特性可以用系統(tǒng)傳遞函數(shù)的分子分母系數(shù)向量b和a來表示。向量b和a的長(zhǎng)度分別為m+1和n+1，對(duì)于一切物理上可實(shí)現(xiàn)的系統(tǒng)，必有n≥m。系統(tǒng)的沖擊響應(yīng)等于傳遞函數(shù)的拉普拉斯反變換，問題歸結(jié)為如何求出這個(gè)反變換。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例如果分母系數(shù)多項(xiàng)式?jīng)]有重根，則可以將兩個(gè)多項(xiàng)式之比分解成n個(gè)一階部分分式之和，即其中p1，p2，…，pn是分母多項(xiàng)式的n個(gè)根，而r1，r2，…，rn則是對(duì)應(yīng)于這n個(gè)根的留數(shù)。一階分式的反變換可以查表得到，所以可得到?jīng)_擊響應(yīng)的公式:第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例可見只要求出根p和留數(shù)r，線性方程的解就得到了。求根是代數(shù)問題，當(dāng)階次高時(shí)，沒有解析解?？上驳氖荕ATLAB提供了用數(shù)值方法求根和留數(shù)的residue函數(shù)，它的調(diào)用方法如下:［r，p］=residue(b，a)只要給出系數(shù)向量b和a，就可以得出根p和留數(shù)r，并得到系統(tǒng)的沖擊響應(yīng)。因此，本題的程序就可方便地編寫如下。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序a=input(′多項(xiàng)式分母系數(shù)向量a=′)

b=input(′多項(xiàng)式分子系數(shù)向量b=′)［r，p］=residue(b，a)，

%求留數(shù)disp(′沖擊響應(yīng)的解析式為h(t)=Σr(i)*exp(p(i)*t)′)k=input(′是否要求波形？是，鍵入1;否，鍵入0′);if

k==1dt=input(′dt=′);tf=input(′tf=′);

%設(shè)定時(shí)間數(shù)組t=0∶dt∶tf;h=zeros(1，length(t));%h的初始化for

i=1∶length(a)-1

%根的數(shù)目等于a的長(zhǎng)度減1h

=

h+

r(i)*exp(p(i)*t);

%疊加各根分量endplot(t，h)，gridelse，end第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序運(yùn)行結(jié)果例如，給出系統(tǒng)傳遞函數(shù)為求沖擊響應(yīng)。根據(jù)程序提問依次輸入:a=poly(［0，-1，-2，-5］)b=［1，7，1］dt=0.05tf=5得出的h(t)如圖9-2所示。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-2高階系統(tǒng)的沖擊響應(yīng)第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序中要的是系數(shù)向量a，而題中給出的是極點(diǎn)向量p=［0，-1，-2，-5］，因此這里用poly函數(shù)來作轉(zhuǎn)換。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例9-1-3】線性時(shí)不變系統(tǒng)的特性可用常系數(shù)線性微分方程表示:求輸入u為0時(shí)，由初始狀態(tài)決定的輸出，即其零輸入響應(yīng)。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例解：建模在零輸入條件下，系統(tǒng)的響應(yīng)取決于微分方程左端特征方程的根，與右端無關(guān)，其通解為其中，p1、p2、…、pn是特征方程a1λn+a2λn-1+…+anλ+an+1=0的根。每個(gè)分量的系數(shù)C1、C2、…、Cn

應(yīng)由y及其各階導(dǎo)數(shù)的初始條件來確定。(Dy0表示y的導(dǎo)數(shù)的初始值)第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例初始條件數(shù)應(yīng)該和待定系數(shù)的數(shù)目相等，構(gòu)成一個(gè)確定C1、…

Cn的線性代數(shù)方程組，寫成V*C=Y0其解為

C=V-1Y0其中C=［C1，C2，…，Cn］′;Y0=［y0，Dy0，…，Dny0］′這種形式的矩陣稱為范德蒙特矩陣。在MATLAB的特殊矩陣庫中提供了專門的函數(shù)vander，給出p向量，就可由V=vander(p)生成范德蒙特矩陣。從help

vander知道，需要將它旋轉(zhuǎn)90°，才與本題的形式相符。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序a=input(′輸入分母系數(shù)向量a=［a1，a2，...］=′);

n=length(a)-1;Y0=input(′輸入初始條件向量Y0=[y0,Dy0,D2y0，...]=

′);p=roots(a);V=rot90(vander(p));

%求根，構(gòu)成范德蒙特矩陣C=V＼Y0′;

%求出Cdt=input(′dt=′);

tf=input(′tf=

′)t=0∶dt∶tf;y=zeros(1，length(t));%給出自變量數(shù)據(jù)組for

k=1∶n

y=y+c(k)*exp(p(k)*t);end%求各分量的時(shí)間函數(shù)并疊加plot(t，y)，grid第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例下面利用這一程序來解一個(gè)三階系統(tǒng)。程序運(yùn)行結(jié)果運(yùn)行此程序并輸入：a=［1，2，9，3］;dt=0.1;tf=5;而初始值Y0分別?。?，0，0］;

［0，1，0］;

［0，0，1］用hold

on語句使三次運(yùn)行生成的圖形畫在一幅圖上，得到圖9-3。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-3三階系統(tǒng)零輸入分量的解第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例

9-1-4】

n級(jí)放大器，每級(jí)的轉(zhuǎn)移函數(shù)均為

，求階躍輸入下的過渡過程，畫出n不同時(shí)的波形及頻率特性進(jìn)行比較。解：建模系統(tǒng)的轉(zhuǎn)移函數(shù)為

，階躍輸入的拉普拉斯變換為

，因此輸出為兩者的乘積，即求Y(s)的拉普拉斯反變換，即可得到輸出過渡過程y(t)。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例這里我們遇到了一個(gè)有多重極點(diǎn)-ωn的H(s)求拉普拉斯反變換的問題。從原理上說，同樣可以用4.3節(jié)中的留數(shù)極點(diǎn)分解法來求。只是在有n重根時(shí)，分解出的部分分式的分母將不再是一次極點(diǎn)，而有(s+ωn)、…、(s+ωn)n-1、(s+ωn)n等項(xiàng)，而(s+ωn)-q的反變換可用如下解析式表示:第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例按照這個(gè)思路，應(yīng)該先求出Y(s)的極點(diǎn)留數(shù)。注意分母中除了有n個(gè)重極點(diǎn)外，還有一個(gè)零極點(diǎn)（即1/s），故共有n+1個(gè)極點(diǎn)。先用poly函數(shù)求出Y(s)分母多項(xiàng)式的系數(shù)向量，并求出極點(diǎn)及相應(yīng)的留數(shù)：by=wn^n;

ay=［poly(-ones(1，n)*wn)，0］［r，p］=residue(by，ay);再把各個(gè)分量疊加在一起。然而，實(shí)際上，這樣編程不僅非常麻煩，而且難以得到正確的結(jié)果。其原因是在重極點(diǎn)處，MATLAB

的residue算法遇到了病態(tài)問題，數(shù)據(jù)中小小的舍入誤差會(huì)使結(jié)果產(chǎn)生很大誤差，即使n取2都得不出正確結(jié)果。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例為了避開重極點(diǎn)問題，可以有意把極點(diǎn)拉開一些。例如，設(shè)n個(gè)極點(diǎn)散布在0.98ωn～1.02ωn之間，那樣就可以全部作為非重

極點(diǎn)來列程序。這種處理在工程上是完全沒有問題的，一般電阻的標(biāo)稱誤差為±5%，電容則更大，在工程實(shí)踐中，使各個(gè)放大器時(shí)常數(shù)完全相同是不可能的，即便要把其誤差控制到±2%以內(nèi)也非易事，所以不必自找麻煩，干脆用非重極點(diǎn)的下列程序來求解。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序clear，clf，N=input(′輸入放大器級(jí)數(shù)N=′);wn=1000;dt=1e-4;tf=0.01;t=0∶dt∶tf;y=zeros(N，length(t));for

n=1∶Np0=-linspace(.95，1.05，n)*wn;%輸出初始化%將H(s)極點(diǎn)分散布置在±5%區(qū)間

ay=poly(［p0，0］);%由Y(s)的極點(diǎn)（比H(s)多一個(gè)零極點(diǎn)）求分母系數(shù)by=prod(abs(p0));

%求Y(s)的分子系數(shù)［r，p］=residue(by，ay);

%求Y(s)的留數(shù)極點(diǎn)第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例for

k=1∶n+1

%把各部分分式對(duì)應(yīng)的時(shí)域分量相加y(n,∶)=y(n,∶)+r(k)*exp(p(t);endfigure(1)，plot(t，y(n,∶));grid，hold

on%繪制過渡過程曲%下面這幾條語句用來繪制波特圖figure(2)，bode(prod(abs(p0))，poly(p0));hold

onbh=by;

ah=poly(p0);w=logspace(2，4);%求H(s)的分子分母系數(shù)%給出頻率范圍和分度H=polyval(bh，j*w)./polyval(ah，j*w);%求H(s)在各頻點(diǎn)的值H(jw)aH=unwrap(angle(H))*180/pi;%求出以度為單位的連續(xù)相角第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例fH=20*log10(abs(H));

%求出以分貝為單位的振幅figure(2)，subplot(2，1，1)，semilogx(w，fH)，grid

on，hold

on%繪幅頻圖subplot(2，1，2)，semilogx(w，aH)，grid

on，hold

on%繪相頻圖end，hold

off第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序運(yùn)行結(jié)果運(yùn)行此程序，設(shè)N=4，可得1級(jí)到4級(jí)放大器的過渡過程如圖9-4所示，從中看出輸出信號(hào)達(dá)到0.6處所需的時(shí)間約為單級(jí)時(shí)常數(shù)乘以級(jí)數(shù)。此程序在N>4時(shí)又會(huì)出現(xiàn)很大誤差，讀者可自己編寫更好的程序。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-4多級(jí)放大器的階躍過渡過程第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例為了畫波特圖，程序可按以下步驟進(jìn)行：求頻率特性。用多項(xiàng)式求值函數(shù)polyval，并且用了元素群運(yùn)算，輸入頻率數(shù)組作為自變量，一次就求出全部的頻率特性。注意頻率特性是復(fù)數(shù)，通常關(guān)心的是它們的振幅和相角。振幅和相位特性橫坐標(biāo)都要用對(duì)數(shù)坐標(biāo)并且上下對(duì)齊。振幅的縱坐標(biāo)單位應(yīng)化為分貝,這里用了20*log10(abs(H))相位的縱坐標(biāo)單位應(yīng)為度，并且應(yīng)連續(xù)變化，不取主角，所以這里加了unwrap命令。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例在控制系統(tǒng)工具箱中有一個(gè)bode命令可以直接完成這些功能，但本書遵循的原則是不用工具箱，以便讓讀者知道工具箱是怎么編程的。當(dāng)然讀者也應(yīng)知道，工具箱中的bode函數(shù)，遠(yuǎn)不止這幾條語句，作為商用軟件庫中的一個(gè)正式函數(shù)，它必須考慮自動(dòng)確定頻率區(qū)間，自動(dòng)檢查輸入有無錯(cuò)誤等，程序要復(fù)雜得多。圖9-5繪出了多級(jí)放大器的頻率特性。其幅頻特性顯示了低通的特點(diǎn)，隨級(jí)數(shù)的增加，通帶減小；相頻特性說明，隨級(jí)數(shù)的增加，負(fù)相移成比例地增加。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-5多級(jí)放大器的頻率特性第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例9-1-5】設(shè)方波信號(hào)的寬度為5

s，信號(hào)持續(xù)期為10

s，試求其在0～20(1/s)頻段間的頻譜特性。如只取0～10(1/s)的頻譜分量（相當(dāng)于通過了一個(gè)低通濾波器），求其輸出波形。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例解：建模設(shè)信號(hào)的時(shí)域波形為f(t)，在0~10

s的區(qū)間外信號(hào)為0，則其傅里葉變換為按MATLAB作數(shù)值計(jì)算的要求，必須把t分成N份，用相加來代替積分，對(duì)于任一給定的ω，可寫成這說明求和的問題可以用f(t)行向量乘以ejωt列向量來實(shí)現(xiàn)。此處的Δt是t的增量，在程序中，將用dt來代替。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例由于要求出一系列不同的ω處的F值都用同一公式，因此可以利用MATLAB中的元素群運(yùn)算功能，把ω設(shè)成一個(gè)行數(shù)組，分別

代入本公式左右端的ω中去，寫成(程序中把ω寫成w)F=f*exp(-j*t′*w)·Δt其中，F(xiàn)是與w等長(zhǎng)的行向量，exp中的t′是列向量，w是行向量，t′*w是一個(gè)矩陣，其行數(shù)與t相同，列數(shù)與w相同。這個(gè)矩陣乘式就完成了傅里葉變換。類似地可以得到傅里葉逆變換表示式。由此得到下面的傅里葉變換程序。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序

clear，tf=10;N=256;t=linspace(0，tf，N);%給出時(shí)間分割w1=linspace(eps，20，N);

dw=20/(N-1);%dw=1/4/tf;w1=［eps∶dw∶(N-1)/4/tf］;

%給出頻率分割f=［ones(1，N/2)，zeros(1，N/2)］;%給出信號(hào)(此處是方波F1=f*exp(-j*t′*w1)*tf/(N-1);w=［-fliplr(w1)，w1(2∶N)］;%求傅里葉變換%補(bǔ)上負(fù)頻率第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例F=［fliplr(F1)，F(xiàn)1(2∶N)］;%補(bǔ)上負(fù)頻率區(qū)的頻譜w2=w(N/2∶3*N/2);F2=F(N/2∶3*N/2);%取出中段頻率%取出中段頻譜subplot(1，2，1)，plot(w，abs(F)，′linewidth′，1.5)，grf1=F2*exp(j*w2′*t)/pi*dw;

%對(duì)中段頻譜求傅里葉逆變換subplot(1，2，2)，plot(t，f，t，f1，′linewidth′，1.5)，第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序運(yùn)行結(jié)果執(zhí)行這個(gè)程序的結(jié)果如圖9-6所示，因?yàn)榉讲ê泻茇S富的高頻分量，要充分恢復(fù)其原來波形需要很寬的頻帶，實(shí)踐中不可能完全做到。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-6方波信號(hào)的頻譜和取|ω|＜10部分頻譜的逆變換波形(a)方波信號(hào)的頻譜;(b)取|ω|<10部分頻譜的逆變換波形第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例9.2

離散信號(hào)和系統(tǒng)信號(hào)可以分為模擬信號(hào)和數(shù)字信號(hào)。模擬信號(hào)用x(t)表示，其中變量t代表時(shí)間。離散信號(hào)用x(n)表示，其中變量n為整數(shù)并代表時(shí)間的離散時(shí)刻，因此它也稱為離散時(shí)間信號(hào)。離散信號(hào)是一個(gè)數(shù)字的序列，并可以表述為其中，向上的箭頭表示在n=0處的取樣。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例在MATLAB中，可以用一個(gè)向量x來表示一個(gè)有限長(zhǎng)度的序

列。然而這樣一個(gè)向量并沒有包含基準(zhǔn)采樣位置的信息。因此，完全地表示x(n)要用x和n兩個(gè)向量。例如序列x(n)={2，1，-1，

，1，4，3，7}(下面的箭頭為第0個(gè)采樣點(diǎn))，在

MATLAB中表示為n

=

[-3,-2,-1,0,1,2,3,4]，

x=[2,1,-1,5,1,4,3,7]當(dāng)不需要采樣位置信息時(shí)，可以只用x向量來表示。由于內(nèi)存有限，因此MATLAB無法表示無限序列。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例9-2-1】編寫MATLAB程序來產(chǎn)生下列基本脈沖序列。單位脈沖序列：起點(diǎn)n0，終點(diǎn)nf，在ns處有一單位脈沖(n0≤ns≤nf)。單位階躍序列：起點(diǎn)n0，終點(diǎn)nf，在ns前為0，在ns處及以后為1(n0≤ns≤nf)。實(shí)數(shù)指數(shù)序列:x3=(0.9)n。復(fù)數(shù)指數(shù)序列:x4=e(-0.2+0.3j)n。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例解:建模這些基本序列的表達(dá)式比較簡(jiǎn)明，編寫程序也不難。對(duì)單位脈沖序列，我們提供了直接賦值和邏輯關(guān)系兩種方法,其中用邏輯關(guān)系的編法比較簡(jiǎn)潔，讀者從中可看到MATLAB編程的靈活性和技巧性。通常用stem語句來繪制離散序列。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序

clear，no=0;nf=10;ns=3;n1=n0∶nf;

x1=［zeros(1，ns-n0)，1，zeros(1，nf-ns)］;%n1

=

n0∶nf;

x1=［(n1-ns)==0］;%顯然，用邏輯式是比較高明的方法n2=n0∶nf;x2=［zeros(1，ns-n0)，ones(1，nf-ns+1)］;%也有類似的用邏輯比較語句產(chǎn)生單位階躍序列的方法，留給讀者思考n3=n0∶nf;x3=(0.9).^n3;

%實(shí)數(shù)指數(shù)序列n4=n0∶nf;x4=exp((-0.2+0.3j)*n3);%復(fù)數(shù)指數(shù)序列第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例subplot(2，2，1)，stem(n1，x1);subplot(2，2，2)，stem(n2，x2);subplot(2，2，3)，stem(n3，x3);subplot(4，2，6)，stem(n4，real(x4));%注意subplot的輸入變?cè)猻ubplot(4，2，8)，stem(n4，imag(x4));line(［0，10］，［0，0］)，程序運(yùn)行結(jié)果程序運(yùn)行結(jié)果如圖9-7所示。%畫橫坐標(biāo)第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-7基本脈沖序列的波形第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例9-2-2】離散傅里葉變換的計(jì)算。解：建模一個(gè)時(shí)間序列x(n)的離散時(shí)間傅里葉變換的定義為如果序列的長(zhǎng)度是有限的，可以把它看做是周期性無限序列中的一個(gè)周期，其長(zhǎng)度為N。對(duì)這個(gè)周期性序列可以用離散傅里葉變換（注意少了“時(shí)間”兩字）進(jìn)行研究，它的定義為其中第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例用例9-1-4中的方法，引入矩陣乘法來實(shí)現(xiàn)求和運(yùn)算。用元素群算法來求不同k時(shí)的X，把n和k都設(shè)成1×N的行數(shù)組，令nk=n′*k，它就成為N×N的方陣，因而

也是N×N方陣。由此得出離散傅里葉變換的算式MATLAB只能處理有限長(zhǎng)度的序列，因此，適合于計(jì)算離散傅里葉變換及其逆變換。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序設(shè)有限信號(hào)序列xn(n)的長(zhǎng)度為Nx，則按定義，求其N點(diǎn)傅里葉變換Xk(k)的程

序?yàn)?xn=input(′x=′);Nx=length(xn);N=Nx

%取N為x的長(zhǎng)度tic，n=［0∶1∶N-1］;k=［0∶1∶N-1］;%設(shè)定n和k的行向量%WN因子%產(chǎn)生一個(gè)含nk值的N×N維矩陣%換算矩陣WN=exp(-j*2*pi/N);nk=n′*k;WNnk=WN

.^

nk;Xk=xn

*

WNnk;

toc%DFT

系數(shù)向量，即離散傅里葉變換的結(jié)果plot(abs(Xk))，grid

%繪幅頻特性圖第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例在N很大時(shí)，這個(gè)程序的運(yùn)算速度比較低。程序中的tic和toc語句用來測(cè)試它們之間的程序運(yùn)行時(shí)間。實(shí)際上MATLAB已提供了快速離散傅里葉變換的函數(shù)fft，可直接調(diào)用。其調(diào)用格式為X

=

fft(x，N)其中，x是輸入的時(shí)間序列，N是傅里葉變換取的點(diǎn)數(shù)。若省略N，則它自動(dòng)把x的長(zhǎng)度作為N。當(dāng)N取2的冪時(shí)，變換速度最快，所以要提高fft函數(shù)的運(yùn)行速度，程序應(yīng)編寫如下：xn=input(′x

=

′);

Nx=length(xn)%取N為大于Nx而最接近于Nx的2的冪

N=pow2(nextpow2(Nx);tic，X=fft(xn，N);tocNx＜N，x長(zhǎng)度不足N的部分，程序會(huì)自動(dòng)補(bǔ)0。要注意X是一

個(gè)長(zhǎng)度為N的復(fù)數(shù)數(shù)組，可以分解出它的振幅和相位，分別繪圖。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序運(yùn)行結(jié)果按程序提示輸入x=sin(0.1*［1∶700］)+randn(1，700);所得的fft的幅度特性是一樣的，如圖9-8所示，其中有效信號(hào)與噪聲可明顯區(qū)分。在作者的計(jì)算機(jī)上，前一程序的時(shí)間測(cè)試結(jié)果約為8

s，后一程序?yàn)?.05

s。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-8信號(hào)的fft的振幅頻率特性第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例9.3

系統(tǒng)函數(shù)【例9-3-1】簡(jiǎn)單信號(hào)流圖模型的矩陣解法[11]。建模信號(hào)流圖是用來表示和分析復(fù)雜系統(tǒng)內(nèi)的信號(hào)變換關(guān)系的工具。其基本概念如下:系統(tǒng)中每個(gè)信號(hào)用圖上的一個(gè)節(jié)點(diǎn)表示。如圖中的u、x1、x2。系統(tǒng)部件對(duì)信號(hào)實(shí)施的變換關(guān)系用有向線段表示，箭尾為輸入信號(hào)，箭頭為輸出信號(hào)，箭身標(biāo)注對(duì)此信號(hào)進(jìn)行變換的乘子。如圖上的G1、G2。如果乘子為1，可以不必標(biāo)注。每個(gè)節(jié)點(diǎn)信號(hào)的值等于所有指向此節(jié)點(diǎn)的箭頭信號(hào)之和，每個(gè)節(jié)點(diǎn)信號(hào)可以向外輸出給多個(gè)部件，其值相同。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例根據(jù)這幾個(gè)概念，可以列出圖9-9的方程如下：x1=u-G2x2，

x2=G1x1第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-9帶反饋的簡(jiǎn)單信號(hào)流圖第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例寫成矩陣方程或x=Qx+Pu移項(xiàng)整理，可以得到求所有未知信號(hào)向量x的公式第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例定義系統(tǒng)的傳遞函數(shù)W為輸出信號(hào)與輸入信號(hào)之比x/u，則W可按下式求得:W=x/u=inv(I-Q)*P因?yàn)榍筮@個(gè)二階矩陣的逆可以直接用下面的公式：若則第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例所以即u到x1的傳遞函數(shù)為

，u到x2的傳遞函數(shù)為。對(duì)于階次高的情況，求逆就必須用軟件工具了。如果信號(hào)流圖中有G1那樣的符號(hào)變量,那么它的求解要用符號(hào)運(yùn)算工具箱。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序

syms

G1

G2Q=［0，-G2;

G1，0］，P=［1;

0］W=inv(eye(2)-Q)*P第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序運(yùn)行結(jié)果程序運(yùn)行的結(jié)果是與前面的結(jié)果相同。這是一個(gè)簡(jiǎn)單的問題，用一些其他的數(shù)學(xué)方法也能得到同樣的結(jié)果，很出名的“梅森公式”就是用圖形拓樸的方法得到信號(hào)流圖的公式，但這個(gè)公式非常繁瑣，而且無法機(jī)械化。到現(xiàn)在為止，我們還沒有見過任何一本書籍用矩陣方法來推導(dǎo)這個(gè)公式。用矩陣代數(shù)方法的最大好處是可以向任意高的階次、任意復(fù)雜的信號(hào)流圖推廣，實(shí)現(xiàn)復(fù)雜系統(tǒng)傳遞函數(shù)推導(dǎo)的自動(dòng)化[11-12]。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例9-3-2】較復(fù)雜信號(hào)流圖模型的矩陣解法。建模圖9-10是一個(gè)較復(fù)雜的信號(hào)流圖。照上述方法列出它的方程如下:x1=-G4x3+ux2=G1x1-G5x4x3=G2x2x4=G3x3第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例將其列為矩陣方程，得到公式W=x/u=inv(I-Q)·P同樣是正確的，不過這里的Q和P分別為4×4和4×1矩陣，用手工求逆是很麻煩的。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-10帶雙重反饋的信號(hào)流圖第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例MATLAB程序

syms

G1

G2

G3

G4

G5Q=［0，0，-G4，0;

G1，0，0，-G5;

0，G2，0，0;

0，0，G3，0］，P=［1;

0;

0;

0］W=inv(eye(4)-Q)*PPretty(W(4))第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例我們關(guān)心的輸出通常是x4，也就是最后那個(gè)傳遞函數(shù)W(4)＝x4)/u，其結(jié)果為當(dāng)系統(tǒng)內(nèi)各個(gè)環(huán)節(jié)都是線性集總參數(shù)，因而它們的傳遞函數(shù)Gi都可以表示為s的多項(xiàng)式有理分式

時(shí)，不管其階次有多高，傳遞函數(shù)W都可以很容易地由計(jì)算機(jī)直接自動(dòng)算出。這個(gè)方法還可以推廣到離散系統(tǒng)，用來計(jì)算任意復(fù)雜的數(shù)字濾波器系統(tǒng)函數(shù)。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例9.4頻譜及其幾何意義頻譜分析是信號(hào)與系統(tǒng)課程中最重要的內(nèi)容之一,很多讀者在學(xué)習(xí)中感到這部分內(nèi)容很抽象,往往只能從數(shù)學(xué)上了解時(shí)域信號(hào)與其頻譜間的變換關(guān)系,而沒有理解它的物理意義,而用MATLAB可以幫助讀者建立形象的幾何概念,真正掌握它。首先來看歐拉公式,它以最簡(jiǎn)明的方式建立了信號(hào)頻域與時(shí)域的關(guān)系:第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例它說明一個(gè)最簡(jiǎn)單的實(shí)余弦信號(hào)可以由正、負(fù)兩個(gè)Ω0頻率分量合成。在復(fù)平面上,正的Ω0對(duì)應(yīng)反時(shí)針旋轉(zhuǎn)的向

量,負(fù)的Ω0對(duì)應(yīng)順時(shí)針旋轉(zhuǎn)的向量,當(dāng)這兩個(gè)向量的幅度相同而相角符號(hào)相反時(shí),就合成為一個(gè)在實(shí)軸上的向量。它的相角為零,大小按正弦變化,形成了實(shí)信號(hào)cosΩ0t(如圖9-11所示)。推而廣之,任何實(shí)周期信號(hào)必然具有正、負(fù)兩組頻率的頻譜成分,正、負(fù)頻率頻譜的幅度對(duì)稱而相位反對(duì)稱,或者說正、負(fù)頻率的頻譜是共軛的。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-11實(shí)序列由對(duì)稱的正負(fù)頻率合成第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例【例9-4-1】設(shè)計(jì)一個(gè)演示程序,要求它能把用戶任意給定的四個(gè)集總頻譜合成,并生成對(duì)應(yīng)的時(shí)域信號(hào)。解:建模按上述多節(jié)桿合成模型來考慮這個(gè)問題。程序設(shè)計(jì)主要包括三部分:(1)各頻譜分量的輸入,包括其幅度和頻率(有正負(fù)號(hào));(2)將各分量當(dāng)做轉(zhuǎn)動(dòng)的桿件首尾相接;(3)記錄多節(jié)桿系末端的軌跡,并畫出圖形。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例程序運(yùn)行結(jié)果運(yùn)行程序并按提示輸入數(shù)據(jù),如果只取兩個(gè)幅度相等而頻率符號(hào)相反的集總頻譜,那么將得到與圖9-11相仿的結(jié)果?，F(xiàn)取四個(gè)集總頻譜,輸入a(1)=1,w(1)=-1;a(2)=1,w(2)=-1;a(3)=0.5,w(3)=3;a(4)=0.5,w(1)=-4;第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-12四個(gè)頻譜向量組成的多節(jié)桿及其端點(diǎn)軌跡第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例圖9-13四個(gè)集總頻譜生成的復(fù)信號(hào)及其實(shí)信號(hào)分量(a)合成的復(fù)平面信號(hào);(b)合成的虛部信號(hào);(c)合成的實(shí)部信號(hào)第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例輸入頻譜的幅度可以是負(fù)數(shù),也可以是虛數(shù),甚至可以是復(fù)數(shù),它不僅反映了頻譜的大小,還反映了該向量的起始相位;頻譜的頻率則只能是可正可負(fù)的實(shí)數(shù),正頻率和負(fù)頻率以及在該頻率上頻譜的意義在此不言自明。讀者可以做各種各樣的試驗(yàn),例如當(dāng)兩組頻率具有倍頻關(guān)系時(shí),得到的是周期信號(hào);如果頻率比是任意小數(shù),那將得出非周期的信號(hào)。另外,這樣的演示適用于集總頻譜,對(duì)于分布的頻譜密度,就要把它想象為若干小的集總頻譜的疊合?？傊?有了這樣的形象演示,便可以大大擴(kuò)展時(shí)域信號(hào)與頻域譜之間關(guān)系的思維空間。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例可見,只有在復(fù)信號(hào)平面上,才能看到頻率的正負(fù)。許多人往往看不到負(fù)頻率頻譜的意義,這是因?yàn)樗麄兛偸峭Ａ粼趯?shí)信號(hào)的范疇來思考問題。要想知道象有沒有鼻子(負(fù)頻率),必須去摸象頭(復(fù)信號(hào)),只摸象腿(實(shí)信號(hào))是得不到真切的認(rèn)識(shí)(信號(hào)理論普遍而科學(xué)的規(guī)則)的。第9章

MATLAB在信號(hào)和系統(tǒng)中的應(yīng)用舉例9.5場(chǎng)的計(jì)算和偏微分方程的數(shù)值解現(xiàn)代科學(xué)和工程中大量遇到偏微分方程,這類方程通常用來表示分布在三維空間的場(chǎng)。它的解析解是很難求得的,特別是在復(fù)雜的邊界條件和初始條件下,不太可能有解析解。實(shí)際工程上只要得到它的數(shù)值解就可以了。此時(shí)的一個(gè)重要手段就是把微分變?yōu)橛邢薏罘?把積分變?yōu)榍蠛?從而把偏微分方程變?yōu)榇鷶?shù)方程組。要得出比較準(zhǔn)確的結(jié)果,差分必須取得很小,因而所得到的代數(shù)方程組的階數(shù)可能很高,沒有極高水平的計(jì)算機(jī)和算法,是無法得出有