MATLAB在信號與系統(tǒng)中的應(yīng)用

MATLAB應(yīng)用篇MATLAB語言應(yīng)用篇

主講教師：任景英電氣信息工程學(xué)院電信教研室

本篇介紹定位：語言應(yīng)用必修課目標(biāo)：以信號與系統(tǒng)為基礎(chǔ)、直接利用教材上的結(jié)論，用MATLAB編程來解決實(shí)際問題特點(diǎn)：在掌握前面MATLAB基礎(chǔ)的情況下，應(yīng)用MATLAB語言解決電子信息領(lǐng)域中的問題。幫助更好的理解所學(xué)課程的原理

學(xué)習(xí)內(nèi)容

一、MATLAB在連續(xù)信號和系統(tǒng)中的應(yīng)用

二、MATLAB在傅立葉分析中的應(yīng)用

三、MATLAB在離散信號和系統(tǒng)中的應(yīng)用

本節(jié)討論用MATLAB表示和分析連續(xù)信號和線性時(shí)不變(LTI)連續(xù)系統(tǒng)的問題。嚴(yán)格說來，只有用符號推理的方法才能分析連續(xù)系統(tǒng)。用數(shù)值方法是不能表示連續(xù)信號的，因?yàn)樗o出的是各個(gè)樣本點(diǎn)的數(shù)據(jù)。只有當(dāng)樣本點(diǎn)取得很密時(shí)才可看成連續(xù)信號。所謂密，是相對于信號變化的快慢而言。以下均假定相對于采樣點(diǎn)密度而言，信號變化足夠慢。第一節(jié)MATLAB在連續(xù)信號和系統(tǒng)中的應(yīng)用

例3.1連續(xù)信號的MATLAB描述

列出單位沖激函數(shù)、單位階躍、復(fù)指數(shù)函數(shù)等連續(xù)信號的MATLAB表達(dá)式。

解:建模單位沖擊函數(shù)無法直接用MATLAB描述,可以把它看作是寬度為△(程序中用dt表示),幅度為1/△的矩形脈沖,即

單位階躍函數(shù):在t=t1處躍升的階躍可寫為u(t-t1),定義為

復(fù)指數(shù)函數(shù):

若w=0它是實(shí)指數(shù)函數(shù),如u=0則為虛指數(shù)函數(shù),其實(shí)部為余弦函數(shù),虛部為正弦函數(shù).本例u=-0.5,w=10%MATLAB實(shí)現(xiàn)的程序如下clear,t0=0;tf=5;dt=0.05;t1=1;t=[t0:dt:tf];%(1)單位脈沖信號%在t1(t0≤t1≤tf)處有一持續(xù)時(shí)間為dt,面積為1的脈沖信號,其余時(shí)間均為零。t=[t0:dt:tf];st=length(t);n1=floor((t1-t0)/dt); %求t1對應(yīng)的樣本序號x1=zeros(1,st); %把全部信號先初始化為零x1(n1)=1/dt; %給出t1處的脈沖信號subplot(2,2,1),stairs(t,x1),gridon %繪圖，注意為何用stairs命令%(2)單位階躍信號，%信號從t0到tf，在t1(t0≤t1≤tf)前為0，到t1處有一躍變，以后為1.% 程序前幾句即求t,st,n1的語句與上同，只把x1處改為x2x2=[zeros(1,n1-1),ones(1,st-n1+1)];%產(chǎn)生階躍信號subplot(2,2,3),stairs(t,x2),grid %繪圖axis([0,5,0,1.1])%(3)復(fù)數(shù)指數(shù)信號u=-0.5;w=10;x3=exp((u+j*w)*t);subplot(2,2,2),plot(t,real(x3)),grid %繪圖,subplot(2,2,4),plot(t,imag(x3)),grid %繪圖,實(shí)驗(yàn)結(jié)果如上圖線性微分方程的解(residue)線性常微分方程的解可用拉普拉斯算子，表示H(s)=B(s)/A(s)其中B(s)和A(s)都是s的多項(xiàng)式，分母多項(xiàng)式的次數(shù)n通常高于分子多項(xiàng)式的次數(shù)m。在時(shí)間域的解y(t)是Y(s)的拉普拉斯反變換。求反變換的重要方法之一是部分分式法，即將上述多項(xiàng)式分解為多個(gè)s的一次分式之和。留數(shù)函數(shù)residue可以完成這一任務(wù)。步驟:(1)用[r,p,k]=residue(b,a)求出Y(s)的極點(diǎn)數(shù)組p和留數(shù)數(shù)組r(設(shè)分母比分子階數(shù)高，故k=[])，因而Y(s)可表為Y(s)=r(1)/(s-p(1))+r(2)/(s-p(2))+r(3)/(s-p(3))+...(2)求它的反變換。得y(t)=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t)+r(3)*exp(p(3)*t)+r(4)*exp(p(4)*t)+...例3.2n階LTI系統(tǒng)的沖激響應(yīng)n階微分方程，寫成系統(tǒng)函數(shù)為：沖擊響應(yīng)就是H(s)的拉普拉斯反變換，可以把H(s)展開為極點(diǎn)留數(shù)式。

其反變換為a=input('多項(xiàng)式分母系數(shù)向量a=');b=input('多項(xiàng)式分子系數(shù)向量b=');[r,p]=residue(b,a), %求留數(shù)disp('解析式h(t)=Σr(i)*exp(p(i)*t)')disp('給出時(shí)間數(shù)組t=[0:dt:tf]')dt=input('dt='); tf=input('tf='); t=0:dt:tf;%給定時(shí)間數(shù)組h=zeros(1,length(t)); %h的初始化fori=1:length(a)-1 %根數(shù)為a的長度減1h=h+r(i)*exp(p(i)*t); %疊加各根分量endplot(t,h),gridset(gcf,'color','w')%設(shè)置圖形背景色為白色例3.3卷積的計(jì)算根據(jù)卷積公式：因此編程的過程為：（1）寫出h(t)的MATLAB表達(dá)式；（2）寫出u(t)的MATLAB表達(dá)式；（3）利用MATLAB的卷積語句y=conv(u，h)求解（4）畫曲線plot(t,y)。（1）調(diào)用conv函數(shù)的簡單程序u=input('輸入u數(shù)祖u='); %輸入u,h序列h=input('輸入h數(shù)祖h=');dt=input('輸入時(shí)間間隔dt=');y=conv(u,h); %調(diào)用conv函數(shù)%注意下面不知t的長度而要用它繪圖時(shí)的t數(shù)組設(shè)定方法plot(dt*([1:length(y)]-1),y),grid例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)的卷積求得。

cleara=input('多項(xiàng)式分母系數(shù)向量a=');b=input(‘多項(xiàng)式分子系數(shù)向量b=');t=input('輸入時(shí)間序列t=[0:dt:tf]');u=input('輸入序列u=');[r,p]=residue(b,a); %用極點(diǎn)留數(shù)法求沖擊函數(shù)h=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t); %沖擊函數(shù)h(t)subplot(2,1,1),plot(t,h);gridx=conv(u,h); %求u和h的卷積subplot(2,1,2),plot(t,x(1:length(t)));grid程序運(yùn)行結(jié)果執(zhí)行這個(gè)程序取a=[1，2，8]，b=1，t=[0:0.1:5]u=3*t+cos(0.1*t),所得的結(jié)果見圖6.5.圖6.5(a)顯示沖擊響應(yīng)h。。圖6.5沖擊響應(yīng)和卷積法求輸出用MATLAB的函數(shù)求解零狀態(tài)響應(yīng)1、連續(xù)時(shí)間系統(tǒng)零狀態(tài)響應(yīng)的求解LTI連續(xù)系統(tǒng)可用線性常系數(shù)微分方程來描述，系統(tǒng)的零狀態(tài)響應(yīng)可通過求解初始狀態(tài)為零的微分方程得到。在MATLAB的工具箱中，提供了一個(gè)用于求解零初始條件微分方程數(shù)值解的函數(shù)lsim，其調(diào)用格式為：y=lsim(sys,f,t)式中，t表示計(jì)算系統(tǒng)響應(yīng)的抽樣點(diǎn)向量，f是系統(tǒng)輸入信號向量。sys是LTI系統(tǒng)模型，用來表示微分方程、差分方程或狀態(tài)方程。其調(diào)用格式為：sys=tf(b,a)式中，b和a分別為微分方程的右端和左端系數(shù)向量。例：已知某LTI系統(tǒng)的微分方程為其中,,求系統(tǒng)的響應(yīng)。解：顯然，這是一個(gè)求系統(tǒng)零狀態(tài)響應(yīng)的問題。其MATLAB計(jì)算程序及響應(yīng)波形如下：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)');用MATLAB函數(shù)求沖激響應(yīng)和階躍響應(yīng)在MATLAB中，對于連續(xù)LTI系統(tǒng)的沖激響應(yīng)和階躍響應(yīng)，可分別用控制系統(tǒng)工具箱提供的函數(shù)impulse和step來求解。其調(diào)用格式為:y=impulse(sys,t)y=step(sys,t)式中，t表示計(jì)算系統(tǒng)響應(yīng)的抽樣點(diǎn)向量，sys是LTI系統(tǒng)模型。例、已知某LTI系統(tǒng)的微分方程為初始狀態(tài)為零，求系統(tǒng)的沖激響應(yīng)和階躍響應(yīng)。解：其MATLAB計(jì)算程序及響應(yīng)波形如下：ts=0;te=5;dt=0.01;sys=tf([10],[1,2,100]);t=ts:dt:te;h=impulse(sys,t);subplot(2,1,1),plot(t,h);g=step(sys,t);subplot(2,1,2),plot(t,g);第二節(jié)MATLAB在傅里葉分析中的應(yīng)用例3.4方波分解為多次正弦波之和圖示的周期性方波，其傅里葉級數(shù)為分別計(jì)算直到9次諧波，并做圖。圖6.7-1輸入周期性方波--xt=0:0.01:2*pi; %設(shè)定一個(gè)時(shí)間數(shù)組，有629個(gè)點(diǎn)y=sin(t);plot(t,y),pause%頻率為w=1(f=1/2π)的正弦基波fork=3:2:19y=y+sin(k*t)/k;%將各次諧波疊加直到第19次諧波plot(t,y),pauseend例3.5(1)周期信號的頻譜

幅度為A，寬度為，重復(fù)周期為T的矩形脈沖周期信號f(t)如下圖所示：當(dāng)A=1,=0.1s,T=0.5s時(shí)畫出其頻譜圖。建模：周期方波的傅立葉系數(shù)為：TtA-/20/2f(t)……建模：周期方波的傅立葉系數(shù)為：其中該頻譜的包絡(luò)線按照抽樣函數(shù)Sa(x)規(guī)律變化，第一個(gè)過零點(diǎn)為,即%matlab在周期信號頻譜分析中的應(yīng)用a=1;tao=0.1;t=0.5;%脈沖幅度、寬度和信號周期n0=t/tao;n=0:2*n0;%諧波次數(shù)（有效頻帶寬度的2倍）fn_p=a*tao/t*(sin(n*pi*tao/t+eps*(n==0)))./(n*pi*tao/t+eps*(n==0));fn_pabs=abs(fn_p);%取模fn_pang=angle(fn_p);%取相位fn_mabs=fliplr(fn_pabs(2:(2*n0)+1));%幅度譜偶對稱fn_mang=-fliplr(fn_pang(2:(2*n0)+1));%相位譜奇對稱fnabs=[fn_mabsfn_pabs];%雙邊幅度譜fnang=[fn_mangfn_pang];%雙邊相位譜subplot(2,1,1),stem((-2*n0:2*n0),fnabs);%輸出雙邊振幅頻譜text(4,0.11,'amplitudespectrum');subplot(2,1,2),stem((-2*n0:2*n0),fnang);%輸出相位頻譜text(-2,2,'phasespectrum');連續(xù)周期信號的頻譜分析例、利用MATLAB畫出下圖所示的周期三角波信號的頻譜圖。經(jīng)計(jì)算，該周期三角波信號的傅里葉級數(shù)系數(shù)為：畫出該信號頻譜的MATLAB源程序如下：N=10;n1=-N:-1;C1=-4*j*sin(n1*pi/2)/pi^2./n1.^2;C0=0;n2=1:N;C2=-4*j*sin(n2*pi/2)/pi^2./n2.^2;Cn=[C1C0C2];n=-N:N;subplot(2,1,1);stem(n,abs(Cn));subplot(2,1,2);stem(n,angle(Cn));(2)連續(xù)非周期信號傅里葉變換

MATLAB的工具箱提供了能直接求解傅里葉變換及逆變換的函數(shù)fourier()及ifourier()。兩者的調(diào)用格式如下:Fourier變換的調(diào)用格式為：

F=fourier(f)：它是符號函數(shù)f的Fourier變換，默認(rèn)返回是關(guān)于的函數(shù)。Fourier逆變換的調(diào)用格式為：

f=ifourier(F)：它是符號函數(shù)F的ifourier變換，默認(rèn)的獨(dú)立變量為，默認(rèn)返回是關(guān)于x的函數(shù)。注意：在調(diào)用函數(shù)fourier()及ifourier()之前，要用syms命令對所用到的變量（如t、）等進(jìn)行說明，即將這些變量說明成符號變量。例：求的傅里葉變換解：輸入下列兩行語句：symstfourier(exp(-2*abs(t)))可得ans=4/(4+w^2)例：求的傅里葉逆變換f(t).symstwifourier(1/(1+w^2),t)可得ans=1/2*exp(-t)*Heaviside(t)+1/2*exp(t)*Heaviside(-t)其中，Heaviside(t)函數(shù)即為單位階躍函數(shù)

本節(jié)討論用MATLAB表示離散信號(序列)和線性時(shí)不變(LTI)離散系統(tǒng)的問題。由于MATLAB數(shù)值計(jì)算的特點(diǎn)，用它來分析離散的信號與系統(tǒng)是很方便的。在MATLAB中，可以用一個(gè)向量來表示一個(gè)有限長度的序列。然而這樣一個(gè)向量并沒有包含采樣位置的信息。因此，完全地表示x(n)要用x和n兩個(gè)向量，例如序列x(n)={2,1,-1,3,1,4,3,7}(下面的箭頭為第0個(gè)采樣點(diǎn))在MATLAB中表示為:n=[-3,-2,-1,0,1,2,3,4];x=[2,1,-1,3,1,4,3,7];當(dāng)不需要采樣位置信息或這個(gè)信息是多余的時(shí)候(例如該序列從n=0開始)，可以只用x向量來表示。

第三節(jié)MATLAB在離散信號與系統(tǒng)中的應(yīng)用↑例3.6離散信號的MATLAB表述

編寫MATLAB程序來產(chǎn)生下列基本脈沖序列:(1)單位脈沖序列，起點(diǎn)n0，終點(diǎn)nf，在ns處有一單位脈沖(n0≤ns≤nf)。遲延的單位脈沖序列(2)單位階躍序列，起點(diǎn)n0，終點(diǎn)nf，在ns前為0，在ns后為1(n0≤ns≤nf)。(3)復(fù)指數(shù)序列MATLAB實(shí)現(xiàn)的程序如下clear,n0=0;nf=10;ns=3;n1=n0:nf;x1=[zeros(1,ns-n0),1,zeros(1,nf-ns)];%單位脈沖序列的產(chǎn)生%n1=n0:nf;x1=[(n1-ns)==0]; %顯然,用邏輯式是比較高明的方法n2=n0:nf;x2=[zeros(1,ns-n0),ones(1,nf-ns+1)];%單位階躍序列的產(chǎn)生%n1=n0:nf;x2=[(n1-ns)>=0]; %用邏輯式產(chǎn)生單位階躍序列n3=n0:nf;x3=exp((-0.2+0.5j)*n3); %復(fù)數(shù)指數(shù)序列subplot(2,2,1),stem(n1,x1);title('單位脈沖序列')subplot(2,2,3),stem(n2,x2);title('單位階躍序列')subplot(2,2,2),stem(n3,real(x3));%畫x3實(shí)部序列l(wèi)ine([0,10],[0,0]) %畫橫坐標(biāo)title('復(fù)指數(shù)序列'),ylabel('實(shí)部')subplot(2,2,4),stem(n3,imag(x3));%畫x3虛部line([0,10],[0,0]),ylabel('虛部')圖3.6幾種基本的離散信號例3.7離散時(shí)間系統(tǒng)的零狀態(tài)響應(yīng)LTI離散系統(tǒng)一般用差分方程來描述：已知差分方程的n個(gè)初始狀態(tài)及輸入f(k)，就可以通過編程由下式迭代算出系統(tǒng)的輸出：在初始狀態(tài)為零的情況下，MATLAB的信號處理工具箱提供了一個(gè)函數(shù)，可求系統(tǒng)的零狀態(tài)響應(yīng)：y=filter(b,a,f)式中，a、b分別是差分方程左、右端的系數(shù)向量，f表示輸入序列，y表示輸出序列。注意輸出序列的長度和輸入序列的長度相同。例：受噪聲干擾的信號為f(k)=s(k)+d(k)，其中是原始信號，d(k)是噪聲。已知M點(diǎn)滑動(dòng)平均系統(tǒng)的輸入輸出關(guān)系為：

試?yán)肕ATLAB編程實(shí)現(xiàn)用M點(diǎn)滑動(dòng)平均系統(tǒng)對受噪聲干擾的信號去噪。解：建模系統(tǒng)的輸入信號f(k)含有有用信號s(k)和噪聲信號d(k)。噪聲信號d(k)可以用MATLAB中的rand函數(shù)產(chǎn)生，將其疊加在有用信號s(k)上，即得到受噪聲干擾的輸入信號f(k)。下面的程序?qū)崿F(xiàn)了對信號去噪，取M=5。MATLAB實(shí)現(xiàn)的程序如下%通過滑動(dòng)平均系統(tǒng)對信號進(jìn)行去噪r=51;%設(shè)定輸入信號的長度d=rand(1,r)-0.5;

%在（-0.5，0.5）產(chǎn)生一個(gè)隨機(jī)信號k=0:r-1;s=2*k.*(0.9.^k);f=s+d;subplot(1,2,1);plot(k,d,'r-.',k,s,'b:',k,f,'k-');m=5;b=ones(1,m)/m;a=1;y=filter(b,a,f);subplot(1,2,2);plot(k,s,'b:',k,y,'r-');例3.8離散時(shí)間系統(tǒng)的沖激響應(yīng)和階躍響應(yīng)在MATLAB中，求解離散系統(tǒng)的單位沖激響應(yīng)，可應(yīng)用信號處理工具箱提供的函數(shù)impz，其調(diào)用格式為：h=impz(b,a,k)式中，a、b分別是差分方程左、右端的系數(shù)向量，k表示輸出序列的取值范圍（可?。?。h就是系統(tǒng)的單位沖激響應(yīng)對于MATLAB6.x及以上版本，在信號處理工具箱中還提供了求解離散系統(tǒng)單位階躍響應(yīng)的函數(shù)s