. 信號與系統及MATLAB實現實驗指導書前 言長期以來，信號與系統課程一直采用單一理論教學方式，同學們依靠做習題來鞏固和理解教學內容，雖然手工演算訓練了計算能力和思維方法，但是由于本課程數學公式推導較多，概念抽象，常需畫各種波形，作題時難免花費很多時間，現在，我們給同學們介紹一種國際上公認的優秀科技應用軟件MATLAB，借助它我們可以在電腦上輕松地完成許多習題的演算和波形的繪制。 MATLAB的功能非常強大，我們此處僅用到它的一部分，在后續課程中我們還會用到它，在未來地科學研究和工程設計中有可能繼續用它，所以有興趣的同學，可以對MATLAB再多了解一些。MATLAB究竟有那些特點呢？1高效的數值計算和符號計算功能，使我們從繁雜的數學運算分析中解脫出來；2完備的圖形處理功能，實現計算結果和編程的可視化；3友好的用戶界面及接近數學表達式的自然化語言，易于學習和掌握；4功能豐富的應用工具箱，為我們提供了大量方便實用的處理工具；MATLAB的這些特點，深受大家歡迎，由于個人電腦地普及，目前許多學校已將它做為本科生必須掌握的一種軟件。正是基于這些背景，我們編寫了這本信號與系統及MATLAB實現指導書，內容包括信號的MATLAB表示、基本運算、系統的時域分析、頻域分析、S域分析、狀態變量分析等。通過這些練習，同學們在學習信號與系統的同時，掌握MATLAB的基本應用，學會應用MATLAB的數值計算和符號計算功能，擺脫煩瑣的數學運算，從而更注重于信號與系統的基本分析方法和應用的理解與思考，將課程的重點、難點及部分習題用MATLAB進行形象、直觀的可視化計算機模擬與仿真實現，加深對信號與系統的基本原理、方法及應用的理解，為學習后續課程打好基礎。另外同學們在進行實驗時，最好事先預習一些MATLAB的有關知識，以便更好地完成實驗，同時實驗中也可利用MATLAB的help命令了解具體語句以及指令的使用方法。實驗一 基本信號在MATLAB中的表示和運算一、實驗目的1 學會用MATLAB表示常用連續信號的方法；2 學會用MATLAB進行信號基本運算的方法；二、實驗原理1 連續信號的MATLAB表示MATLAB提供了大量的生成基本信號的函數，例如指數信號、正余弦信號。表示連續時間信號有兩種方法，一是數值法，二是符號法。數值法是定義某一時間范圍和取樣時間間隔，然后調用該函數計算這些點的函數值，得到兩組數值矢量，可用繪圖語句畫出其波形；符號法是利用MATLAB的符號運算功能，需定義符號變量和符號函數，運算結果是符號表達的解析式，也可用繪圖語句畫出其波形圖。例1-1指數信號 指數信號在MATLAB中用exp函數表示。如，調用格式為 ft=A*exp(a*t) 程序是 A=1; a=-0.4;t=0:0.01:10; %定義時間點ft=A*exp(a*t); %計算這些點的函數值plot(t,ft); %畫圖命令，用直線段連接函數值表示曲線grid on; %在圖上畫方格例1-2 正弦信號 正弦信號在MATLAB中用 sin 函數表示。調用格式為 ft=A*sin(w*t+phi)A=1; w=2*pi; phi=pi/6;t=0:0.01:8; %定義時間點ft=A*sin(w*t+phi); %計算這些點的函數值plot(t,ft); %畫圖命令grid on; %在圖上畫方格例1-3 抽樣信號 抽樣信號Sa(t)=sin(t)/t在MATLAB中用 sinc 函數表示。定義為 t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft); grid on; axis(-10,10,-0.5,1.2); %定義畫圖范圍，橫軸，縱軸 title(抽樣信號) %定義圖的標題名字例1-4 三角信號 三角信號在MATLAB中用 tripuls 函數表示。調用格式為 ft=tripuls(t,width,skew)，產生幅度為1，寬度為width，且以0為中心左右各展開width/2大小，斜度為skew的三角波。width的默認值是1，skew的取值范圍是-1+1之間。一般最大幅度1出現在t=(width/2)*skew的橫坐標位置。 t=-3:0.01:3; ft=tripuls(t,4,0.5); plot(t,ft); grid on; axis(-3,3,-0.5,1.5);例1-5 虛指數信號 調用格式是f=exp(j*w)*t) t=0:0.01:15;w=pi/4;X=exp(j*w*t);Xr=real(X); %取實部 Xi=imag(X); %取虛部Xa=abs(X); %取模Xn=angle(X); %取相位subplot(2,2,1),plot(t,Xr),axis(0,15,-(max(Xa)+0.5),max(Xa)+0.5),title(實部);subplot(2,2,3),plot(t,Xi),axis(0,15,-(max(Xa)+0.5),max(Xa)+0.5),title(虛部);subplot(2,2,2), plot(t,Xa),axis(0,15,0,max(Xa)+1),title(模);subplot(2,2,4),plot(t,Xn),axis(0,15,-(max(Xn)+1),max(Xn)+1),title(相角); %subplot(m,n,i) 命令是建立m行n列畫圖窗口，并指定畫圖位置i例1-6 復指數信號 調用格式是f=exp(a+j*b)*t) t=0:0.01:3;a=-1;b=10;f=exp(a+j*b)*t); subplot(2,2,1),plot(t,real(f),title(實部)subplot(2,2,3),plot(t,imag(f),title(虛部) subplot(2,2,2),plot(t,abs(f),title(模)subplot(2,2,4),plot(t,angle(f),title(相角)例1-7 矩形脈沖信號 矩形脈沖信號可用rectpuls函數產生，調用格式為y=rectpuls(t,width),幅度是1，寬度是width，以t=0為對稱中心。t=-2:0.01:2;width=1;ft=2*rectpuls(t,width);plot(t,ft)grid on;例1-8 單位階躍信號 單位階躍信號u(t)用“t=0”產生，調用格式為ft=(t=0) t=-1:0.01:5; ft=(t=0); plot(t,ft); grid on; axis(-1,5,-0.5,1.5);例1-9 正弦信號符號算法syms t %定義符號變量ty=sin(pi/4*t) %符號函數表達式ezplot(y,-16,16) %符號函數畫圖命令或者f=sym(sin(pi/4*t) %定義符號函數表達式ezplot(f,-16,16)例1-10 單位階躍信號 MATTLAB符號數學函數Heaviside表示階躍信號，但要畫圖需在工作目錄創建Heaviside的M文件function f=Heaviside(t)f=(t0); 保存，文件名是Heaviside ，調用該函數即可畫圖，例t=-1:0.01:3;f=heaviside(t);plot(t,f)axis(-1,3,-0.2,1.2) 或者y=sym(Heaviside(t);ezplot(y,-1,5);grid on 2 信號基本運算的MATLAB實現 信號基本運算是乘法、加法、尺度、反轉、平移、微分、積分，實現方法有數值法和符號法例1-11 以f(t)為三角信號為例，求f(2t) , f(2-2t) t=-3:0.001:3; ft=tripuls(t,4,0.5); subplot(3,1,1); plot(t,ft); grid on; title (f(t); ft1= tripuls(2*t,4,0.5); subplot(3,1,2); plot(t,ft1); grid on; title (f(2t); ft2= tripuls(2-2*t,4,0.5); subplot(3,1,3); plot(t,ft2); grid on; title (f(2-2t);例1-12 已知f1(t)=sinwt , f2(t)=sin8wt , w=2pi , 求f1(t)+f2(t)和f1(t)f2(t) 的波形圖 w=2*pi;t=0:0.01:3;f1=sin(w*t);f2=sin(8*w*t);subplot(211)plot(t,f1+1,:,t,f1-1,:,t,f1+f2)grid on,title(f1(t)+f2(t)subplot(212)plot(t,f1,:,t,-f1,:,t,f1.*f2)grid on,title(f1(t)*f2(t)符號算法也可實現上述運算，以信號的微積分運算為例說明符號算法應用微分的調用格式為 diff(function,variable,n)積分的調用格式為 int(function,variable,a,b)式中function表示要微分或積分的函數，variable表示運算變量，n表示求導階數，默認值是求一階導數，a是積分下限，b是積分上限，a b默認是求不定積分。例1-13 求一階導數的例題，已知，clear syms a x y1 y2 %定義符號變量a， x ，y1， y2 y1=sin(a*x2); %符號函數y1 y2=x*sin(x)*log(x); %符號函數y2 dy1=diff(y1,x) %無分號直接顯示結果 dy2=diff(y2) %無分號直接顯示結果例1-14 求積分的例題， clearsyms a x y3 y4y3=x5-a*x2+sqrt(x)/2;y4=(x*exp(x)/(1+x)2;iy3=int(y3,x)iy4=int(y4,0,1)三、上機實驗內容1 驗證實驗原理中程序2 畫出信號波形（1） （2）3信號，求、波形實驗二 離散信號與系統的時域分析一、實驗目的1學會用MATLAB表示常用離散信號的方法；2學會用MATLAB實現離散信號卷積的方法；3 學會用MATLAB求解離散系統的單位響應；4 學會用MATLAB求解離散系統的零狀態響應；二、實驗原理1離散信號的MATLAB表示表示離散時間信號f(k)需要兩個行向量，一個是表示序號k= ，一個是表示相應函數值f= ，畫圖命令是stem。例2-1正弦序列信號 正弦序列信號可直接調用MATLAB函數cos，例，當是整數或分數時，才是周期信號。畫，波形程序是：k=0:40;subplot(2,1,1)stem(k,cos(k*pi/8),filled)title(cos(k*pi/8)subplot(2,1,2)stem(k,cos(2*k),filled)title(cos(2*k) 例2-2 單位序列信號 本題先建立一個畫單位序列的M函數文件，畫圖時調用。M文件建立方法：file / new / m-file 在文件編輯窗輸入程序，保存文件名用函數名。function dwxulie(k1,k2,k0) % k1 , k2 是畫圖時間范圍，k0是脈沖位置k=k1:k2;n=length(k);f=zeros(1,n);f(1,-k0-k1+1)=1; stem(k,f,filled) axis(k1,k2,0,1.5)title(單位序列(k) 保存文件名dwxulie.m畫圖時在命令窗口調用，例：dwxulie(-5,5,0)例2-3 單位階躍序列信號 本題也可先建立一個畫單位階躍序列的M函數文件，畫圖時調用。 function jyxulie(k1,k2,k0)k=k1:-k0-1;kk=-k0:k2;n=length(k);nn=length(kk)u=zeros(1,n);uu=ones(1,nn); stem(kk,uu,filled)hold onstem(k,u,filled)hold offtitle(單位階躍序列)axis(k1 k2 0 1.5)保存文件名jyxulie.m畫圖時在命令窗口調用，例：jyxulie(-3,8,0)例2-4 實指數序列信號 ，c、 a是實數。建立一個畫實指數序列的M函數文件，畫圖時調用。function dszsu(c,a,k1,k2)%c：指數序列的幅度%a：指數序列的底數%k1：繪制序列的起始序號%k2：繪制序列的終止序號k=k1:k2;x=c*(a.k);stem(k,x,filled)hold onplot(k1,k2,0,0)hold off調用該函數畫信號： ，波形。dszsu(1,5/4,0,40)dszsu(1,-3/4,0,40)2 .離散信號的卷積和兩個有限長序列f1，f2卷積可調用MATLAB函數conv，調用格式是f=conv(f1,f2), f是卷積結果，但不顯示時間序號，可自編一個函數dconv給出f和k，并畫圖。function f,k=dconv(f1,f2,k1,k2) %The function of compute f=f1*f2% f: 卷積和序列f(k)對應的非零樣值向量% k： 序列f(k)的對應序號向量% f1: 序列f1(k)非零樣值向量% f2: 序列f2(k)的非零樣值向量% k1: 序列f1(k)的對應序號向量% k2: 序列f2(k)的對應序號向量f=conv(f1,f2) %計算序列f1與f2的卷積和fk0=k1(1)+k2(1);%計算序列f非零樣值的起點位置k3=length(f1)+length(f2)-2;%計算卷積和f的非零樣值的寬度k=k0:k0+k3 %確定卷積和f非零樣值的序號向量subplot(2,2,1)stem(k1,f1)%在子圖1繪序列f1(k)時域波形圖title(f1(k)xlabel(k)ylabel(f1(k)subplot(2,2,2)stem(k2,f2)%在圖2繪序列f2(k)時波形圖title(f1(k)xlabel(k)ylabel(f2(k)subplot(2,2,3)stem(k,f);%在子圖3繪序列f(k)的波形圖title(f(k)f1(k)與f2(k)的卷積和f(k)xlabel(k)ylabel(f(k)h=get(gca,position);h(3)=2.5*h(3);set(gca,position,h)%將第三個子圖的橫坐標范圍擴為原來的2.5倍例2-5求卷積和， f1=1 2 1;k1=-1 0 1;f2=ones(1,5);k2=-2:2;f, k=dconv(f1,f2,k1,k2) 由運行結果知，f的長度等于f1和f2長度之和減一， f的起點是f1和f2的起點之和，f的終點是f1和f2的終點之和。3 離散系統的單位響應MATLAB提供畫系統單位響應函數impz，調用格式是impz(b,a) 式中b和a是表示離散系統的行向量；impz(b,a,n) 式中b和a是表示離散系統的行向量，時間范圍是0n；impz(b,a,n1,n2) 時間范圍是n1n2 ；y=impz(b,a,n1,n2) 由y給出數值序列；例2-6已知 求單位響應。a=1,-1,0.9;b=1;impz(b,a)impz(b,a,60)impz(b,a,-10:40)4 離散系統的零狀態響應MATLAB提供求離散系統零狀態響應數值解函數filter，調用格式為filter(b,a,x),式中b和a是表示離散系統的向量，x是輸入序列非零樣值點行向量，輸出向量序號同x一樣。例2-7 已知 ， 求零狀態響應， 范圍020。a=1 -0.25 0.5;b=1 1;t=0:20;x=(1/2).t;y=filter(b,a,x)subplot(2,1,1)stem(t,x)title(輸入序列)subplot(2,1,2)stem(t,y)title(響應序列)三、上機實驗內容1驗證實驗原理中程序2已知，畫單位響應波形。3已知，輸入，畫輸出波形，范圍015。實驗三 連續時間LTI系統的時域分析一、實驗目的1學會用MATLAB求解連續系統的零狀態響應；2. 學會用MATLAB求解沖激響應及階躍響應； 3學會用MATLAB實現連續信號卷積的方法；二、實驗原理1連續時間系統零狀態響應的數值計算 我們知道，LTI連續系統可用如下所示的線性常系數微分方程來描述， 在MATLAB中，控制系統工具箱提供了一個用于求解零初始條件微分方程數值解的函數lsim。其調用格式y=lsim(sys,f,t)式中，t表示計算系統響應的抽樣點向量，f是系統輸入信號向量，sys是LTI系統模型，用來表示微分方程，差分方程或狀態方程。其調用格式sys=tf(b,a)式中，b和a分別是微分方程的右端和左端系數向量。例如，對于以下方程: 可用 獲得其LTI模型。注意，如果微分方程的左端或右端表達式中有缺項，則其向量a或b中的對應元素應為零，不能省略不寫，否則出錯。例3-1 已知某LTI系統的微分方程為 y(t)+ 2y(t)+100y(t)=f(t) 其中，求系統的輸出y(t).解：顯然，這是一個求系統零狀態響應的問題。其MATLAB計算程序如下： ts=0;te=5;dt=0.01; sys=tf(1,1,2,100); t=ts:dt:te; f=10*sin(2*pi*t); y=lsim(sys,f,t); plot(t,y); xlabel(Time(sec); ylabel(y(t);2連續時間系統沖激響應和階躍響應的求解 在MATLAB中，對于連續LTI系統的沖激響應和階躍響應，可分別用控制系統工具箱提供的函數impluse和step來求解。其調用格式為 y=impluse(sys,t) y=step(sys,t)式中，t表示計算系統響應的抽樣點向量，sys是LTI系統模型。例3-2已知某LTI系統的微分方程為 y(t)+ 2y(t)+100y(t)=10f(t)求系統的沖激響應和階躍響應的波形.解：ts=0;te=5;dt=0.01; sys=tf(10,1,2,100);t=ts:dt:te;h=impulse(sys,t);figure;plot(t,h); xlabel(Time(sec); ylabel(h(t);g=step(sys,t);figure;plot(t,g); xlabel(Time(sec);ylabel(g(t);3. 用MATLAB實現連續時間信號的卷積信號的卷積運算有符號算法和數值算法，此處采用數值計算法，需調用MATLAB的conv( )函數近似計算信號的卷積積分。連續信號的卷積積分定義是 如果對連續信號和進行等時間間隔均勻抽樣，則和分別變為離散時間信號和。其中，為整數。當足夠小時，和既為連續時間信號和。因此連續時間信號卷積積分可表示為 采用數值計算時，只求當時卷積積分的值，其中，n為整數，既 其中，實際就是離散序列和的卷積和。當足夠小時，序列就是連續信號的數值近似，既 上式表明，連續信號和的卷積，可用各自抽樣后的離散時間序列的卷積再乘以抽樣間隔。抽樣間隔越小，誤差越小。例3-3用數值計算法求與的卷積積分。解：因為是一個持續時間無限長的信號，而計算機數值計算不可能計算真正的無限長信號，所以在進行的抽樣離散化時，所取的時間范圍讓衰減到足夠小就可以了，本例取。程序是dt=0.01; t=-1:dt:2.5;f1=Heaviside(t)-Heaviside(t-2);f2=exp(-3*t).*Heaviside(t);f=conv(f1,f2)*dt; n=length(f); tt=(0:n-1)*dt-2;subplot(221), plot(t,f1), grid on;axis(-1,2.5,-0.2,1.2); title(f1(t); xlabel(t)subplot(222), plot(t,f2), grid on;axis(-1,2.5,-0.2,1.2); title(f2(t); xlabel(t)subplot(212), plot(tt,f), grid on;title(f(t)=f1(t)*f2(t); xlabel(t)由于和的時間范圍都是從t=-1開始，所以卷積結果的時間范圍從 t=-2開始，增量還是取樣間隔，這就是語句tt=(0:n-1)*dt-2的由來。三、上機實驗內容1 驗證實驗原理中所述的相關程序2 已知描述系統的微分方程和激勵信號f(t)如下，試用解析法求系統的零狀態響應y(t)，并用MATLAB繪出系統零狀態響應的時域仿真波形，驗證結果是否相同y(t)+ 4y(t)+4y(t)=f(t)+3f(t) f(t)= exp(-t)u(t) 3已知描述系統的微分方程如下，試用MATLAB求系統在010秒范圍內沖激響應和階躍響應的數值解，并用繪出系統沖激響應和階躍響應的時域波形y(t)+3y(t)+2y(t)=f(t)y(t)+ 2y(t)+2y(t)=f(t) 4畫出信號卷積積分的波形，實驗四 傅里葉變換、系統的頻域分析一、 實驗目的1、學會用MATLAB實現連續時間信號傅里葉變換2、學會用MATLAB分析LTI系統的頻域特性3、學會用MATLAB分析LTI系統的輸出響應二、實驗原理1傅里葉變換的MATLAB求解 MATLAB的symbolic Math Toolbox 提供了直接求解傅里葉變換及逆變換的函數fourier()及ifourier()兩者的調用格式如下。Fourier 變換的調用格式F=fourier(f)：它是符號函數f的fourier變換默認返回是關于w的函數。F=fourier(f，v):它返回函數F是關于符號對象v的函數，而不是默認的w，即 Fourier逆變換的調用格式f=ifourier(F):它是符號函數F的fourier逆變換，默認的獨立變量為w，默認返回是關于x的函數。f=ifourier(f,u):它的返回函數f是u的函數，而不是默認的x.注意：在調用函數fourier()及ifourier()之前，要用syms命令對所用到的變量（如t,u,v,w）進行說明,即將這些變量說明成符號變量。例4-1 求的傅立葉變換解: 可用MATLAB解決上述問題：syms tFw=fourier(exp(-2*abs(t)例4-2 求的逆變換f(t)解： 可用MATLAB解決上述問題syms t wft=ifourier(1/(1+w2),t)2連續時間信號的頻譜圖例4-3 求調制信號的頻譜，式中解：MATLAB程序如下所示ft=sym(4*cos(2*pi*6*t)*(Heaviside(t+1/4)-Heaviside(t-1/4);Fw=simplify(fourier(ft)subplot(121)ezplot(ft,-0.5 0.5),grid onsubplot(122)ezplot(abs(Fw),-24*pi 24*pi),grid 用MATLAB符號算法求傅里葉變換有一定局限，當信號不能用解析式表達時，會提示出錯，這時用MATLAB的數值計算也可以求連續信號的傅里葉變換，計算原理是當足夠小時，近似計算可滿足要求。若信號是時限的，或當時間大于某個給定值時，信號已衰減的很厲害，可以近似地看成時限信號時，n的取值就是有限的，設為N，有 是頻率取樣點時間信號取樣間隔應小于奈奎斯特取樣時間間隔，若不是帶限信號可根據計算精度要求確定一個頻率 W0為信號的帶寬。例4-4 用數值計算法求信號的傅里葉變換解，信號頻譜是，第一個過零點是，一般將此頻率視為信號的帶寬，若將精度提高到該值的50倍，既W0=50，據此確定取樣間隔，R=0.02;t=-2:R:2;f=Heaviside(t+1)-Heaviside(t-1);W1=2*pi*5;N=500;k=0:N;W=k*W1/N;F=f*exp(-j*t*W)*R;F=real(F);W=-fliplr(W),W(2:501);F=fliplr(F),F(2:501);subplot(2,1,1);plot(t,f);xlabel(t);ylabel(f(t);title(f(t)=u(t+1)-u(t-1);subplot(2,1,2);plot(W,F);xlabel(w);ylabel(F(w);title(f(t)的付氏變換F(w);3用MATLAB分析LTI系統的頻率特性 當系統的頻率響應H（jw）是jw的有理多項式時，有 MATLAB信號處理工具箱提供的freqs函數可直接計算系統的頻率響應的數值解。其調用格式如下H=freqs(b,a,w)其中，a和b分別是H(jw)的分母和分子多項式的系數向量，w為形如w1:p:w2的向量，定義系統頻率響應的頻率范圍，w1為頻率起始值，w2為頻率終止值,p為頻率取樣間隔。H返回w所定義的頻率點上，系統頻率響應的樣值。例如，運行如下命令，計算02pi頻率范圍內以間隔0.5取樣的系統頻率響應的樣值a=1 2 1;b=0 1;h=freqs(b,a,0:0.5:2*pi)例 4-5 三階歸一化的butterworth 低通濾波器的頻率響應為 試畫出該系統的幅度響應和相位響應。解 其MATLAB程序及響應的波形如下w=0:0.025:5;b=1;a=1,2,2,1;H=freqs(b,a,w);subplot(2,1,1);plot(w,abs(H);grid;xlabel(omega(rad/s);ylabel(|H(jomega)|);title(H(jw)的幅頻特性);subplot(2,1,2);plot(w,angle (H);grid;xlabel(omega(rad/s);ylabel(phi(omega);title(H(jw)的相頻特性);4用MATLAB分析LTI系統的輸出響應例 4-6已知一RC電路如圖所示 系統的輸入電壓為f(t)，輸出信號為電阻兩端的電壓y(t).當RC=0.04，f(t)=cos5t+cos100t, 試求該系統的響應y(t)-+-+f(t)y(t)RC解 由圖可知 ，該電路為一個微分電路，其頻率響應為 由此可求出余弦信號通過LTI系統的響應為 計算該系統響應的MATLAB程序及響應波形如下RC=0.04;t=linspace(-2,2,1024);w1=5;w2=100;H1=j*w1/(j*w1+1/RC);H2=j*w2/(j*w2+1/RC);f=cos(5*t)+cos(100*t);y=abs(H1)*cos(w1*t+angle(H1)+ abs(H2)*cos(w2*t+angle(H2);subplot(2,1,1);plot(t,f);ylabel(f(t);xlabel(Time(s);subplot(2,1,2);plot(t,y);ylabel(y(t);xlabel(Time(s);三、 上機實驗內容1.驗證實驗原理中所述的相關程序；2.試用MATLAB求單邊指數數信號的傅立葉變換，并畫出其波形；3.設，試用MATLAB畫出該系統的幅頻特性和相頻特性，并分析系統具有什么濾波特性。實驗五 信號抽樣與恢復一、實驗目的學會用MATLAB實現連續信號的采樣和重建二、實驗原理 1抽樣定理若是帶限信號，帶寬為, 經采樣后的頻譜就是將的頻譜 在頻率軸上以采樣頻率為間隔進行周期延拓。因此，當時，不會發生頻率混疊；而當 時將發生頻率混疊。2信號重建經采樣后得到信號經理想低通則可得到重建信號，即：=*其中：=所以：=*=* =上式表明，連續信號可以展開成抽樣函數的無窮級數。利用MATLAB中的來表示，有 ，所以可以得到在MATLAB中信號由重建的表達式如下：=我們選取信號=作為被采樣信號，當采樣頻率=2時，稱為臨界采樣。我們取理想低通的截止頻率=。下面程序實現對信號=的采樣及由該采樣信號恢復重建：例5-1 Sa(t)的臨界采樣及信號重構；wm=1; %信號帶寬wc=wm; %濾波器截止頻率Ts=pi/wm; %采樣間隔ws=2*pi/Ts; %采樣角頻率 n=-100:100; %時域采樣電數nTs=n*Ts %時域采樣點f=sinc(nTs/pi);Dt=0.005;t=-15:Dt:15;fa=f*Ts*wc/pi*sinc(wc/pi)*(ones(length(nTs),1)*t-nTs*ones(1,length(t); %信號重構t1=-15:0.5:15;f1=sinc(t1/pi);subplot(211);stem(t1,f1);xlabel(kTs);ylabel(f(kTs);title(sa(t)=sinc(t/pi)的臨界采樣信號);subplot(212);plot(t,fa)xlabel(t);ylabel(fa(t);title(由sa(t)=sinc(t/pi)的臨界采樣信號重構sa(t);grid;例5-2 Sa(t)的過采樣及信號重構和絕對誤差分析程序和例4-1類似，將采樣間隔改成Ts=0.7*pi/wm , 濾波器截止頻率該成wc=1.1*wm ，添加一個誤差函數wm=1;wc=1.1*wm; Ts=0.7*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Tsf=sinc(nTs/pi);Dt=0.005;t=-15:Dt:15;fa=f*Ts*wc/pi*sinc(wc/pi)*(ones(length(nTs),1)*t-nTs*ones(1,length(t);error=abs(fa-sinc(t/pi); %重構信號與原信號誤差t1=-15:0.5:15;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel(kTs);ylabel(f(kTs);title(sa(t)=sinc(t/pi)的采樣信號);subplot(312);plot(t,fa)xlabel(t);ylabel(fa(t);title(由sa(t)=sinc(t/pi)的過采樣信號重構sa(t);grid;subplot(313);plot(t,error);xlabel(t);ylabel(error(t);title(過采樣信號與原信號的誤差error(t);例5-3 Sa(t)的欠采樣及信號重構和絕對誤差分析程序和例4-2類似，將采樣間隔改成Ts=1.5*pi/wm , 濾波器截止頻率該成wc=wm=1三、上機實驗內容1驗證實驗原理中所述的相關程序；2設f(t)=0.5*(1+cost)*(u(t+pi)-u(t-pi) ，由于不是嚴格的頻帶有限信號，但其頻譜大部分集中在0，2之間，帶寬wm可根據一定的精度要求做一些近似。試根據以下兩種情況用 MATLAB實現由f(t)的抽樣信號fs(t)重建f(t) 并求兩者誤差，分析兩種情況下的結果。(1) wm=2 , wc=1.2wm , Ts=1;(2) wm=2 , wc=2 , Ts=2.5實驗六 信號與系統復頻域分析一、實驗目的1.學會用MATLAB進行部分分式展開；2.學會用MATLAB分析LTI系統的特性；3.學會用MATLAB進行Laplace正、反變換。4.學會用MATLAB畫離散系統零極點圖；5.學會用MATLAB分析離散系統的頻率特性；二、實驗原理及內容1用MATLAB進行部分分式展開用MATLAB函數residue可以得到復雜有理分式F(s)的部分分式展開式，其調用格式為 其中，num,den分別為F(s)的分子和分母多項式的系數向量，r為部分分式的系數，p為極點，k為F(s)中整式部分的系數，若F(s)為有理真分式，則k為零。例6-1 用部分分式展開法求F(s)的反變換 解:其MATLAB程序為format rat;num=1,2;den=1,4,3,0;r,p=residue(num,den)程序中format rat是將結果數據以分數形式顯示F(s)可展開為 所以，F(s)的反變換為 2用MATLAB分析LTI系統的特性系統函數H（s）通常是一個有理分式，其分子和分母均為多項式。計算H（s）的零極點可以應用MATLAB中的roots函數，求出分子和分母多項式的根，然后用plot命令畫圖。在MATLAB中還有一種更簡便的方法畫系統函數H（s）的零極點分布圖，即用pzmap函數畫圖。其調用格式為 pzmap(sys)sys表示LTI系統的模型，要借助tf函數獲得，其調用格式為sys=tf(b,a)式中，b和a分別為系統函數H（s）的分子和分母多項式的系數向量。如果已知系統函數H（s），求系統的單位沖激響應h(t)和頻率響應可以用以前介紹過的impulse和freqs函數。例6-2 已知系統函數為 試畫出其零極點分布圖，求系統的單位沖激響應h(t)和頻率響應，并判斷系統是否穩定。解：其MATLAB程序如下：num=1;den=1,2,2,1;sys=tf(num,den);figure(1);pzmap(sys);t=0:0.02:10;h=impulse(num,den,t);figure(2);plot(t,h)title(Impulse Response)H,w=freqs(num,den);figure(3);plot(w,abs(H)xlabel(omega)title(Magnitude Response) 3用MATLAB進行Laplace正、反變換MATLAB的符號數學工具箱提供了計算Laplace正、反變換的函數Laplace和ilaplace,其調用格式為 上述兩式右端的f和F分別為時域表示式和s域表示式的符號表示，可以應用函數sym實現，其調用格式為 S=sym(A)式中，A為待分析表示式的字符串，S為符號數字或變量。例6-3 試分別用Laplace和ilaplace函數求（1）的Laplace變換；（2）的Laplace反變換。解：（1）其程序為f=sym(exp(-t)*sin(a*t);F=laplace(f)或syms a tF=laplace(exp(-t)*sin(a*t)（2）其程序為F=sym(s2/(s2+1);ft=ilaplace(F)或syms sft= ilaplace(s2/(s2+1)4離散系統零極點圖離散系統可以用下述差分方程描述：Z變換后可得系統函數：用MATLAB提供的root函數可分別求零點和極點，調用格式是p=a0,a1an,q=b0,b1bm,0,00, 補0使二者維數一樣。畫零極點圖的方法有多種，可以用MATLAB函數z,p,k=tf2zp(b,a)和zplane(q,p)，也可用plot命令自編一函數ljdt.m,畫圖時調用。function ljdt(A,B)% The function to draw the pole-zero diagram for discrete systemp=roots(A)； %求系統極點q=roots(B)； %求系統零點p=p； %將極點列向量轉置為行向量q=q； %將零點列向量轉置為行向量x=max(abs(p q 1)；%確定縱坐標范圍x=x+0.1；y=x；%確定橫坐標范圍clfhold onaxis(-x x -y y)%確定坐標軸顯示范圍w=0:pi/300:2*pi；t=exp(i*w)；plot(t)%畫單位園axis(square)plot(-x x,0 0)%畫橫坐標軸plot(0 0,-y y)%畫縱坐標軸text(0.1,x,jImz)text(y,1/10,Rez)plot(real(p),imag(p),x)%畫極點plot(re