第七講MATLAB在電路中的應(yīng)用

利用MATLAB進(jìn)行信號與系統(tǒng)分析

信號的MATLAB表示信號的基本運(yùn)算信號的能量和功率利用MATLAB對信號進(jìn)行時域分析利用MATLAB對信號進(jìn)行頻域分析傅里葉(Fourier)變換基本信號的MATLAB表示

指數(shù)信號Aeat、指數(shù)序列ak、抽樣函數(shù)Sa(t)、正弦型信號、矩形脈沖信號、三角脈沖信號一、信號的MATLAB表示一、信號的MATLAB表示1、連續(xù)時間信號的MATLAB表示連續(xù)時間信號：時間變化連續(xù)如y=x(t)離散時間信號(序列)：時間離散如x(nT)=x(t)|t=nT例：用MATLAB繪制連續(xù)時間信號關(guān)于t的曲線，t的取值范圍為0-30s，并以0.1s遞增。%file65.mt=0:0.1:30;x=exp(-0.707*t).*sin(2*t/3);plot(t,x);gridon;xlabel('Time(sec)');ylabel('x(t)');函數(shù)名功能函數(shù)名功能sawtooth(t,width)產(chǎn)生鋸齒波或三角波信號(周期為2π，幅值－1~1)pulstran(t,d,’func’)產(chǎn)生沖激串square(t,duty)產(chǎn)生方波信號(周期為2π，幅值－1~1)rectpuls(t,width)產(chǎn)生非周期的方波信號sinc(t)產(chǎn)生抽樣函數(shù)波形tripuls(t,width,skew)產(chǎn)生非周期的三角波信號chirp(t,f0,t1,f1,’method’,phi)產(chǎn)生調(diào)頻余弦掃頻信號diric(x,n)產(chǎn)生Dirichlet或周期sinc函數(shù)gauspuls(T,FC,BW,BWR)產(chǎn)生高斯正弦脈沖信號gmonopuls(t,fc)產(chǎn)生高斯單脈沖信號voc(x,fc,fs)電壓控制振蕩器2、工具箱中的信號產(chǎn)生函數(shù)

舉例（1）%產(chǎn)生周期為0.02的三角波file66.mt=0:0.0001:1;x1=sawtooth(2*pi*50*t,0.8);plot(t,x1);axis([0,0.2,-1,1])%產(chǎn)生頻率為50HZ，占空比為50%的方波file67.mt=0:0.0001:0.08;x1=square(2*pi*50*t,50);plot(t,x1);axis([0,0.08,-1,1.5])舉例（2）%產(chǎn)生線性調(diào)頻信號：file68.mt=0:0.001:2;y=chirp(t,0,1,150);figure(1);plot(t,y);axis([0,0.5,0,1])%產(chǎn)生非周期三角波信號：file69.mt=0:0.01:1;y1=tripuls(t);y2=tripuls(t,0.6);subplot(2,1,1);plot(t,y1);subplot(2,1,2);plot(t,y2);3、離散時間信號的表示

在MATLAB中，離散時間信號x(n)的表示：需要用一個向量x表示序列幅值，用另一個等長的向量來定位時間變量n來完整地表示一個序列。例：繪制離散時間信號的棒狀圖。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=0,x(4)=-1。%離散時間信號的表示file70.mn=-3:5;%定位時間變量x=[0,0,-1,1,2,1,-1,0,0];stem(n,x);grid;%繪制火柴桿圖line([-3,5],[0,0]);%畫x軸線xlabel('n');ylabel('x[n]')4、幾種常用離散時間信號的表示

1.單位脈沖序列x=zeros(1,N);x(1,n0)=1;2.單位階躍序列n=[ns:nf];x=[(n-n0)>=0]3.實(shí)指數(shù)序列n=[ns:nf];x=a.^n4.復(fù)指數(shù)序列n=[ns:nf];x=exp((sigema+jw)*n)5.正(余)弦序列n=[ns:nf];x=cos(w*n+sita)%decayingexponentialt=0:001:10;%file71.mA=1;a=-0.4;ft=A*exp(a*t);plot(t,ft)t=0:0.5:10;A=1;a=-0.4;ft=A*exp(a*t);stem(t,ft)舉例（1）一、信號的MATLAB表示%rectpulsfile72.mt=0:0.001:4;T=1;ft=rectpuls(t-2*T,T);plot(t,ft)axis([0,4,-0.5,1.5])一、信號的MATLAB表示舉例（2）%tripuls

file73.mt=-3:0.001:3;ft=tripuls(t,4,0.5);plot(t,ft)ft1=tripuls(t,4,1);plot(t,ft1)一、信號的MATLAB表示舉例（3）%unitimpulssequence

file75.mk=-50:50;delta=[zeros(1,50),1,zeros(1,50)];stem(k,delta)function[f,k]=impseq(k0,k1,k2)%產(chǎn)生

f[k]=delta(k-k0)；k1<=k<=k2k=[k1:k2];f=[(k-k0)==0];k0=0;k1=-50;k2=50;%file74.m[f,k]=impseq(k0,k1,k2);stem(k,f)一、信號的MATLAB表示舉例（4）%unitstepsequencek=-50:50;%file76.muk=[zeros(1,50),ones(1,51)];stem(k,uk)function[f,k]=stepseq(k0,k1,k2)%產(chǎn)生

f[k]=u(k-k0);k1<=k<=k2k=[k1:k2];f=[(k-k0)>=0];k0=0;k1=-50;k2=50;[f,k]=stepseq(k0,k1,k2);stem(k,f)function[f,k]=stepseq(k0,k1,k2)%產(chǎn)生

f[k]=u(k-k0);k1<=k<=k2k=[k1:k2];f=[(k-k0)>=0];k0=0;k1=-50;k2=50;%file77.m[f,k]=stepseq(k0,k1,k2);stem(k,f)一、信號的MATLAB表示舉例（5）求的實(shí)部序列，虛部序列，幅度序列，相位序列

fsxl(0.1,pi/3,-10,10)file78.mfunctionfsxl(delta,omig,n1,n2)e=delta+omig.*j;n=[n1:n2];x=exp(e*n);x_real=real(x);x_imag=imag(x);x_magnitude=abs(x);x_phase=(180/pi)*angle(x)subplot(2,2,1);stem(n,x_real);xlabel('n');ylabel('x_real');title('realpart');subplot(2,2,2);stem(n,x_imag);xlabel('n');ylabel('x_imag');title('imaginarypart');subplot(2,2,3);stem(n,x_magnitude);xlabel('n');ylabel('x_magnitude');title('amplitudepart');subplot(2,2,4);stem(n,x_phase);xlabel('n');ylabel('x_phase');title('phasepart');二、信號的基本運(yùn)算%信號的翻轉(zhuǎn)和時移file79.mt=-3:0.001:3;ft1=tripuls(2*t,4,0.5);subplot(2,1,1)plot(t,ft1)title('f(2t)')ft2=tripuls((2-2*t),4,0.5);subplot(2,1,2)plot(t,ft2)title('f(2-2t)')1.

信號的尺度變換、翻轉(zhuǎn)、時移（平移）

已知三角波f(t)，用MATLAB畫出的f(2t)和f(2-2t)波形

二、信號基本運(yùn)算2.

信號的相加與相乘相加用算術(shù)運(yùn)算符“+”實(shí)現(xiàn)相乘用數(shù)組運(yùn)算符“.*”實(shí)現(xiàn)

例:畫信號Aeatcos(w0t+f)的波形

t=0:0.001:8;%file80.mA=1;a=-0.4;w0=2*pi;phi=0;ft1=A*exp(a*t).*sin(w0*t+phi);plot(t,ft1)二、信號基本運(yùn)算3.

序列移位與周期延拓運(yùn)算序列移位：y(n)=x(n-m)MATLAB實(shí)現(xiàn)：y=x;ny=nx-m

file85.m序列周期延拓：y(n)=x((n))MMATLAB實(shí)現(xiàn)：ny=nxs:nxf；y=x(mod(ny,M)+1)file86.mn1=[-5:4];%序列x1(n)的起始及終止位置file81.mn1s=-5;n1f=4;x1=[231-13421-5-3];%序列x1(n)不同時間的幅度n2=[0:9];%序列x2(n)的起始及終止位置n2s=0;n2f=9;x2=[1111111111];%序列x2(n)不同時間的幅度ns=min(n1s,n2s);nf=max(n1f,n2f);%新信號的時間起始終止位置n=ns:nf;y1=zeros(1,length(n));%延拓序列初始化y2=zeros(1,length(n));y1(find((n>=n1s)&(n<=n1f)==1))=x1;%給延拓序列y1賦值x1y2(find((n>=n2s)&(n<=n2f)==1))=x2;%給延拓序列y2賦值x2ya=y1+y2;%逐點(diǎn)相加yp=y1.*y2;%逐點(diǎn)相乘subplot(4,1,1);stem(n,y1,'.');subplot(4,1,2);stem(n,y2,'.');subplot(4,1,3);stem(n,ya,'.');subplot(4,1,4);stem(n,yp,'.');%file82.m周期延拓N=24;M=8;m=3;n=0:N-1;x1=(0.8).^n;%產(chǎn)生指數(shù)序列x2=[(n>=0)&(n<=M)];%形成矩形序列Rm(n)x=x1.*x2;%截取操作形成新序列x(n)xm=zeros(1,N);fork=m+1:m+Mxm(k)=x(k-m);%產(chǎn)生移位序列x(n-3)Endxc=x(mod(n,M)+1);%產(chǎn)生x(n)的周期延拓x(n)8xcm=x(mod(n-m,M)+1);%產(chǎn)生x(n)的周期延拓x(n-3)8subplot(4,1,1);stem(n,x,'.');ylabel('x(n)');subplot(4,1,2);stem(n,xm,'.');ylabel('x(n-3)');subplot(4,1,3);stem(n,xc,'.');ylabel('x(n)_8');subplot(4,1,4);stem(n,xcm,'.');ylabel('x(n-3)_8')4.序列翻褶與序列累加運(yùn)算序列翻褶數(shù)學(xué)描述：y(n)=x(-n)。MATLAB可實(shí)現(xiàn):y=fliplr(x)序列累加的數(shù)學(xué)描述為：MATLAB實(shí)現(xiàn)：y=cumsum(x)5.兩序列的卷積運(yùn)算兩序列卷積運(yùn)算：MATLAB實(shí)現(xiàn)：y=conv(x1,x2)。序列x1(n)和x2(n)必須長度有限。二、信號基本運(yùn)算6.兩序列的相關(guān)運(yùn)算兩序列相關(guān)運(yùn)算：MATLAB實(shí)現(xiàn)：y=xcorr(x1,x2)。例：求序列的翻褶序列y(n)=x(-n)和y(n)=x(-n+3)及累加和。file83.m二、信號基本運(yùn)算n=0:10;x=3*exp(-0.2*n);y=fliplr(x);n1=-fliplr(n);n2=-fliplr(n-3);%在指定位置m=3處的時間序列翻褶s=cumsum(x);subplot(2,2,1);stem(n,x);xlabel('n');ylabel('x(n)');subplot(2,2,2);stem(n2,y);xlabel('n');ylabel('y(n)=x(-3+n)');subplot(2,2,3);stem(n1,y);xlabel('n');ylabel('y(n=)x(-n)');subplot(2,2,4);stem(n,s);xlabel('n');ylabel('s(n)');1.信號能量數(shù)字定義：MATLAB實(shí)現(xiàn):E=sum(x.*conj(x));

或E=sum(abs(x).^2);數(shù)字定義：2.信號功率MATLAB實(shí)現(xiàn):P=sum(x.*conj(x))/N;

或E=sum(abs(x).^2)/N;非周期三角波信號能量的MATLAB計算：file84.mdt=0.0001;t=0:dt:1;x=tripuls(t);E=sum(abs(x).^2*dt);余弦信號的平均功率MATLAB計算：file85.mdt=0.0001;t=0:dt:2*pi;x=cos(t);P=sum(abs(x).^2*dt)./(2*pi);一個信號處理就是將輸入信號變換成輸出信號的運(yùn)算過程，如圖所示。在此過程中，輸出的信號稱為系統(tǒng)對輸入信號作出的響應(yīng)，輸入信號稱為系統(tǒng)的激勵信號。當(dāng)輸入信號為連續(xù)時間信號x(t)時，該系統(tǒng)稱為連續(xù)時間處理系統(tǒng)，可描述為：y(t)=T[x(t)]當(dāng)輸入信號為離散時間信號x[n]時，該系統(tǒng)稱為離散時間處理系統(tǒng)，可描述為：y[n]=T[x[n]]當(dāng)一個系統(tǒng)具有可加性和齊次性，則該系統(tǒng)稱為線性系統(tǒng)。如果系統(tǒng)響應(yīng)與激勵加于系統(tǒng)的時刻無關(guān)，則該系統(tǒng)為時不變系統(tǒng)。1.系統(tǒng)的描述(1)．常系數(shù)線性微分/差分方程(2)．系統(tǒng)傳遞函數(shù)連續(xù)系統(tǒng)：離散系統(tǒng)：系統(tǒng)傳遞函數(shù)的MATLAB表示：注意：在MATLAB中，傳遞函數(shù)由分子、分母兩個多項(xiàng)式的系數(shù)來表示，系數(shù)為降冪排列。例如：可以表示為：num=[1,0.2,1];den=[1,0.5,0.7]如果有空項(xiàng)的話，需要用零來補(bǔ)齊。sys=tf(num,den)%獲得連續(xù)系統(tǒng)的傳遞函數(shù)形式sys=tf(num,den,0,1,’variable’,’z^-1’)%獲得離散系統(tǒng)的傳遞函數(shù)形式(3)．零－極點(diǎn)增益模型zplane連續(xù)系統(tǒng)：離散系統(tǒng)：(4)．極點(diǎn)留數(shù)模型離散系統(tǒng)：連續(xù)系統(tǒng)：在MATLAB中，增益系數(shù)、零點(diǎn)向量、極點(diǎn)向量用三個列向量表示：增益系數(shù)(Gain)：k零點(diǎn)向量(Zero)：z=[z1z2

…zn]’極點(diǎn)向量(Pole)：p=[p1p2

…pn]’sys=zpk(z,p,k)%獲得零-極點(diǎn)模型表達(dá)式(6)．狀態(tài)空間模型ss連續(xù)系統(tǒng)：離散系統(tǒng)：(5)．二次分式模型連續(xù)系統(tǒng)：離散系統(tǒng)：狀態(tài)向量：x輸出向量：y激勵向量(輸入向量)：u在MATLAB中，用矩陣A，B，C，D表示系統(tǒng)的狀態(tài)空間模型。系統(tǒng)模型的MATLAB表示系統(tǒng)模型表示名稱參數(shù)傳遞函數(shù)模型tfnum-傳遞函數(shù)的分子系數(shù)向量den-傳遞函數(shù)的分母系統(tǒng)向量零極點(diǎn)增益模型zpz-系統(tǒng)零點(diǎn)向量;p-系統(tǒng)極點(diǎn)向量;k-增益系數(shù)二次分式模型sossos-二次分?jǐn)?shù)的系數(shù)矩陣;g-比例系數(shù)狀態(tài)空間模型ss[A,B,C,D]-狀態(tài)空間模型參數(shù)部分分式模型residuezR-分子向量;P-分母系數(shù)向量;K-常數(shù)在MATLAB中，用sos、ss、tf、zp分別表示二次分式模型、狀態(tài)空間模型、傳遞函數(shù)模型和零-極點(diǎn)增益模型。其中sos表示二次分式，g為比例系數(shù)，sos為L×6的矩陣。2.系統(tǒng)模型的轉(zhuǎn)換函數(shù)(1)．ss2tf函數(shù)

格式：[num,den]=ss2tf(A,B,C,D,iu)功能：將指定輸入量iu的線性系統(tǒng)(A,B,C,D)轉(zhuǎn)換為傳遞函數(shù)模型[num,den]。(2)．zp2tf函數(shù)格式：[num,den]=zp2tf(z,p,k)功能：將給定系統(tǒng)的零－極點(diǎn)增益模型轉(zhuǎn)換為傳遞函數(shù)模型，z、p、k分別為零點(diǎn)列向量、極點(diǎn)列向量和增益系數(shù)。(3).tf2sos函數(shù)格式：[sos,g]=tf2sos(num,den)功能：將給定系統(tǒng)的傳遞函數(shù)模型[num,den]轉(zhuǎn)換為二次分式模型sos,g為增益系數(shù)。線性系統(tǒng)模型的轉(zhuǎn)換函數(shù)函數(shù)名功能說明函數(shù)名功能說明ss2tf狀態(tài)空間模型轉(zhuǎn)換為傳遞函數(shù)模型tf2ss傳遞函數(shù)模型轉(zhuǎn)換為狀態(tài)空間模型ss2zp狀態(tài)空間模型轉(zhuǎn)換為零－極點(diǎn)增益模型zp2ss零－極點(diǎn)增益模型轉(zhuǎn)換為狀態(tài)空間模型zp2tf零－極點(diǎn)增益模型轉(zhuǎn)換為傳遞函數(shù)模型tf2zp傳遞函數(shù)模型轉(zhuǎn)換為零－極點(diǎn)增益模型residue傳遞函數(shù)模型轉(zhuǎn)換為極點(diǎn)留數(shù)模型residue極點(diǎn)留數(shù)模型轉(zhuǎn)換為傳遞函數(shù)模型ss2sos狀態(tài)空間模型轉(zhuǎn)換為二次分式模型sos2ss二次分式模型轉(zhuǎn)換為狀態(tài)空間模型tf2sos傳遞函數(shù)模型轉(zhuǎn)換為二次分式模型sos2tf二次分式模型轉(zhuǎn)換為傳遞函數(shù)模型zp2sos零－極點(diǎn)增益模型轉(zhuǎn)換為二次分式模型sos2zp二次分式模型轉(zhuǎn)換為零－極點(diǎn)增益模型例：將系統(tǒng)轉(zhuǎn)換為狀態(tài)空間模型[A,B,C,D]在命令窗口輸入：>>z=[3];p=[1,2];k=2;>>[A,B,C,D]=zp2ss(z,p,k);

屏幕顯示為:A=3.0000-1.41421.41420B=10C=2.0000-4.2426D=0例：求離散時間系統(tǒng)的零、極點(diǎn)向量和增益系數(shù)。在命令窗口輸入：>>num=[2,3];den=[1,0.4,1];>>[num,den]=eqtflength(num,den);

%使長度相等>>[z,p,k]=tf2zp(num,den)屏幕顯示為z=

0-1.5000p=-0.2000+0.9798i

-0.2000-0.9798ik=23.系統(tǒng)互聯(lián)與系統(tǒng)結(jié)構(gòu)MATLAB實(shí)現(xiàn)函數(shù)series()

格式：[A,B,C,D]=series(A1,B1,C1,D1,A2,B2,C2,D2)或

[num,den]=series(num1,den1,num2,den2)將系統(tǒng)1、系統(tǒng)2級聯(lián)，可得到級聯(lián)連接的傳遞函數(shù)形式為：(1).系統(tǒng)的級聯(lián)(2).系統(tǒng)的并聯(lián)MATLAB實(shí)現(xiàn)函數(shù)parallel()格式：[A,B,C,D]=parallel(A1,B1,C1,D1,A2,B2,C2,D2)或

[num,den]=parallel(num1,den1,num2,den2)將系統(tǒng)1、系統(tǒng)2并聯(lián)，可得到并聯(lián)連接的傳遞函數(shù)形式為：(3).兩個系統(tǒng)的反饋連接函數(shù)feedback格式：[A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign)或

[num,den]=feedback(num1,den1,num2,den2,sign)將系統(tǒng)1和系統(tǒng)2進(jìn)行反饋連接，sign表示反饋方式(默認(rèn)值為-1)；sig=+1時表示正反饋；當(dāng)sig=-1時表示負(fù)反饋。例：求兩個單輸入單輸出子系統(tǒng)的級聯(lián)、并聯(lián)和反饋后系統(tǒng)的傳遞函數(shù)。file90.mMATLAB源程序?yàn)椋篺ile90.mnum1=1;den1=[1,1];%系統(tǒng)1num2=2;den2=[1,2];%系統(tǒng)2[nums,dens]=series(num1,den1,num2,den2)%實(shí)現(xiàn)兩個系統(tǒng)級聯(lián)[nump,denp]=parallel(num1,den1,num2,den2)%實(shí)現(xiàn)兩個系統(tǒng)并聯(lián)[numf,denf]=feedback(num1,den1,num2,den2)%實(shí)現(xiàn)兩個系統(tǒng)反饋程序運(yùn)行結(jié)果為：nums=002；dens=132nump=034；denp=132numf=012；denf=134因此，各系統(tǒng)的傳遞函數(shù)分別為：

而在實(shí)際應(yīng)用中，也可以把一個復(fù)雜的線性時不變系統(tǒng)分解為幾個簡單系統(tǒng)的組合結(jié)構(gòu)，即直接型結(jié)構(gòu)、級聯(lián)型結(jié)構(gòu)和并聯(lián)型結(jié)構(gòu)。MATLAB所提供的系統(tǒng)模型變換函數(shù)，實(shí)質(zhì)就是給出了這幾種系統(tǒng)結(jié)構(gòu)的互換關(guān)系。系統(tǒng)的傳遞函數(shù)對應(yīng)于系統(tǒng)的直接型結(jié)構(gòu)，二次分式(sos)模型對應(yīng)于級聯(lián)型結(jié)構(gòu)，系統(tǒng)傳遞函數(shù)的部分分式(residue或residuez)形式對應(yīng)于并聯(lián)型結(jié)構(gòu)。例：已知某數(shù)字濾波器的傳遞函數(shù)為：請求出其級聯(lián)型結(jié)構(gòu)和格型結(jié)構(gòu)。file91.m%file87.mclearall;b=[2,13/12,5/4,2/3];a=[1];%設(shè)定參數(shù)fprintf('級聯(lián)型結(jié)構(gòu)系數(shù)：');[sos,g]=tf2sos(b,a);%直接型到級聯(lián)型轉(zhuǎn)換fprintf('格型結(jié)構(gòu)系數(shù)(反射系數(shù))：');[K]=tf2latc(b,a);%直接型到格型轉(zhuǎn)換sos=1.00000.536001.0000001.00000.00570.62191.000000g=2格型結(jié)構(gòu)系數(shù)（反射系數(shù)）：K=0.25000.50000.3333根據(jù)級聯(lián)結(jié)構(gòu)系數(shù)，其傳遞函數(shù)可以改寫為：例：已經(jīng)某數(shù)字濾波器的傳遞函數(shù)為

請求出其級聯(lián)型結(jié)構(gòu)和并聯(lián)型結(jié)構(gòu)%file88.m

b=[1,-3,11,-27,18];a=[16,12,2,-4,-1];disp('級聯(lián)型結(jié)構(gòu)系數(shù)');[sos,g]=tf2sos(b,a);disp('并聯(lián)型結(jié)構(gòu)系數(shù)');[R,P,K]=residuez(b,a);級聯(lián)型結(jié)構(gòu)系數(shù)：sos=1.0000-3.00002.00001.0000-0.2500-0.12501.00000.00009.00001.00001.00000.5000g=0.0625并聯(lián)型結(jié)構(gòu)系數(shù)：R=-5.0250-1.0750i-5.0250+1.0750i0.925027.1875P=

-0.5000+0.5000i-0.5000-0.5000i0.5000-0.2500

K=

-18根據(jù)級聯(lián)結(jié)構(gòu)系數(shù)，其傳遞函數(shù)可以改寫為：根據(jù)并聯(lián)結(jié)構(gòu)系數(shù)，其傳遞函數(shù)可改寫為：模型轉(zhuǎn)換函數(shù)示意圖(1)．連續(xù)LTI系統(tǒng)的響應(yīng)(2)．離散LTI系統(tǒng)的響應(yīng)用MATLAB中的卷積函數(shù)conv()來實(shí)現(xiàn)。用MATLAB中的卷積函數(shù)conv()來實(shí)現(xiàn)。%file89.mdt=input('輸入時間間隔dt')x=ones(1,fix(10/dt));h=exp(-0.1*[0:fix(10/dt)]*dt);y=conv(x,h);t=dt*([1:length(y)]-1);plot(t,y);grid;例：某LT1系統(tǒng)的單位沖激響應(yīng)，輸入初始條件為零，求系統(tǒng)的響應(yīng)y(t)。%file90.mx=ones(1,10);n=[0:14];h=0.5.^n;y=conv(x,h);stem(y);xlabel('n');ylabel('y[n]');例：已知某LT1離散系統(tǒng)的單位沖激響應(yīng)為：求輸入信號序列的系統(tǒng)響應(yīng)。nx=-5:4;x=ones(1,10);nh=[0:14];h=0.5.^n;y=conv(x,h);n0=nx(1)+nh(1);%求卷積序列y起始時間位置N=length(nx)+length(nh)-2;%求卷積序列y序列長度ny=n0:n0+N;%求卷積序列y的時間向量subplot(2,2,1);stem(nx,x);title('x[n]');xlabel('n');xlabel('x[n]');subplot(2,2,2);stem(nh,h);title('h[n]');xlabel('n');xlabel('h[n]');subplot(2,2,3);stem(ny,y);title('x[n]*h[n]');xlabel('n');xlabel('y[n]');h=get(gca,'position');h(3)=2.5*h(3)set(gca,'position',h)改寫為：file91.m(3)．時域響應(yīng)函數(shù)(a)連續(xù)系統(tǒng)的單位沖激響應(yīng) 調(diào)用格式：[Y,T]＝impulse(sys,t)或impulse(sys)功能：返回系統(tǒng)的響應(yīng)Y和時間向量T，自動選擇仿真的時間范圍。其中sys可為系統(tǒng)傳遞函數(shù)、零極增益模型或狀態(tài)空間模型。%file92.mb=[23];a=[10.41];sys=tf(b,a);[y,t,x]=impulse(sys);plot(t,y)%file93.ma=[-0.55,-0.78;0.78,0];b=[1;0];c=[1.96,6.45];d=[0];impulse(a,b,c,d);title('LTI系統(tǒng)的沖激響應(yīng)')例：有一個二階系統(tǒng)：求系統(tǒng)的單位沖激響應(yīng)。(b)連續(xù)系統(tǒng)的單位階躍響應(yīng) 調(diào)用格式：格式：[Y,T]=step(sys,t)功能：返回系統(tǒng)的單位階躍響應(yīng)Y和仿真時間T，仿真時間范圍為t，或不設(shè)則可自動選擇。其中sys可為系統(tǒng)傳遞函數(shù)(tf)、零極增益模型(zpk)或狀態(tài)空間模型(ss)。%file94.mb=[23];a=[10.41];step(b,a)(c)連續(xù)系統(tǒng)的零輸入響應(yīng) 調(diào)用格式：[y,t,x]=initial(sys,t,x0)計算出連續(xù)時間系統(tǒng)LT1由于初始狀態(tài)x0所引起的零輸入響應(yīng)y，其中x為狀態(tài)記錄，t為仿真所用的采樣時間向量。(d)對任意輸入的連續(xù)LTI系統(tǒng)響應(yīng)函數(shù)lsim() 格式：[y,x]=lsim(a,b,c,d,u,t)功能：返回連續(xù)LTI系統(tǒng)對任意輸入時，系統(tǒng)的輸出響應(yīng)y和狀態(tài)記錄x，其中u給出每個輸入的時序列，一般情況下u為一個矩陣；t用于指定仿真的時間軸，它應(yīng)為等間隔。%file95.mnum=[2,5,1];den=[1,2,3];t=0:0.1:10peiod=4;u=(rem(t,peiod)>=peiod/2)lsim(num,den,u,t);title('方波響應(yīng)');

例：有一個二階系統(tǒng)：求當(dāng)輸入是周期為4s的方波時的輸出響應(yīng)。(e)離散系統(tǒng)(數(shù)字濾波器)的單位沖激響應(yīng) 調(diào)用格式：[y,x]＝impz(num,den,n)功能：返回單位沖激響應(yīng)y向量和時間狀態(tài)歷史記錄x向量。%file96.mnum=[0.18080.10470.31060.1047];den=[1.0000-1.14801.5107-0.6991]tf(num,den,0.1,’variable’,’z^-1’)impz(num,den,50);%取50個采樣點(diǎn)%也可用以下方法實(shí)現(xiàn)%file97.mn=50;imp=[1zeros(1,n-1)];y=filter(num,den,imp);stem(y);%file98.mclearall;num=[2,-3.5,1.5];den=[1,-1.7,0.3];dimpulse(num,den,10);title('離散LTI系統(tǒng)的沖激響應(yīng)')例：有一個二階系統(tǒng)：求系統(tǒng)的單位沖激響應(yīng)。(f)離散系統(tǒng)(數(shù)字濾波器)的單位階躍響應(yīng) 調(diào)用格式：[y,x]=stepz(num,den)/dstep功能：返回多項(xiàng)式傳遞函數(shù)G(z)=num(z)/den(z)表示的系統(tǒng)單位階躍響應(yīng)。%file99.mb=[12];a=[121];stepz(b,a)(g)離散系統(tǒng)(數(shù)字濾波器)的零輸入響應(yīng) 調(diào)用格式：[y,x,n]=dinitial(a,b,c,d,x0)功能：計算離散時間LTI系統(tǒng)由初始狀態(tài)x0所引起的零輸入響應(yīng)y和狀態(tài)響應(yīng)響應(yīng)x，取樣點(diǎn)數(shù)由函數(shù)自動選取。n為仿真所用的點(diǎn)數(shù)。(h)對任意輸入的離散LTI系統(tǒng)響應(yīng)函數(shù)dlsim()調(diào)用格式：[y,x]=dlsim(a,b,c,d,u)功能：返回離散LTI系統(tǒng)對輸入序列u的響應(yīng)y和狀態(tài)記錄x。%file100.mnum=[2,-3.4,5.5];den=[1,-1.2,0.8];u=randn(1,100);dlsim(num,den,u);title('隨機(jī)噪聲響應(yīng)')

例：有一個二階系統(tǒng)：求系統(tǒng)對100點(diǎn)隨機(jī)噪聲的響應(yīng)曲線。(i)濾波函數(shù)filter()格式：y=filter(B,A,x)功能：對向量x中的數(shù)據(jù)進(jìn)行濾波處理，即差分方程求解，產(chǎn)生輸出序列向量y。B和A分別為數(shù)字濾波器系統(tǒng)函數(shù)H(z)的分子和分母多項(xiàng)式系數(shù)向量。連續(xù)系統(tǒng)的頻率響應(yīng)H(jω)的定義與物理意義：

Yf(jω)=H(jω)F(jω)LT1系統(tǒng)把頻譜為F(jω)的輸入改變成頻譜為H(jω)F(jω)的響應(yīng)，改變的規(guī)律完全由H(jω)決定。H(jω)稱為系統(tǒng)的頻率響應(yīng)，定義為：幅度響應(yīng)相位響應(yīng)H(jω)與h(t)的關(guān)系：H(jω)反映了系統(tǒng)對輸入信號不同頻率分量的傳輸特性。由H(jω)的定義，顯然有

H(jω)=F[h(t)]即H(jω)等于系統(tǒng)沖激響應(yīng)h(t)的Fourier變換。1．求模擬濾波器Ha(s)的頻率響應(yīng)函數(shù)freqs()格式一：H＝freqs(B,A,W)

功能：計算由向量W(rad/s)指定的頻率點(diǎn)上模擬濾器系統(tǒng)函數(shù)Ha(s)的頻率響應(yīng)Ha(jΩ)，結(jié)果存于H向量中。格式二：[H,w]＝freqs(B,A,M)功能：計算出M個頻率點(diǎn)上的頻率響應(yīng)存于向量H中,M個頻率存放在向量w中。freqs自動將這M個頻率點(diǎn)設(shè)置在適當(dāng)?shù)念l率范圍。默認(rèn)w和M時,freqs自動選取200個頻率點(diǎn)計算。求該模擬濾波器的頻率響應(yīng)。MATLAB源程序如下：%file101.mB=1;A=[12.61313.41422.61311];W=0:0.1:2*pi*5;freqs(B,A,W)已知某模擬濾波器的系統(tǒng)函數(shù)

已知某濾波器的系統(tǒng)函數(shù)為2．求數(shù)字濾波器H(z)的頻率響應(yīng)函數(shù)freqz()格式一：H=freqz(B,A,W)功能：計算由向量W(rad)指定的數(shù)字頻率點(diǎn)上(通常指在H(z)的頻率響應(yīng)H(ejw)。

格式二：[H,W]=freqz(B,A,N)格式三：[H,W]=freqz(B,A,N,‘whole’)求該濾波器的頻率響應(yīng)。MATLAB源程序?yàn)椋?file102.mB=[10000000-1];A=1;freqz(B,A)該程序運(yùn)行所繪出的幅頻與相頻性曲線如圖所示。例：已知一系統(tǒng)傳遞函數(shù)為：求系統(tǒng)的單位沖激響應(yīng)h(n),以及h(n)的幅頻和相頻響應(yīng)。%file103.mN=128;x=[1,zeros(1,N-1)];%產(chǎn)生單位沖激函數(shù)num=[0.01-0.030.06-0.030.01];den=[12.531.80.5];y=dlsim(num,den,x);%計算單位沖激響應(yīng)figure(1);n=1:N;stem(n,y);gridon;title('單位沖激響應(yīng)');figure(2)freqz(num,den,N);%繪制幅頻和相頻響應(yīng)曲線gridon;dimpulse(num,den)單位沖激響應(yīng)函數(shù)，不帶輸出變量時，直接繪出單位沖激響應(yīng)幅頻和相頻響應(yīng)沖激響應(yīng)3．濾波函數(shù)filter從頻域角度，無論是連續(xù)時間LTI系統(tǒng)還是離散時間LTI系統(tǒng)，系統(tǒng)對輸入信號的響應(yīng)實(shí)質(zhì)上就是對輸入信號頻譜進(jìn)行不同選擇處理的過程，這個過程稱為濾波。格式：y=filter(B,A,x)功能：對向量x中的數(shù)據(jù)進(jìn)行濾波處理，即差分方程求解，產(chǎn)生輸出序列向量y。B和A分別為數(shù)字濾波器系統(tǒng)函數(shù)H(z)的分子和分母多項(xiàng)式系數(shù)向量。設(shè)系統(tǒng)的差分方程為求該系統(tǒng)對信號的響應(yīng)。MATLAB源程序?yàn)椋篺ile104.mB=1;A=[1,-0.8];n=0:31;x=0.8.^n;y=filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,2);stem(y)該程序運(yùn)行所得結(jié)果如圖所示時域信號特性頻域特性變換名稱非周期連續(xù)信號非周期連續(xù)頻譜傅里葉變換周期性連續(xù)信號非周期離散頻譜傅里葉級數(shù)非周期離散信號周期連續(xù)頻譜序列傅里葉變換周期離散信號周期離散頻譜離散傅里葉級數(shù)離散信號（有限樣本點(diǎn)）周期離散頻譜離散傅里葉變換傅里葉變換：是建立以時間為自變量的“信號”與以頻率為自變量的“頻譜函數(shù)”之間的某種對應(yīng)關(guān)系。遵循以下原則：時域連續(xù)周期離散非周期非周期離散周期連續(xù)頻域1.連續(xù)時間、連續(xù)頻率－傅里葉變換正變換：

逆變換：

這說明求和的問題可以用x(t)行向量乘以列向量來實(shí)現(xiàn)。式中是t的增量，在程序中用dt表示.例：求如圖矩形脈沖信號f(t)在-40~40rad/s區(qū)間的頻譜。%file105.mclc;clearall;N=input('取時間分隔的點(diǎn)數(shù)N=');dt=10/N;t=[1:N]*dt;%給出時間分割f=[ones(1,N/2),zeros(1,N/2)];%給出信號(方波)wf=input('需求的頻譜寬度wf=');Nf=input('需求的頻譜點(diǎn)數(shù)Nf=');w1=linspace(0,wf,Nf);F1=f*exp(-j*t'*w1)*dt;%求傅里葉變換w=[-fliplr(w1),w1(2:Nf)];%補(bǔ)上負(fù)頻率F=[fliplr(F1),F1(2:Nf)];%補(bǔ)上負(fù)頻率對應(yīng)的頻譜subplot(121)plot(t,f,'linewidth',1.5)gridon;axis([0,10,0,1.1]);subplot(122)plot(w,abs(F),'linewidth',1.5);gridon;取時間分隔的點(diǎn)數(shù)N=64需求的頻譜寬度wf=40需求的頻譜點(diǎn)數(shù)Nf=256取時間分隔的點(diǎn)數(shù)N=256需求的頻譜寬度wf=40需求的頻譜點(diǎn)數(shù)Nf=64采樣較密采樣稀，有頻率泄漏2連續(xù)時間、離散頻率－傅里葉級數(shù)正變換：

逆變換：

式中為離散頻率相鄰兩譜線間的角頻率間隔，k為諧波序號。3離散時間、離散頻率－離散傅里葉級數(shù)正變換：

逆變換：

5離散時間、離散頻率－離散傅里葉變換(DFT)正變換：

逆變換：

4時間離散、連續(xù)頻率－序列傅里葉變換正變換：

逆變換：

注意，對于序列傅里葉變換，如果x(n)為無限長，那么就不能用MATLAB直接利用上式計算，只可以用它對表達(dá)式在頻率點(diǎn)上求值，再畫出它的幅度和相位。如果x(n)為有限長，就可直接用MATLAB計算5離散時間、離散頻率－離散傅里葉變換(DFT)在MATLAB中，旋轉(zhuǎn)因子矩陣可以由dftmtx函數(shù)實(shí)現(xiàn)。dftmtx函數(shù):調(diào)用格式：W=dftmtx(n)功能：函數(shù)調(diào)用后，返回一個n*n的離散傅里葉變換DFT的旋轉(zhuǎn)因子矩陣W。這樣對于給定長度為n的行向量信號x(n),利用Xk=x*W就可以獲得x(n)的傅里葉變換X(k)。其實(shí)，旋轉(zhuǎn)因子矩陣也可以自己通過表達(dá)式構(gòu)造。例如：編寫函數(shù)，首先利用MATLAB計算序列x(n)的DFT和IDFT函數(shù).functionxn=idft(xk)N=length(xk);n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);WNnk=WN.^(-n'*k);xn=xk*WNnk/N;IDFT:functionxk=dft(xn)N=length(xn);WNnk=dftmtx(N);xk=xn*WNnk;DFT:WNnk=conj(dftmtx(N))DFT參數(shù)對頻率分辨率的影響：頻率分辨率就是，在頻率軸上能分辨出的兩個頻率點(diǎn)的最小間隔要能區(qū)分兩個頻率點(diǎn)f1,f2,，有效數(shù)據(jù)長度N必須滿足：可以對序列進(jìn)行補(bǔ)零操作，以增加數(shù)據(jù)的長度，這樣做雖然不能提高頻率分辨率，但補(bǔ)零后對原來的X(k)起到了插值作用，可以平滑頻譜的包絡(luò)。但是，如果增加有效數(shù)據(jù)的長度，則可以提高頻率分辨率。利用上面編寫的DFT函數(shù)，求的頻譜特性，其中取fs=20Hz%file106.mN=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);k=(0:N/2-1)*fs/Nsubplot(211);plot(n,xn);gridon;title('x(n)(0=<n<128)')subplot(212);plot(k,abs(xk(1:N/2)));gridon;%file107.mN=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);N1=256;n1=0:N1-1;xn1=[xn,zeros(1,N1-N)];xk1=dft(xn1);m=(0:N1/2-1)*fs/N1subplot(211);plot(n1,xn1);gridon;title('x(n)補(bǔ)零后');subplot(212);plot(m,abs(xk1(1:N1/2)));gridon;title(‘x(n)補(bǔ)零后’)%file108.mN=402;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);k=(0:N/2-1)*fs/Nsubplot(211);plot(n,xn);gridon;title('x(n)(0=<n<401)')subplot(212);plot(k,abs(xk(1:N/2)));gridon;

%file109.mN=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);subplot(211)plot(abs(xk));gridon;subplot(212)k=(0:N-1)*fs/N;plot(k,abs(xk));gridon;(1)．一維快速正傅里葉變換函數(shù)fft格式：X=fft(x,N)功能：采用FFT算法計算序列向量x的N點(diǎn)DFT變換，當(dāng)N缺省時，fft函數(shù)自動按x的長度計算DFT。當(dāng)N為2整數(shù)次冪時，fft按基-2算法計算，否則用混合算法。(2)．一維快速逆傅里葉變換函數(shù)ifft格式：x=ifft(X,N)功能：采用FFT算法計算序列向量X的N點(diǎn)IDFT變換。(3).將零頻分量移至頻譜中心的函數(shù)格式：Y=fftshift(X)功能：用來重新排列X=fft(x)的輸出，把X的左右兩半進(jìn)行交換，從而將零頻分量移至頻譜中心。%file110.mfs=1000;

%采樣頻率t=0:1/fs:0.5;%時間軸取值范圍及間隔x=cos(160*pi*t);N=512;

%FFT點(diǎn)數(shù)subplot(311)；plot(t,x);gridon;y=fft(x,N);

%N點(diǎn)FFT[x(t)]xk=fftshift(abs(y));

%求幅度譜并將零頻分量移至頻譜中心位置f=(-N/2:N/2-1)*fs/512;

%頻率范圍subplot(312)plot(f,xk);%做頻譜圖gridon;%用IFFT恢復(fù)原始信號z=real(ifft(y,N));ti=[0:length(z)-1]/fs;subplot(313)plot(ti,z);gridon;例：求的FFT序列和IFFT序列用快速傅里葉變換FFT計算下面兩個序列的卷積。MATLAB程序(部分)：file111.m%線性卷積xn=sin(0.4*[1:15]);

%對序列x(n)賦值,M=15hn=0.9.^(1:20);

%對序列h(n)賦值,N=20ticyn=conv(xn,hn); %直接調(diào)用函數(shù)conv計算卷積tocsubplot(2,2,1);stem(xn);subplot(2,2,2);stem(hn);subplot(2,2,3);stem(yn);h=get(gca,'position');h(3)=2.3*h(3)set(gca,'position',h)%園周卷積L=pow2(nextpow2(M+N-1));

ticXk=fft(xn,L);

Hk=fft(hn,L);

Yk=Xk.*Hk;

yn=ifft(Yk,L);

toc

x(n),h(n)及其線性卷積波形Elapsedtimeis0.005074seconds.Elapsedtimeis0.000072seconds.例題部分例1求系統(tǒng)

y"(t)+2y'(t)+100y(t)=10f(t)的零狀態(tài)響應(yīng)，已知f(t)=sin(2pt)u(t)。%file112.mts=0;te=5;dt=0.01;sys=tf([1],[12100]);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求系統(tǒng)

y"(t)+2y'(t)+100y(t)=10f(t)的零狀態(tài)響應(yīng)，已知f(t)=d(t)。%連續(xù)時間系統(tǒng)的沖激響應(yīng)file113.m

ts=0;te=5;dt=0.01;sys=tf([10],[12100]);t=ts:dt:te;y=impulse(sys,t);plot(t,y);xlabel('Time(sec)')ylabel('h(t)')例3

分析噪聲干擾的信號f[k]=s[k]+d[k]通過M點(diǎn)滑動平均系統(tǒng)的響應(yīng)，

其中s[k]=(2k)0.9k是原始信號，d[k]是噪聲。File114.mR=51;d=rand(1,R)-0.5;k=0:R-1;s=2*k.*(0.9.^k);f=s+d;figure(1);plot(k,d,'r-.',k,s,'b--',k,f,'g-');M=5;b=ones(M,1)/M;a=1;y=filter(b,a,f);figure(2);plot(k,s,'b--',k,y,'r-');噪聲干擾信號f[k]=s[k]+d[k]通過M點(diǎn)滑動平均系統(tǒng)的響應(yīng)例4求系統(tǒng)y[k]+3y[k-1]+2y[k-2]=10f[k]的單位脈沖響應(yīng)。%

離散系統(tǒng)的單位脈沖響應(yīng)file115.m

k=0:10;a=[132];b=[1];h=impz(b,a,k);stem(k,h)根據(jù)卷積公式：因此編程的過程為：（1）寫出h(t)＝e-0.1t的MATLAB表達(dá)式；（2）寫出u(t)的MATLAB表達(dá)式；（3）利用MATLAB的卷積語句y=conv(u，h)求解（4）畫曲線plot(t,y)。例5計算x[k]*y[k]并畫出卷積結(jié)果，已知x[k]={1,2,3,4;k=0,1,2,3}，

y[k]={1,1,1,1,1;k=0,1,2,3,4}。%file116.mx=[1,2,3,4];y=[1,1,1,1,1];z=conv(x,y);N=length(z);stem(0:N-1,z);例6方波分解為多次正弦波之和圖示的周期性方波，其傅里葉級數(shù)為分別計算直到9次諧波，并做圖。file118.m輸入周期性方波--x利用MATLAB進(jìn)行信號的頻域分析

周期信號頻譜的MATLAB實(shí)現(xiàn)

用數(shù)值積分分析非周期信號頻譜一、周期信號頻譜的MATLAB實(shí)現(xiàn)

頻譜Fn一般為復(fù)數(shù)，可分別利用abs和angle函數(shù)獲得其幅頻特性和相頻特性。其調(diào)用格式分別為x=abs(Fn)y=angle(Fn)

周期信號的頻譜Fn為離散信號，可以用stem例1試用MATLAB畫出圖示周期三角波信號的頻譜。解：周期信號的頻譜為畫三角波信號頻譜的MATLAB程序%file119.mN=8;n1=-N:-1;%計算n=-N到-1的Fourier系數(shù)

c1=-4*j*sin(n1*pi/2)/pi^2./n1.^2;c0=0;%計算n=0時的Fourier系數(shù)n2=1:N;%計算n=1到N的Fourier系數(shù)

c2=-4*j*sin(n2*pi/2)/pi^2./n2.^2;cn=[c1c0c2];n=-N:N;subplot(2,1,1);stem(n,abs(cn));ylabel('Cn的幅度');subplot(2,1,2);stem(n,angle(cn));ylabel('Cn的相位');xlabel('\omega/\omega0');程序運(yùn)行結(jié)果例2求周期矩形脈沖的Fourier級數(shù)表示式。并用MATLAB求出由前N項(xiàng)Fourier級數(shù)系數(shù)得出的信號近似波形。取A=1,T=2,t=1,w0=p

%Gibbsphenomenon%file120.mt=-2:0.001:2;N=input('Numberofharmonics=');c0=0.5;cN=c0*ones(1,length(t));

%dccomponent

forn=0:2:N-1%evenharmonicsarezero

cN=cN+cos(pi*n*t).*sinc(n*pi/2);

endplot(t,cN);%GibbsphenomenonN=5N=15N=50N=500二、用數(shù)值積分分析非周期信號頻譜數(shù)值函數(shù)積分quad8可用來計算非周期信號頻譜F

是一個字符串，它表示被積函數(shù)的文件名。a,b

分別表示定積分的下限和上限quad8的返回是用自適應(yīng)Simpson算法得出的積分值y=quad8('F',a,b)例3試用數(shù)值方法近似計算三角波信號的頻譜F(jw)=Sa2(w

/2)解:圖示三角波可表示為三角波信號頻譜的理論值為例3試用數(shù)值方法近似計算三角波信號的頻譜functiony=sf1(t,w);

y=(t>=-1&t<=1).*(1-abs(t)).*exp(-j*w*t);%file121.mw=linspace(-6*pi,6*pi,512);N=length(w);F=zeros(1,N);fork=1:NF(k)=quad('sf1',-1,1,[],[],w(k));end

figure(1);plot(w,real(F));title('')xlabel('\omega');ylabel('F(j\omega)');figure(2);plot(w,real(F)-sinc(w/2/pi).^2);xlabel('\omega');title('計算誤差');例3試用數(shù)值方法近似計算三角波信號的頻譜運(yùn)行結(jié)果利用MATLAB進(jìn)行系統(tǒng)頻域分析連續(xù)系統(tǒng)頻響特性的計算周期信號通過系統(tǒng)的響應(yīng)離散系統(tǒng)頻響特性的計算一、連續(xù)系統(tǒng)頻響特性的計算b

分子多項(xiàng)式系數(shù)

a

分母多項(xiàng)式系數(shù)

w

需計算的H(jw)的抽樣點(diǎn)

(數(shù)組w中少需包含兩個w的抽樣點(diǎn))。計算頻響的MATLAB函數(shù)H=freqs(b,a,w)一、連續(xù)系統(tǒng)頻響特性的計算例1三階歸一化的Butterworth低通濾波器的系統(tǒng)函數(shù)為%file122.mw=linspace(0,5,200);b=[1];a=[1221];h=freqs(b,a,w);subplot(2,1,1);plot(w,abs(h));subplot(2,1,2);plot(w,angle(h));

試畫出|H(jw)|和(w)。一、連續(xù)系統(tǒng)頻響特性的計算三階Butterworth低通濾波器的幅度響應(yīng)和相位響應(yīng)二、周期信號通過系統(tǒng)的響應(yīng)例2周期方波通過RC系統(tǒng)的響應(yīng)。

二、周期信號通過系統(tǒng)的響應(yīng)例2周期方波通過RC系統(tǒng)的響應(yīng)。%p5_2PeriodicsignalpassLTIsystemfile123.mT=4;w0=2*pi/T;RC=0.1;t=-6:0.01:6;N=51;c0=0.5;xN=c0*ones(1,length(t));

%dcforn=1:2:N%evenharmonicsarezero

H=abs(1/(1+j*RC*w0*n));phi=angle(1/(1+j*RC*w0*n));

xN=xN+H*cos(w0*n*t+phi)*sinc(n*0.5);

endplot(t,xN);xlabel(['timeRC=',num2str(RC)]);grid;set(gca,'xtick',[-5-3-10135]);二、周期信號通過系統(tǒng)的響應(yīng)例2周期方波通過RC系統(tǒng)的響應(yīng)。三、離散系統(tǒng)頻率響應(yīng)的計算計算頻率響應(yīng)的MATLAB函數(shù)b

分子的系數(shù)

a

分母系數(shù)w

抽樣的頻率點(diǎn)(至少2點(diǎn))，

w在0~2p之間幅頻特性:abs,相頻特性:angle

h=freqz(b,a,w)

b=[1];w=linspace(0,2*pi,512);h2=freqz(b,a2,w);plot(w/pi,abs(h1),w/pi,abs(h2),':');a1=[1-0.9];a2=[10.9];h1=freqz(b,a1,w);legend('\alpha=0.9','\alpha=-0.9');%file124.m三、離散系統(tǒng)頻響特性的計算利用MATLAB進(jìn)行連續(xù)系統(tǒng)的復(fù)頻域分析部分分式展開的MATLAB實(shí)現(xiàn)H(s)的零極點(diǎn)與系統(tǒng)特性的MATLAB計算一、部分分式展開的MATLAB實(shí)現(xiàn)[r,p,k]=residue(num,den)num,den分別為F(s)分子多項(xiàng)式和分母多項(xiàng)式的系數(shù)向量。r為部分分式的系數(shù)，p為極點(diǎn)，k為多項(xiàng)式二、H(s)的零極點(diǎn)與系統(tǒng)特性的MATLAB計算計算多項(xiàng)式根roots的函數(shù)可用于計算H(s)的零極點(diǎn)。r=roots(N)%計算多項(xiàng)式N的根H(s)零極點(diǎn)分布圖可用pzmap函數(shù)畫出，調(diào)用形式為pzmap(sys)表示畫出sys所描述系統(tǒng)的零極點(diǎn)圖。例1試畫出系統(tǒng)

的零極點(diǎn)分布圖，求其單位沖激響應(yīng)h(t)和頻率響應(yīng)H(j)，并判斷系統(tǒng)是否穩(wěn)定。%file125.mnum=[1];den=[1221];sys=tf(num,den);poles=roots(den)figure(1);pzmap(sys);t=0:0.02:10;h=impulse(num,den,t);figure(2);plot(t,h)title('ImpulseRespone')[H,w]=freqs(num,den);figure(3);plot(w,abs(H))xlabel('\omega')title('Magnitu')

利用MATLAB進(jìn)行離散系統(tǒng)的Z域分析部分分式展開的MATLAB實(shí)現(xiàn)H(z)的零極點(diǎn)與系統(tǒng)特性的MATLAB計算一、部分分式展開的MATLAB實(shí)現(xiàn)[r,p,k]=residuez(num,den)num,den分別為F(z)分子多項(xiàng)式和分母多項(xiàng)式的系數(shù)向量。r為部分分式的系數(shù)，p為極點(diǎn)，k為多項(xiàng)式的系數(shù)。若為真分式，則k為零。二、H(z)的零極點(diǎn)與系統(tǒng)特性的MATLAB計算利用t