matlab在信號處理中的應(yīng)用課件

第4章MATLAB在信號處理中的應(yīng)用4.1信號及其表示4.2信號的基本運算4.3信號的能量和功率4.4線性時不變系統(tǒng)4.5線性時不變系統(tǒng)的響應(yīng)4.6線性時不變系統(tǒng)的頻率響應(yīng)

4.7傅里葉(Fourier)變換4.8IIR數(shù)字濾波器的設(shè)計方法4.9FIR數(shù)字濾波器設(shè)計1可編輯課件PPT4.1信號及其表示4.1.1連續(xù)時間信號的表示

在MATLAB中，對連續(xù)時間信號用采樣點的數(shù)據(jù)來表示，當(dāng)采樣點很密時可看成連續(xù)信號。連續(xù)時間信號：時間變化連續(xù)。如y=x(t)離散時間信號(序列)：時間離散，如x(nT)=x(t)|t=nT.當(dāng)T足夠小時，可以認(rèn)為是連續(xù)信號。例如：用MATLAB命令繪出關(guān)于t的曲線，t的范圍為0~30s,并以0.1遞增。2可編輯課件PPTt=0:0.1:30; %對時間變量賦值x=exp(-0.707*t).*sin(2/3.*t);%計算變量所對應(yīng)的函數(shù)值plot(t,x);%繪制函數(shù)曲線grid;%圖形窗口加網(wǎng)格xlabel('time(sec)');%標(biāo)注x軸,y軸ylabel('x(t)');axis([-0.05,30,-0.05,0.35]);%確定x軸，y軸顯示范圍3可編輯課件PPT函數(shù)名格式功能sawtoothSawtooth（t,width）產(chǎn)生鋸齒波或三角波信號（周期為2，幅值-1~1）squareSquare(t,duty)產(chǎn)生方波信號（周期為2，幅值-1~1）sincSinc(x)產(chǎn)生sinc函數(shù)波形chirpchirp(t,f0,t1,f1,’method’,phi)產(chǎn)生線性調(diào)頻掃頻信號gauspulsgauspuls（T,FC,BW,BWR）產(chǎn)生高斯正弦脈沖信號vcovoc(x,fc,fs)電壓控制振蕩器pulstranpulstran(t,d,’func’)產(chǎn)生沖激串rectpulerectpule(t,w)產(chǎn)生非周期的方波信號tripulstripuls(t,w)產(chǎn)生非周期的三角波信號diricdiric(x,n)產(chǎn)生Dirichlet或周期sinc函數(shù)gmonopulsgmonopuls(t,fc)產(chǎn)生高斯單脈沖信號4.1.2工具箱中的信號產(chǎn)生函數(shù)4可編輯課件PPTt=0:0.0001:1;x1=sawtooth(2*pi*50*t,0.8)plot(t,x1)axis([0,0.2,-1,1])t=0:0.0001:0.05;x1=square(2*pi*50*t,50)plot(t,x1)例如：產(chǎn)生周期為0.02的三角波產(chǎn)生頻率為50Hz，占空比為50%的方波5可編輯課件PPTt=0:0.001:2y=chirp(t,0,1,150)figure(1)plot(t,y)axis([0,0.5,0,1])產(chǎn)生線性調(diào)頻信號：6可編輯課件PPTt=0:0.01:1;y1=tripuls(t);y2=tripuls(t,0.6);subplot(211);plot(t,y1);subplot(212);plot(t,y2);產(chǎn)生非周期三角波信號：7可編輯課件PPT4.1.3離散時間信號的表示在MATLAB中，離散時間信號x(n)的表示：需用一個向量x表示序列幅值，用另一個等長的定位時間變量n，才能完整地表示一個序列。

[例4-10]繪制離散時間信號的棒狀圖。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=-1,x(4)=0。MATLAB源程序為：n=-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]')運行結(jié)果如圖4.11所示。圖4.11離散時間信號圖形8可編輯課件PPT4.1.4幾種常用離散時間信號的表示1．單位脈沖序列直接實現(xiàn)：x=zeros(1,N);x(1,n0)=1;2．單位階躍序列

直接實現(xiàn)：n=[ns:nf];x=[(n-n0)>=0];函數(shù)實現(xiàn)：function

[x,n]=stepseq(n0,ns,nf)n=[ns:nf];x=[(n-n0)>=0];9可編輯課件PPT3．實指數(shù)序列直接實現(xiàn)：n=[ns:nf];x=a.^n;4．復(fù)指數(shù)序列直接實現(xiàn)：n=[ns:nf];x=exp((sigema+jw)*n);5．正(余)弦序列直接實現(xiàn)：n=[ns:nf];x=cos(w*n+sita);復(fù)數(shù)求實部：real(x)復(fù)數(shù)求虛部：imag(x)復(fù)數(shù)求幅度：abs(x)復(fù)數(shù)求相角：angel(x)復(fù)數(shù)運算的常用函數(shù)：10可編輯課件PPTfunctionfsxl(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);%生成相位序列fsxl(0.1,pi/3,-10,10)11可編輯課件PPT4.2信號的基本運算4.2.1信號的相加與相乘

y(n)=x1(n)+x2(n)y(n)=x1(n)×x2(n)MATLAB實現(xiàn)：y=x1+x2;y=x1.*x2兩信號相加和相乘的MATLAB實現(xiàn)：首先將兩序列時間變量延拓到相同長度，然后再逐點相加相乘4.2.2序列移位與周期延拓運算序列移位：y(n)=x(n-m)。MATLAB實現(xiàn)：y=x;ny=nx-m序列周期延拓：y(n)=x((n))M，MATLAB實現(xiàn)：ny=nxs:nxf；y=x(mod(ny,M)+1)12可編輯課件PPTn1=[-5:4];%序列x1(n)的起始及終止位置n1s=-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;%逐點相加yp=y1.*y2;%逐點相乘subplot(411);stem(n,y1,'.');subplot(412);stem(n,y2,'.');subplot(413);stem(n,ya,'.');subplot(414);stem(n,yp,'.');13可編輯課件PPT14可編輯課件PPTN=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(411);stem(n,x,'.');ylabel('x(n)');subplot(412);stem(n,xm,'.');ylabel('x(n-3)');subplot(413);stem(n,xc,'.');ylabel('x((n))_8');subplot(414);stem(n,xcm,'.');ylabel('x((n-3))_8');15可編輯課件PPT16可編輯課件PPT4.2.4兩序列的卷積運算兩序列卷積運算：

MATLAB實現(xiàn)：y=conv(x1,x2)。序列x1(n)和x2(n)必須長度有限。

4.2.5兩序列的相關(guān)運算兩序列相關(guān)運算：

MATLAB實現(xiàn)：y=xcorr(x1,x2)。4.2.3序列翻褶與序列累加運算序列翻褶：y(n)=x(-n)。MATLAB可實現(xiàn):y=fliplr(x)序列累加的數(shù)學(xué)描述為：

MATLAB實現(xiàn)：y=cumsum(x)17可編輯課件PPT例如：求序列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(221);stem(n,x);xlabel('n');ylabel('x(n)');subplot(223);stem(n1,y);xlabel('n');ylabel('y(n)=x(-n)');subplot(222);stem(n2,y);xlabel('n');ylabel('y(n)=x(-3+n)');subplot(224);stem(n,s);xlabel('n');ylabel('s(n)');18可編輯課件PPT4.3信號的能量和功率1.信號能量數(shù)字定義：MATLAB實現(xiàn):E=sum(x.*conj(x));或E=sum(abs(x).^2);數(shù)字定義：2.信號功率MATLAB實現(xiàn):

P=sum(x.*conj(x))/N;或E=sum(abs(x).^2)/N;19可編輯課件PPT非周期三角波信號能量的MATLAB計算：dt=0.0001;t=0:dt:1;x=tripuls(t);E=sum(abs(x).^2*dt)dt=0.001;t=0:dt:2*pi;x=cos(t);P=sum(abs(x).^2*dt)./(2*pi)E=0.1667P=0.5001余弦信號的平均功率MATLAB計算：20可編輯課件PPT4.4線性時不變系統(tǒng)一個信號處理系統(tǒng)就是將輸入信號變換成輸出信號的運算過程，如圖所示。在此過程中，輸出的信號稱為系統(tǒng)對輸入信號作出的響應(yīng)，輸入信號稱為系統(tǒng)的激勵信號。當(dāng)一個系統(tǒng)具有可加性和齊次性，則該系統(tǒng)稱為線性系統(tǒng)。如果系統(tǒng)響應(yīng)與激勵加于系統(tǒng)的時刻無關(guān)，則稱該系統(tǒng)為時不變系統(tǒng)。21可編輯課件PPT4.4線性時不變系統(tǒng)4.4.1系統(tǒng)的描述1．常系數(shù)線性微分/差分方程2．系統(tǒng)傳遞函數(shù)連續(xù)系統(tǒng)：

離散系統(tǒng)：

連續(xù)系統(tǒng)：

離散系統(tǒng)：

22可編輯課件PPT注意：在MATLAB中，傳遞函數(shù)由分子、分母兩個多項式的系數(shù)來表示，系數(shù)為降冪排列。例如：系統(tǒng)傳遞函數(shù)的MATLAB表示：可以表示為：num=[1,0.2,1];den=[1,0.5,1]3．零－極點增益模型連續(xù)系統(tǒng)：

離散系統(tǒng)：

23可編輯課件PPT4．極點留數(shù)模型離散系統(tǒng)：

連續(xù)系統(tǒng)：

5．二次分式模型連續(xù)系統(tǒng)：

離散系統(tǒng)：

6．狀態(tài)空間模型連續(xù)系統(tǒng)：

離散系統(tǒng)：

24可編輯課件PPT系統(tǒng)模型的MATLAB表示系統(tǒng)模型表示名稱參數(shù)傳遞函數(shù)模型tfnum---傳遞函數(shù)的分子系數(shù)向量den---傳遞函數(shù)的分母系數(shù)向量零極點增益模型zpz---系統(tǒng)零點向量；p---系統(tǒng)極點向量k---增益系數(shù)二次分式模型sossos---二次分式的系數(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ù)模型和零－極點增益模型。其中sos表示二次分式，g為比例系數(shù)，sos為L×6的矩陣，即

25可編輯課件PPT4.4.2系統(tǒng)模型的轉(zhuǎn)換函數(shù)1．ss2tf函數(shù)

格式：[num,den]=ss2tf(A,B,C,D,iu)功能：將指定輸入量iu的線性系統(tǒng)(A,B,C,D)狀態(tài)空間模型轉(zhuǎn)換為傳遞函數(shù)模型[num,den]。2．zp2tf函數(shù)格式：[num,den]=zp2tf(z,p,k)功能：將給定系統(tǒng)的零－極點增益模型轉(zhuǎn)換為傳遞函數(shù)模型，z、p、k分別為零點列向量、極點列向量和增益系數(shù)。3．tf2sos函數(shù)格式：[sos,g]=tf2sos

(num,den)功能：將給定系統(tǒng)的傳遞函數(shù)模型[num,den]轉(zhuǎn)換為二次分式模型sos,g為增益系數(shù)。26可編輯課件PPT線性系統(tǒng)模型的變換函數(shù)函數(shù)名功能說明函數(shù)名功能說明ss2tf狀態(tài)空間模型轉(zhuǎn)換為傳遞函數(shù)模型zp2tf零－極點增益模型轉(zhuǎn)換為傳遞函數(shù)模型ss2zp狀態(tài)空間模型轉(zhuǎn)換為零－極點增益模型zp2ss零－極點增益模型轉(zhuǎn)換為狀態(tài)空間模型ss2sos狀態(tài)空間模型轉(zhuǎn)換為二次分式模型zp2sos零－極點增益模型轉(zhuǎn)換為二次分式模型tf2ss傳遞函數(shù)模型轉(zhuǎn)換為狀態(tài)空間模型sos2tf二次分式模型轉(zhuǎn)換為傳遞函數(shù)模型tf2zp傳遞函數(shù)模型轉(zhuǎn)換為零－極點增益模型sos2zp二次分式模型轉(zhuǎn)換為零－極點增益模型tf2sos傳遞函數(shù)模型轉(zhuǎn)換為二次分式模型sos2ss二次分式模型轉(zhuǎn)換為狀態(tài)空間模型27可編輯課件PPT在命令窗口輸入：>>z=[3];p=[1,2];k=2;>>[A,B,C,D]=zp2ss(z,p,k)屏幕顯示結(jié)果：A=3.0000-1.41421.41420B=10C=2.0000-4.2426D=0[例4-18]將系統(tǒng)轉(zhuǎn)換為狀態(tài)空間模型[A,B,C,D]28可編輯課件PPT[例4-19]求離散時間系統(tǒng)的零、極點向量和增益系數(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=229可編輯課件PPT4.4.3系統(tǒng)互聯(lián)與系統(tǒng)結(jié)構(gòu)MATLAB實現(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)30可編輯課件PPTMATLAB實現(xiàn)函數(shù)parallel()格式：[A,B,C,D]=parallel(A1,B1,C1,D1,A2,B2,C2,D2)或

[num,den]=parallel(num1,den1,num2,den2)2.系統(tǒng)的并聯(lián)將系統(tǒng)1、系統(tǒng)2并聯(lián)，可得到并聯(lián)連接的傳遞函數(shù)形式為：31可編輯課件PPT3.兩個系統(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)；當(dāng)sig=+1時表示正反饋；當(dāng)sig=-1時表示負(fù)反饋。32可編輯課件PPT[例4-20]求兩個單輸入單輸出子系統(tǒng)的級聯(lián)、并聯(lián)和反饋后系統(tǒng)的傳遞函數(shù)。MATLAB源程序為：

num1=1;den1=[1,1];%系統(tǒng)1

num2=2;den2=[1,2];%系統(tǒng)2

[nums,dens]=series(num1,den1,num2,den2)%實現(xiàn)兩個系統(tǒng)級聯(lián)

[nump,denp]=parallel(num1,den1,num2,den2)%實現(xiàn)兩個系統(tǒng)并聯(lián)

[numf,denf]=feedback(num1,den1,num2,den2)%實現(xiàn)兩個系統(tǒng)反饋程序運行結(jié)果為：nums=002；dens=132nump=034；denp=132numf=012；denf=134因此，各系統(tǒng)的傳遞函數(shù)分別為：33可編輯課件PPT

在實際應(yīng)用中，也可以把一個復(fù)雜的線性時不變（LTI）系統(tǒng)分解為幾個簡單系統(tǒng)的組合結(jié)構(gòu)，即直接型結(jié)構(gòu)、級聯(lián)型結(jié)構(gòu)和并聯(lián)型結(jié)構(gòu)。Matlab所提供的系統(tǒng)模型變換函數(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)。[例4-20]已知FIR數(shù)字濾波器的傳遞函數(shù)為：求出其級聯(lián)型結(jié)構(gòu)和格型結(jié)構(gòu)clear;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)換MATLAB源程序：34可編輯課件PPT級聯(lián)型結(jié)構(gòu)系數(shù)：sos=1.00000.536001.0000001.00000.00570.62191.000000g=2格型結(jié)構(gòu)系數(shù)（反射系數(shù)）：K=0.25000.50000.333335可編輯課件PPTb=[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)MATLAB源程序：36可編輯課件PPT級聯(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=-1837可編輯課件PPT38可編輯課件PPT4.5線性時不變系統(tǒng)的響應(yīng)4.5.1線性時不變系統(tǒng)的時域響應(yīng)1．連續(xù)LTI系統(tǒng)的響應(yīng)2．離散LTI系統(tǒng)的響應(yīng)用MATLAB中的卷積函數(shù)conv()來實現(xiàn)。用MATLAB中的卷積函數(shù)conv()來實現(xiàn)。39可編輯課件PPTdt=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;40可編輯課件PPTx=ones(1,10);n=[0:14];h=0.5.^n;y=conv(x,h);stem(y);xlabel('n');ylabel('y[n]');41可編輯課件PPTnx=-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(221);stem(nx,x);title('x[n]');xlabel('n');xlabel('x[n]');subplot(222);stem(nh,h);title('h[n]');xlabel('n');xlabel('h[n]');subplot(223);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)改寫為：42可編輯課件PPT格式：[y,x]=lsim(a,b,c,d,u,t)功能：返回連續(xù)LTI系統(tǒng)

(2)對任意輸入的離散LTI系統(tǒng)響應(yīng)函數(shù)dlsim()格式：[y,x]=dlsim(a,b,c,d,u)功能：返回離散LTI系統(tǒng)

對任意輸入時系統(tǒng)的輸出響應(yīng)y和狀態(tài)記錄x，其中u給出每個輸入的時序列，一般情況下u為一個矩陣；t用于指定仿真的時間軸，它應(yīng)為等間隔。對輸入序列u的響應(yīng)y和狀態(tài)記錄x。3．時域響應(yīng)函數(shù)(1)對任意輸入的連續(xù)LTI系統(tǒng)響應(yīng)函數(shù)lsim()

43可編輯課件PPTnum=[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)');

44可編輯課件PPT45可編輯課件PPTnum=[2,-3.4,5.5];den=[1,-1.2,0.8]；u=randn(1,100);dlsim(num,den,u);title('隨機(jī)噪聲響應(yīng)')

46可編輯課件PPT4.5.2LTI系統(tǒng)的單位沖激響應(yīng)1.求連續(xù)LTI系統(tǒng)的單位沖激響應(yīng)函數(shù)impulse()格式：[Y,T]=impulse(sys)或impulse(sys)功能：返回系統(tǒng)的響應(yīng)Y和時間向量T，自動選擇仿真的時間范圍。其中sys可為系統(tǒng)傳遞函數(shù)、零極增益模型或狀態(tài)空間模型。2.求離散系統(tǒng)的單位沖激響應(yīng)函數(shù)dimpulse()格式：[y,x]=dimpulse(num,den)功能：返回多項式傳遞函數(shù)的單位沖激響應(yīng)y向量和時間狀態(tài)歷史記錄x向量。47可編輯課件PPTa=[-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)')48可編輯課件PPTnum=[2,-3.5,1.5];den=[1,-1.7,0.3];dimpulse(num,den);title('離散LTI系統(tǒng)的沖激響應(yīng)')49可編輯課件PPT4.5.3時域響應(yīng)的其它函數(shù)1.求連續(xù)LTI系統(tǒng)的零輸入響應(yīng)函數(shù)initial()格式：[y,t,x]=initial(a,b,c,d,x0)功能：計算出連續(xù)時間LTI系統(tǒng)由于初始狀態(tài)x0所引起的零輸入響應(yīng)y。其中x為狀態(tài)記錄，t為仿真所用的采樣時間向量。2.求離散系統(tǒng)的零輸入響應(yīng)函數(shù)dinitial()格式：[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，取樣點數(shù)由函數(shù)自動選取。n為仿真所用的點數(shù)。3.求連續(xù)系統(tǒng)的單位階躍響應(yīng)函數(shù)step()格式：[Y,T]=step(sys)

功能：返回系統(tǒng)的單位階躍響應(yīng)Y和仿真所用的時間向量T，自動選擇仿真的時間范圍。其中sys可為系統(tǒng)傳遞函數(shù)(TF)、零極增益模型(ZPK)或狀態(tài)空間模型(SS)。4.求離散系統(tǒng)的單位階躍響應(yīng)函數(shù)dstep()格式：[y,x]=dstep(num,den)功能：返回多項式傳遞函數(shù)G(z)=num(z)/den(z)表示的系統(tǒng)單位階躍響應(yīng)。50可編輯課件PPT4.6線性時不變系統(tǒng)的頻率響應(yīng)51可編輯課件PPT52可編輯課件PPT4.6線性時不變系統(tǒng)的頻率響應(yīng)1．求模擬濾波器Ha(s)的頻率響應(yīng)函數(shù)freqs()格式一：H＝freqs(B,A,W)

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

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

格式二：[H,W]=freqz(B,A,N)格式三：[H,W]=freqz(B,A,N,‘whole’)求該濾波器的頻率響應(yīng)。MATLAB源程序為：B=[10000000–1];A=1;freqz(B,A)56可編輯課件PPT圖4.31濾波器幅度和相位曲線

該程序運行所繪出的幅頻與相頻性曲線如圖4.31所示。57可編輯課件PPT例：已知一系統(tǒng)傳遞函數(shù)為：求系統(tǒng)的單位沖激響應(yīng)h(n),以及h(n)的幅頻和相頻響應(yīng)。N=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)58可編輯課件PPT單位沖激響應(yīng)幅頻和相頻響應(yīng)沖激響應(yīng)59可編輯課件PPT3．濾波函數(shù)filter從頻域角度，無論是連續(xù)時間LTI系統(tǒng)還是離散時間LTI系統(tǒng)，系統(tǒng)對輸入信號的響應(yīng)實質(zhì)上就是對輸入信號頻譜進(jìn)行不同選擇處理的過程，這個過程稱為濾波。格式：y=filter(B,A,x)功能：對向量x中的數(shù)據(jù)進(jìn)行濾波處理，即差分方程求解，產(chǎn)生輸出序列向量y。B和A分別為數(shù)字濾波器系統(tǒng)函數(shù)H(z)的分子和分母多項式系數(shù)向量。[例4-33]設(shè)系統(tǒng)差分方程為求該系統(tǒng)對信號的響應(yīng)。MATLAB源程序為：B=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)60可編輯課件PPT圖4.32系統(tǒng)對信號的響應(yīng)

該程序運行所得結(jié)果如圖4.32所示61可編輯課件PPT時域信號特性頻域特性變換名稱非周期連續(xù)信號非周期連續(xù)頻譜傅里葉變換周期性連續(xù)信號非周期離散頻譜傅里葉級數(shù)非周期離散信號周期連續(xù)頻譜序列傅里葉變換周期離散信號周期離散頻譜離散傅里葉級數(shù)離散信號（有限樣本點）周期離散頻譜離散傅里葉變換傅里葉變換：是建立以時間為自變量的“信號”與以頻率為自變量的“頻譜函數(shù)”之間的某種對應(yīng)關(guān)系。遵循以下原則：時域連續(xù)周期離散非周期非周期離散周期連續(xù)頻域4.7傅里葉(Fourier)變換62可編輯課件PPT63可編輯課件PPT64可編輯課件PPT65可編輯課件PPT4.7傅里葉(Fourier)變換4.7.1連續(xù)時間、連續(xù)頻率－傅里葉變換正變換：

逆變換：

這說明求和的問題可以用x(t)行向量乘以列向量來實現(xiàn)。式中是t的增量，在程序中用dt表示.例：求如圖矩形脈沖信號f(t)在-40~40rad/s區(qū)間的頻譜。66可編輯課件PPTclc;clearall;N=input('取時間分隔的點數(shù)N=');dt=10/N;t=[1:N]*dt;

%給出時間分割f=[ones(1,N/2),zeros(1,N/2)];

%給出信號(方波)wf=input('需求的頻譜寬度wf=');Nf=input('需求的頻譜點數(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;67可編輯課件PPT取時間分隔的點數(shù)N=64需求的頻譜寬度wf=40需求的頻譜點數(shù)Nf=256取時間分隔的點數(shù)N=256需求的頻譜寬度wf=40需求的頻譜點數(shù)Nf=64采樣較密采樣稀，有頻率泄漏68可編輯課件PPT4.7.2連續(xù)時間、離散頻率－傅里葉級數(shù)正變換：

逆變換：

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

逆變換：

69可編輯課件PPT4.7.5離散時間、離散頻率－離散傅里葉變換(DFT)正變換：

逆變換：

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

逆變換：

注意，對于序列傅里葉變換，如果x(n)為無限長，那么就不能用MATLAB直接利用上式計算，只可以用它對表達(dá)式在頻率點上求值，再畫出它的幅度和相位。如果x(n)為有限長，就可直接用MATLAB計算70可編輯課件PPT4.7.5離散時間、離散頻率－離散傅里葉變換(DFT)在MATLAB中，旋轉(zhuǎn)因子矩陣可以由dftmtx函數(shù)實現(xiàn)。71可編輯課件PPTdftmtx函數(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)。其實，旋轉(zhuǎn)因子矩陣也可以自己通過表達(dá)式構(gòu)造。例如：編寫函數(shù)，首先利用MATLAB計算序列x(n)的DFT和IDFT函數(shù).72可編輯課件PPTfunctionxn=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))73可編輯課件PPTDFT參數(shù)對頻率分辨率的影響：頻率分辨率就是，在頻率軸上能分辨出的兩個頻率點的最小間隔要能區(qū)分兩個頻率點f1,f2,，有效數(shù)據(jù)長度N必須滿足：可以對序列進(jìn)行補(bǔ)零操作，以增加數(shù)據(jù)的長度，這樣做雖然不能提高頻率分辨率，但補(bǔ)零后對原來的X(k)起到了插值作用，可以平滑頻譜的包絡(luò)。但是，如果增加有效數(shù)據(jù)的長度，則可以提高頻率分辨率。74可編輯課件PPT利用上面編寫的DFT函數(shù)，求的頻譜特性，其中取fs=20Hz75可編輯課件PPTN=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;76可編輯課件PPTN=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ǔ)零后’)77可編輯課件PPTN=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;

78可編輯課件PPTN=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;79可編輯課件PPT1．一維快速正傅里葉變換函數(shù)fft格式：X=fft(x,N)功能：采用FFT算法計算序列向量x的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點IDFT變換。3.將零頻分量移至頻譜中心的函數(shù)格式：Y=fftshift(X)功能：用來重新排列X=fft(x)的輸出，把X的左右兩半進(jìn)行交換，從而將零頻分量移至頻譜中心。,80可編輯課件PPTfs=1000;

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

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

%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序列81可編輯課件PPT82可編輯課件PPT[例4-36]用快速傅里葉變換FFT計算下面兩個序列的卷積。圖4.35快速卷積框圖83可編輯課件PPTMATLAB程序(部分)：%線性卷積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計算卷積toc%園周卷積L=pow2(nextpow2(M+N-1));

ticXk=fft(xn,L);

Hk=fft(hn,L);

Yk=Xk.*Hk;

yn=ifft(Yk,L);

toc

圖4.36x(n),h(n)及其線性卷積波形Elapsedtimeis0.005074seconds.Elapsedtimeis0.000072seconds.84可編輯課件PPT練習(xí)求x(t)在有噪聲和無噪聲條件下的FFT.(提示：噪聲由randn函數(shù)生成)1.85可編輯課件PPT4.8IIR數(shù)字濾波器的設(shè)計方法1.數(shù)字濾波器的頻率響應(yīng)函數(shù)幅度響應(yīng)：相位響應(yīng)：圖4.37理想低通、高通、帶通、帶阻數(shù)字濾波器幅度特性

86可編輯課件PPT2.濾波器的技術(shù)指標(biāo)幅度響應(yīng)指標(biāo)、相位響應(yīng)指標(biāo)

圖4.38數(shù)字低通濾波器的幅度特性

通帶要求：

阻帶要求：

通帶最大衰減：阻帶最小衰減：87可編輯課件PPT4.8.1沖激響應(yīng)不變法2.MATLAB信號處理工箱中的專用函數(shù)impinvar()：格式：[BZ,AZ]=impinvar(B,A,Fs)功能：把具有[B,A]模擬濾波器傳遞函數(shù)模型轉(zhuǎn)換成采樣頻率為Fs(Hz)的數(shù)字濾波器的傳遞函數(shù)模型[BZ,AZ]。采樣頻率Fs的默認(rèn)值為Fs=1。1.沖激響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的基本原理：[例4-37]MATLAB源程序如下：num=[1];%模擬濾波器系統(tǒng)函數(shù)的分子den=[1,sqrt(5),2,sqrt(2),1];%模擬濾波器系統(tǒng)函數(shù)的分母[num1,den1]=impinvar(num,den)%求數(shù)字低通濾波器的系統(tǒng)函數(shù)程序的執(zhí)行結(jié)果如下：num1=-0.00000.09420.21580.0311den1=1.0000-2.00321.9982-0.76120.106988可編輯課件PPTMATLAB信號處理工具箱中的專用雙線性變換函數(shù)bilinear()格式：[numd,dend]＝bilinear(num,den,Fs)功能：把模擬濾波器的傳遞函數(shù)模型轉(zhuǎn)換成數(shù)字濾波器的傳遞函數(shù)模型。4.8.2雙線性變換法雙線性變換利用頻率變換關(guān)系：

[例4-38]MATLAB源程序如下：

num=1;%模擬濾波器系統(tǒng)函數(shù)的分子

den=[1,sqrt(3),sqrt(2),1];%模擬濾波器系統(tǒng)函數(shù)的分母

[num1,den1]=bilinear(num,den,1)%求數(shù)字濾波器的傳遞函數(shù)運算的結(jié)果如下：num1=0.05330.15990.15990.0533den1=1.0000-1.33820.9193-0.154689可編輯課件PPT4.8.3IIR數(shù)字濾波器的頻率變換設(shè)計法1.IIR數(shù)字濾波器的頻率變換設(shè)計法的基本原理根據(jù)濾波器設(shè)計要求，設(shè)計模擬原型低通濾波器，然后進(jìn)行頻率變換，將其轉(zhuǎn)換為相應(yīng)的模擬濾波器(高通、帶通等)，最后利用沖激響應(yīng)不變法或雙線性變換法，將模擬濾波器數(shù)字化成相應(yīng)的數(shù)字濾波器。

圖4.39IIR數(shù)字濾波器MATLAB設(shè)計步驟流程圖

90可編輯課件PPT1．MATLAB的典型設(shè)計利用在MATLAB設(shè)計IIR數(shù)字濾波器可分以下幾步來實現(xiàn)(1)按一定規(guī)則將數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換為模擬低通濾波器的技術(shù)指標(biāo)；(2)根據(jù)轉(zhuǎn)換后的技術(shù)指標(biāo)使用濾波器階數(shù)函數(shù)，確定濾波器的最小階數(shù)N和截止頻率Wc；(3)利用最小階數(shù)N產(chǎn)生模擬低通濾波原型；(4)利用截止頻率Wc把模擬低通濾波器原型轉(zhuǎn)換成模擬低通、高通、帶通或帶阻濾波器；(5)利用沖激響應(yīng)不變法或雙線性不變法把模擬濾波器轉(zhuǎn)換成數(shù)字濾波器。[例4-39]設(shè)計一個數(shù)字信號處理系統(tǒng)，它的采樣率為Fs＝100Hz，希望在該系統(tǒng)中設(shè)計一個Butterworth型高通數(shù)字濾波器，使其通帶中允許的最小衰減為0.5dB，阻帶內(nèi)的最小衰減為40dB，通帶上限臨界頻率為30Hz，阻帶下限臨界頻率為40Hz。

91可編輯課件PPTMATLAB源程序設(shè)計如下：%把數(shù)字濾波器的頻率特征轉(zhuǎn)換成模擬濾波器的頻率特征

wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100;

[N,Wc]=buttord(wp,ws,rp,rs,'s');%選擇濾波器的最小階數(shù)

[Z,P,K]=buttap(N);%創(chuàng)建Butterworth低通濾波器原型

[A,B,C,D]=zp2ss(Z,P,K);%零－極點增益模型轉(zhuǎn)換為狀態(tài)空間模型

[AT,BT,CT,DT]=lp2hp(A,B,C,D,Wc);%實現(xiàn)低通向高通的轉(zhuǎn)變

[num1,den1]=ss2tf(AT,BT,CT,DT);%狀態(tài)空間模型轉(zhuǎn)換為傳遞函數(shù)模型%運用雙線性變換法把模擬濾波器轉(zhuǎn)換成數(shù)字濾波器

[num2,den2]=bilinear(num1,den1,100);

[H,W]=freqz(num2,den2);%求頻率響應(yīng)

plot(W*Fs/(2*pi),abs(H));grid;%繪出頻率響應(yīng)曲線

xlabel('頻率/Hz'); ylabel('幅值')程序運行結(jié)果如圖4.40所示。

92可編輯課件PPT2．MATLAB的直接設(shè)計圖4.39IIR數(shù)字濾波器MATLAB設(shè)計步驟流程圖

[例4-41]試設(shè)計一個帶阻IIR數(shù)字濾波器，其具體的要求是：通帶的截止頻率：wp1＝650Hz、wp2＝850Hz；阻帶的截止頻率：ws1＝700Hz、ws2＝800Hz；通帶內(nèi)的最大衰減為rp＝0.1dB；阻帶內(nèi)的最小衰減為rs＝50dB；采樣頻率為Fs＝2000Hz。MATLAB源程序設(shè)計如下：

wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000;

wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);%利用Nyquist頻率頻率歸一化

[N,wc]=ellipord(wp,ws,rp,rs,'z');%求濾波器階數(shù)