MATLAB實(shí)驗(yàn)二(修改)

1、實(shí)驗(yàn)二 信號(hào)的表示及其基本運(yùn)算一、實(shí)驗(yàn)?zāi)康?、掌握連續(xù)信號(hào)及其MATLAB實(shí)現(xiàn)方法；2、掌握離散信號(hào)及其MATLAB實(shí)現(xiàn)方法3、掌握離散信號(hào)的基本運(yùn)算方法，以及MATLAB實(shí)現(xiàn)4 熟悉應(yīng)用MATLAB實(shí)現(xiàn)求解系統(tǒng)響應(yīng)的方法4、了解離散傅里葉變換的MATLAB實(shí)現(xiàn)5、了解IIR數(shù)字濾波器設(shè)計(jì)6、了解FIR數(shù)字濾波器設(shè)計(jì)1二、實(shí)驗(yàn)設(shè)備計(jì)算機(jī)，Matlab軟件三、實(shí)驗(yàn)內(nèi)容（一）、 連續(xù)信號(hào)及其MATLAB實(shí)現(xiàn)1、 單位沖激信號(hào)例1.1：單位沖擊信號(hào)的MATLAB實(shí)現(xiàn)程序如下：t1=-4;t2=4;t0=0;dt=0.01;t=t1:dt:t2;n=length(t);x=zeros(1,n);x(1

2、,(-t0-t1)/dt+1)=1/dt;stairs(t,x);axis(t1,t2,0,1.2/dt);2、 任意函數(shù)例1.2：用MATLAB畫出如下表達(dá)式的脈沖序列3 單位階躍函數(shù)例1.3：用MATLAB實(shí)現(xiàn)單位階躍函數(shù)clear all;t=-0.5:0.001:1;t0=0;u=stepfun(t,t0);plot(t,u)axis(-0.5 1 -0.2 1.2)4 斜坡函數(shù)例1.4：用MATLAB實(shí)現(xiàn)g(t)=3(t-1)clear all;t=0:0.01:3;B=3;t0=1;u=stepfun(t,t0);n=length(t);for i=1:n u(i)=B*u(i)*

3、(t(i)-t0);endplot(t,u)axis(-0.2 3.1 -0.2 6.2)5 抽樣信號(hào) 抽樣信號(hào)Sa(t)=sin(t)/t在MATLAB中用 sinc 函數(shù)表示。定義為 t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft); grid on; axis(-10,10,-0.5,1.2); %定義畫圖范圍，橫軸，縱軸 title(抽樣信號(hào)) %定義圖的標(biāo)題名字6 指數(shù)函數(shù) 例1.5：用MATLAB實(shí)現(xiàn)7 正弦函數(shù)例1.6：用MATLAB實(shí)現(xiàn)正弦函數(shù)f(t)=3cos(10t+1)8 虛指數(shù)信號(hào)例 虛指數(shù)信號(hào) 調(diào)用格式是f=exp(j*w)

4、*t) t=0:0.01:15;w=pi/4;X=exp(j*w*t);Xr=real(X); %取實(shí)部 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(實(shí)部);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),ti

5、tle(模);subplot(2,2,4),plot(t,Xn),axis(0,15,-(max(Xn)+1),max(Xn)+1),title(相角); %subplot(m,n,i) 命令是建立m行n列畫圖窗口，并指定畫圖位置i9 復(fù)指數(shù)信號(hào)例 復(fù)指數(shù)信號(hào) 調(diào)用格式是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(實(shí)部)subplot(2,2,3),plot(t,imag(f),title(虛部) subplot(2,2,2),plot(t,abs(f),ti

6、tle(模)subplot(2,2,4),plot(t,angle(f),title(相角)（二）、離散信號(hào)及其MATLAB實(shí)現(xiàn)1、 單位沖激序列例2.1：用MATLAB產(chǎn)生64點(diǎn)的單位沖激序列clear all;N=64;x=zeros(1,N);x(1)=1;xn=0:N-1;stem(xn,x)axis(-1 65 0 1.1)2、 任意序列例2.2：用MATLAB畫出如下表達(dá)式的脈沖序列3、 單位階躍序列例2.3：用MATLAB實(shí)現(xiàn)單位階躍函數(shù)4、 斜坡序列例2.4：用MATLAB實(shí)現(xiàn)g(n)=3(n-4)點(diǎn)數(shù)為32的斜坡序列clear all;N=32;k=4B=3;t0=1;x=z

7、eros(1,k) ones(1,N-k);for i=1:N x(i)=B*x(i)*(i-k);endxn=0:N-1;stem(xn,x)axis(-1 32 0 90)5、 正弦序列例2.5：用MATLAB實(shí)現(xiàn)幅度A=3，頻率f=100，初始相位=1.2，點(diǎn)數(shù)為32的正弦信號(hào)6、 實(shí)指數(shù)序列例2.6：用MATLAB實(shí)現(xiàn)，點(diǎn)數(shù)為32的實(shí)指數(shù)序列clear all;N=32;A=3;a=0.7;xn=0:N-1;x=A*a.xn;stem(xn,x)7、 復(fù)指數(shù)序列例2.7：用MATLAB實(shí)現(xiàn)幅度A=3，a=0.7，角頻率=314，點(diǎn)數(shù)為32的實(shí)指數(shù)序列clear all;N=32;A=3

8、;a=0.7;w=314;xn=0:N-1;x=A*exp(a+j*w)*xn);stem(xn,x)8、 隨機(jī)序列利用MATLAB產(chǎn)生兩種隨機(jī)信號(hào)：rand(1,N)在區(qū)間上產(chǎn)生N點(diǎn)均勻分布的隨機(jī)序列randn(1,N)產(chǎn)生均值為0，方差為1的高斯隨機(jī)序列，即白噪聲序列例2.8：用MATLAB產(chǎn)生點(diǎn)數(shù)為32的均勻分布的隨機(jī)序列與高斯隨機(jī)序列clear all;N=32;x_rand=rand(1,N);x_randn=randn(1,N);xn=0:N-1;figure(1);stem(xn,x_rand)figure(2);stem(xn,x_randn)（三）、離散信號(hào)的基本運(yùn)算1、 信

9、號(hào)的延遲給定離散信號(hào)x(n)，若信號(hào)y(n)定義為：y(n)=x(n-k)，那么y(n)是信號(hào)x(n)在時(shí)間軸上右移k個(gè)抽樣周期得到的新序列。例3.1：正弦序列y(n)=sin(100n)右移3個(gè)抽樣周期后所得的序列，MATLAB程序如下：clear all;N=32;w=100;k=3;x1=zeros(1,k);xn=0:N-1;x2=sin(100*xn);figure(1)stem(xn,x2)x=x1 x2;axis(-1 N -1.1 1.1)N=N+k;xn=0:N-1;figure(2)stem(xn,x)axis(-1 N -1.1 1.1)利用for循環(huán)語句實(shí)現(xiàn)周期延遲.2

10、、 信號(hào)相加若信號(hào)，值得注意的是當(dāng)序列和的長度不相等或者位置不對(duì)應(yīng)時(shí)，首先應(yīng)該使兩者的位置對(duì)齊，然后通過zeros函數(shù)左右補(bǔ)零使其長度相等后再相加例3.2：用MATLAB實(shí)現(xiàn)兩序列相加clear all;n1=0:3x1=2 0.5 0.9 1;figure(1)stem(n1,x1)axis(-1 8 0 2.1 )n2=0:7x2= 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7;figure(2)stem(n2,x2)axis(-1 8 0 0.8 )n=0:7;x1=x1 zeros(1,8-length(n1);x2= zeros(1,8-length(n2),x2;x=

11、x1+x2;figure(3)stem(n,x)axis(-1 8 0 2.1)已知f1(t)=sinwt , f2(t)=sin8wt , w=2pi , 求f1(t)+f2(t)和f1(t)f2(t) 的波形圖3、 信號(hào)相乘信號(hào)序列和相乘所得信號(hào)的表達(dá)式為：這是樣本與樣本之間的點(diǎn)乘運(yùn)算，在MATLAB中可采用“.*”來實(shí)現(xiàn)，但是在信號(hào)序列相乘之前，應(yīng)對(duì)其做與相加運(yùn)算一樣的操作。例3.3：用MATLAB實(shí)現(xiàn)上例中兩序列相乘clear all;n1=0:3x1=2 0.5 0.9 1;figure(1)stem(n1,x1)axis(-1 8 0 2.1 )n2=0:7x2= 0 0.1 0.

12、2 0.3 0.4 0.5 0.6 0.7;figure(2)stem(n2,x2)axis(-1 8 0 0.8 )n=0:7;x1=x1 zeros(1,8-length(n1);x2= zeros(1,8-length(n2),x2;x=x1.*x2;figure(3)stem(n,x)axis(-1 8 0 0.35)4、 信號(hào)翻轉(zhuǎn)信號(hào)翻轉(zhuǎn)的表達(dá)式為：y(n)=x(-n)，在MATLAB中可以用fliplr函數(shù)實(shí)現(xiàn)此操作例3.4：用MATLAB實(shí)現(xiàn)“信號(hào)相加”中的序列翻轉(zhuǎn)clear all;n=0:3x1=2 0.5 0.9 1;x=fliplr(x1);stem(n,x)axis(-

13、1 4 0 2.1 )5、 信號(hào)和對(duì)于N點(diǎn)信號(hào)，其和的定義為：例3.5：用MATLAB實(shí)現(xiàn)“信號(hào)相加”中的序列和clear all;n=0:3x1=2 0.5 0.9 1;x=sum(x1)6、 信號(hào)積對(duì)于N點(diǎn)信號(hào)，其積的定義為：例3.5：用MATLAB實(shí)現(xiàn)“信號(hào)相加”中的序列積clear all;n=0:3x1=2 0.5 0.9 1;x=prod(x1)7 卷積1、 完成與兩函數(shù)的卷積運(yùn)算其中：在一個(gè)圖形窗口中，畫出、以及卷積結(jié)果。要求每個(gè)坐標(biāo)系有標(biāo)題、坐標(biāo)軸名稱。 p=0.1;t=0:p:10;f1=exp(-2*t).*u(t); f2=u(t)-u(t-4);f=conv(f1,f2

14、);subplot(1,3,1);plot(t,f1,r);title(f1(t)=e-2*t*u(t);xlabel(t(sec);ylabel(f1(t);subplot(1,3,2);plot(t,f2,g);title(f2(t)=u(t)-u(t-4);xlabel(t(sec);ylabel(f2(t);subplot(1,3,3);plot(f);title(f(t)=f1(t)*f2(t);xlabel(t(sec);ylabel(f(t);四 MATLAB在信號(hào)與系統(tǒng)中的應(yīng)用該實(shí)驗(yàn)用MATLAB中庫函數(shù)，如tf2zp(b,a)，ss2zp(A,B,C,D)，zplane(z,

15、p)，freqz(b,a)等。例如：1傳遞函數(shù)為，求其零極點(diǎn)圖。程序如下：num=1 0.5 2； 分子系數(shù)，按降冪順序排列den=1 0.4 1； 分母系數(shù)，按降冪順序排列z,p=tf2zp(num,den)； 用tf2zp函數(shù)求出其零點(diǎn)z和極點(diǎn)pzplane(z,p) 作出零極點(diǎn)圖 2若給出的是濾波器的輸入與輸出的狀態(tài)方程，如：，求其零極點(diǎn)圖。程序如下： A=1，0；1，-3； B=1；0； C=-，1； D=0； z,p=ss2zp(A,B,C,D)；求出其零極點(diǎn)z,p zplane(z,p)在連續(xù)時(shí)間系統(tǒng)中，當(dāng)極點(diǎn)在虛軸的右半平面時(shí)，系統(tǒng)不穩(wěn)定，在虛軸上，單階極點(diǎn)系統(tǒng)穩(wěn)定；若零點(diǎn)均處于

16、左半平面內(nèi)，則系統(tǒng)為最小相位系統(tǒng)。在離散系統(tǒng)中，極點(diǎn)在單位圓外系統(tǒng)不穩(wěn)定，在單位圓上，單階極點(diǎn)系統(tǒng)穩(wěn)定；零點(diǎn)在單位圓內(nèi)，系統(tǒng)為最小相位系統(tǒng)。對(duì)一濾波器，我們不僅要知道它的零點(diǎn)和極點(diǎn)，還要了解它的頻率特性，本實(shí)驗(yàn)可求其頻率特性。對(duì)模擬濾波器，可用freqs函數(shù)求得其頻率特性，對(duì)數(shù)字濾波器，則用freqz函數(shù)求得。3 已知模擬濾波器的傳遞函數(shù)為，求其頻率特性。 程序如下： num=0.2 0.3 1; den=1 0.4 1; w=logspace(-1,1); 頻率范圍 freqs(num,den,w)例4數(shù)字濾波器，取樣點(diǎn)數(shù)為128點(diǎn)，求其頻率特性。 程序如下： num=0.2 0.3 1;

17、den=1 0.4 1; freqz(num,den,128)實(shí)驗(yàn)內(nèi)容:1已知下列傳遞函數(shù)H(s)或H(z)，求其零極點(diǎn)，并畫出零極點(diǎn)圖。a. b.c.d. 2已知下列H(s)或H(z)，求其頻響。a. b. c. d. 1： 若某連續(xù)系統(tǒng)的輸入為e(t)，輸出為r(t)，系統(tǒng)的微分方程為：求該系統(tǒng)的單位沖激響應(yīng)h(t)及其單位階躍響應(yīng)g(t)。若 求出系統(tǒng)的零狀態(tài)響應(yīng)y(t)分析： 求沖激響應(yīng)及階躍響應(yīng)的MATLAB程序：a=1 5 6;b=3 2;subplot(2,1,1), impulse(b,a,4)subplot(2,1,2), step(b,a,4)運(yùn)行結(jié)果如右： 求零狀態(tài)響應(yīng)的

18、MATLAB程序：a=1 5 6;b=3 2;p1=0.01; %定義取樣時(shí)間間隔為0.01t1=0:p1:5; %定義時(shí)間范圍x1=exp(-2*t1); %定義輸入信號(hào)lsim(b,a,x1,t1), %對(duì)取樣間隔為0.01時(shí)系統(tǒng)響應(yīng)進(jìn)行仿真hold on; %保持圖形窗口以便能在同一窗口中繪制多條曲線p2=0.5; %定義取樣間隔為0.5t2=0:p2:5; %定義時(shí)間范圍x2=exp(-2*t2); %定義輸入信號(hào)lsim(b,a,x2,t2), hold off %對(duì)取樣間隔為0.5時(shí)系統(tǒng)響應(yīng)進(jìn)行仿真并解除保持運(yùn)行結(jié)果如下：2 已知一個(gè)過阻尼二階系統(tǒng)的狀態(tài)方程和輸出方程分別為：，r(

19、t)=0 1X(t)。 若系統(tǒng)初始狀態(tài)為X(0)=4 -5T， 求系統(tǒng)在作用下的全響應(yīng)。求全響應(yīng)程序如下：A=0 1 ; -2 -3 ;B=0 2;C=0 1;D=0;X0=4 -5; %定義系統(tǒng)初始狀態(tài)t=0: 0.01:10; E=3*exp(-4*t).*ones(size(t);%定義系統(tǒng)激勵(lì)信號(hào)r , x=lsim(A,B,C,D,E,t,X0);%求出系統(tǒng)全響應(yīng)的數(shù)值解plot(t,r) %繪制系統(tǒng)全響應(yīng)波形運(yùn)行結(jié)果如右。3 已知描述離散系統(tǒng)的差分方程為：，且已知系統(tǒng)輸入序列為， 求出系統(tǒng)的單位函數(shù)響應(yīng)h(k)在-3 10離散時(shí)間范圍內(nèi)響應(yīng)波形。 求出系統(tǒng)零狀態(tài)響應(yīng)在0 15區(qū)間上

20、的樣值；并畫出輸入序列的時(shí)域波形以及系統(tǒng)零狀態(tài)響應(yīng)的波形 分析：求系統(tǒng)的單位函數(shù)響應(yīng)的MATLAB程序： a=1,-0.25,0.5; b=1,1,0; impz(b,a,-3:10), title(單位響應(yīng)) %繪出單位函數(shù)響應(yīng)在-3 10區(qū)間上的波形 運(yùn)行結(jié)果如圖a。求零狀態(tài)響應(yīng)的MATLAB程序：a=1,-0.25,0.5;b=1,1,0k=0:15; %定義輸入序列取值范圍x=(1/2).k; %定義輸入序列表達(dá)式y(tǒng)=filter(b,a,x) %求解零狀態(tài)響應(yīng)樣值subplot(2,1,1),stem(k,x) %繪制輸入序列的波形 title(輸入序列)subplot(2,1,2)

21、,stem(k,y) %繪制零狀態(tài)響應(yīng)的波形title(輸出序列)運(yùn)行結(jié)果如下：y = Columns 1 through 10 1.0000 1.7500 0.6875 -0.3281 -0.2383 0.1982 0.2156 -0.0218 -0.1015 -0.0086 Columns 11 through 16 0.0515 0.0187 -0.0204 -0.0141 0.0069 0.0088圖a. 運(yùn)行結(jié)果 圖b. 運(yùn)行結(jié)果實(shí)驗(yàn)內(nèi)容1. 已知描述系統(tǒng)的微分方程和激勵(lì)信號(hào)e(t) 分別如下，并用MATLAB繪出系統(tǒng)單位沖激響應(yīng)和系統(tǒng)零狀態(tài)響應(yīng)的波形. ；如下圖所示的電路中，已知，

22、且兩電感上初始電流分別為，如果以電阻上電壓作為系統(tǒng)輸出，請(qǐng)求出系統(tǒng)在激勵(lì)（v）作用下的全響應(yīng)。2. 請(qǐng)用MATLAB分別求出下列差分方程所描述的離散系統(tǒng)，在020時(shí)間范圍內(nèi)的單位函數(shù)響應(yīng)、階躍響應(yīng)和系統(tǒng)零狀態(tài)響應(yīng)的數(shù)值解，并繪出其波形。另外，請(qǐng)將理論值與MATLAB仿真結(jié)果在對(duì)應(yīng)點(diǎn)上的值作比較，并說出兩者的區(qū)別和產(chǎn)生誤差的原因。 ； ； ；一帶通濾波器可由下列差分方程描述：，其中為系統(tǒng)輸入， 為系統(tǒng)輸出。請(qǐng)求出當(dāng)激勵(lì)為（選取適當(dāng)?shù)膎值）時(shí)濾波器的穩(wěn)態(tài)輸出。（四）、傅里葉變換的MATLAB實(shí)現(xiàn)在MATLAB中實(shí)現(xiàn)傅里葉變換的函數(shù)為：l F=fourier( f ) 對(duì)f(t)進(jìn)行傅里葉變換，其結(jié)

23、果為F(w)l F=fourier(f,v) 對(duì)f(t)進(jìn)行傅里葉變換，其結(jié)果為F(v)l F=fourier( f,u,v ) 對(duì)f(u)進(jìn)行傅里葉變換，其結(jié)果為F(v)傅里葉反變換l f=ifourier( F ) 對(duì)F(w)進(jìn)行傅里葉反變換，其結(jié)果為f(x)l f=ifourier(F,U) 對(duì)F(w)進(jìn)行傅里葉反變換，其結(jié)果為f(u)l f=ifourier( F,v,u ) 對(duì)F(v)進(jìn)行傅里葉反變換，其結(jié)果為f(u)例 求門函數(shù)的傅里葉變換，并畫出幅度頻譜圖MATLAB程序如下：syms t w %定義兩個(gè)符號(hào)變量t,wGt=sym(Heaviside(t+1)-Heaviside

24、(t-1); %產(chǎn)生門寬為2的門函數(shù)Fw=fourier(Gt,t,w); %對(duì)門函數(shù)作傅氏變換求F(jw)FFw=maple(convert,Fw,piecewise); %數(shù)據(jù)類型轉(zhuǎn)換，轉(zhuǎn)為分段函數(shù)，此處可以去掉FFP=abs(FFw); %求振幅頻譜| F(jw)|ezplot(FFP,-10*pi 10*pi);grid; %繪制函數(shù)圖形，并加網(wǎng)格axis(-10*pi 10*pi 0 2.2) %限定坐標(biāo)軸范圍運(yùn)行結(jié)果：Fw= exp(i*w)*(pi*Dirac(w)-i/w)-exp(-i*w)*(pi*Dirac(w)-i/w) % Dirac(w)為()，即傅立葉變換結(jié)果中含

25、有奇異函數(shù)，故繪圖前需作函數(shù)類型轉(zhuǎn)換FFw= -i*exp(i*w)/w+i*exp(-i*w)/w % FFw為復(fù)數(shù) FFP= abs(-i*exp(i*w)/w+i*exp(-i*w)/w) %求FFw的模值例 求函數(shù)的傅里葉反變換f(t) MATLAB程序如下：syms t w %定義兩個(gè)符號(hào)變量t,wFw=sym(1/(1+w2); %定義頻譜函數(shù)F(jw)ft=ifourier(Fw,w,t); %對(duì)頻譜函數(shù)F(jw)進(jìn)行傅氏反變換運(yùn)行結(jié)果：ft = 1/2*exp(-t)*Heaviside(t)+1/2*exp(t)*Heaviside(-t)2、傅里葉變換的數(shù)值計(jì)算實(shí)現(xiàn)法嚴(yán)格說

26、來，如果不使用symbolic工具箱，是不能分析連續(xù)時(shí)間信號(hào)的。采用數(shù)值計(jì)算方法實(shí)現(xiàn)連續(xù)時(shí)間信號(hào)的傅里葉變換，實(shí)質(zhì)上只是借助于MATLAB的強(qiáng)大數(shù)值計(jì)算功能，特別是其強(qiáng)大的矩陣運(yùn)算能力而進(jìn)行的一種近似計(jì)算。傅里葉變換的數(shù)值計(jì)算實(shí)現(xiàn)法的原理如下：對(duì)于連續(xù)時(shí)間信號(hào)f(t)，其傅里葉變換為：其中為取樣間隔，如果f(t)是時(shí)限信號(hào)，或者當(dāng)|t|大于某個(gè)給定值時(shí)，f(t)的值已經(jīng)衰減得很厲害，可以近似地看成是時(shí)限信號(hào)，則上式中的n取值就是有限的，假定為N，有： 若對(duì)頻率變量進(jìn)行取樣，得：通常?。海渲惺且〉念l率范圍，或信號(hào)的頻帶寬度。采用MATLAB實(shí)現(xiàn)上式時(shí)，其要點(diǎn)是要生成f(t)的N個(gè)樣本值的向量

27、，以及向量，兩向量的內(nèi)積（即兩矩陣的乘積），結(jié)果即完成上式的傅里葉變換的數(shù)值計(jì)算。注意：時(shí)間取樣間隔的確定，其依據(jù)是必須小于奈奎斯特（Nyquist）取樣間隔。如果f(t)不是嚴(yán)格的帶限信號(hào)，則可以根據(jù)實(shí)際計(jì)算的精度要求來確定一個(gè)適當(dāng)?shù)念l率為信號(hào)的帶寬。例 用數(shù)值計(jì)算法實(shí)現(xiàn)上面門函數(shù)的傅里葉變換，并畫出幅度頻譜圖. 分析： 該信號(hào)的頻譜為，其第一個(gè)過零點(diǎn)頻率為，一般將此頻率認(rèn)為是信號(hào)的帶寬。但考慮到的形狀（為抽樣函數(shù)），假如將精度提高到該值的50倍，即取，則據(jù)此確定的Nyquist取樣間隔為：。MATLAB程序如下：R=0.02; %取樣間隔=0.02t=-2:R:2; % t為從-2到2，間

28、隔為0.02的行向量,有201個(gè)樣本點(diǎn)ft=zeros(1,50),ones(1,101),zeros(1,50);% 產(chǎn)生f(t)的樣值矩陣（即f(t)的樣本值組成的行向量）W1=10*pi; %取要計(jì)算的頻率范圍M=500; k=0:M; w=k*W1/M; %頻域采樣數(shù)為M, w為頻率正半軸的采樣點(diǎn)Fw=ft*exp(-j*t*w)*R; %求傅氏變換F(jw) FRw=abs(Fw); %取振幅W=-fliplr(w),w(2:501) ; %由信號(hào)雙邊頻譜的偶對(duì)稱性，利用fliplr(w)形成負(fù)半軸的點(diǎn)，% w(2:501)為正半軸的點(diǎn)，函數(shù)fliplr(w)對(duì)矩陣w行向量作180度

29、反轉(zhuǎn)FW=fliplr(FRw),FRw(2:501); %形成對(duì)應(yīng)于2M+1個(gè)頻率點(diǎn)的值Subplot(2,1,1) ; plot(t,ft) ;grid; %畫出原時(shí)間函數(shù)f(t)的波形，并加網(wǎng)格xlabel(t) ; ylabel(f(t); %坐標(biāo)軸標(biāo)注title(f(t)=u(t+1)-u(t-1); %文本標(biāo)注subplot(2,1,2) ; plot(W,FW) ;grid on;%畫出振幅頻譜的波形,并加網(wǎng)格xlabel (W) ; ylabel (F(W); %坐標(biāo)軸標(biāo)注title(f(t)的振幅頻譜圖); %文本標(biāo)注運(yùn)行結(jié)果如下：三、 實(shí)驗(yàn)內(nèi)容1.編程實(shí)現(xiàn)求下列信號(hào)的幅度頻

30、譜(1) 求出的頻譜函數(shù)F1(j)，請(qǐng)將它與上面門寬為2的門函數(shù)的頻譜進(jìn)行比較，觀察兩者的特點(diǎn)，說明兩者的關(guān)系。(2) 三角脈沖 (3) 單邊指數(shù)信號(hào) (4) 高斯信號(hào)2.利用ifourier( ) 函數(shù)求下列頻譜函數(shù)的傅氏反變換(1) (2) 例：若是一個(gè)N=32的有限序列，利用MATLAB計(jì)算它的DFT并畫出圖形。N=32;n=0:N-1;xn=cos(pi*n/6);k=0:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;figure(1)stem(n,xn)figure(2)stem(k,abs(Xk)在MATLAB中，可以直接利用

31、內(nèi)部函數(shù)fft來實(shí)現(xiàn)FFT算法，該函數(shù)是機(jī)器語言，而不是MATLAB指令寫成的，執(zhí)行速度很快。常用格式為：y=fft(x)y=fft(x,N)42 是一個(gè)N16的有限序列，用MATLAB求其DFT的結(jié)果，并畫出其結(jié)果圖，如圖31所示。圖 31 有限長序列的DFT結(jié)果圖程序N=16;n=0:1:N-1; %時(shí)域采樣xn=sin(n*pi/8)+sin(n*pi/4);k=0:1:N-1; %頻域采樣WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;subplot(2,1,1)stem(n,xn);subplot(2,1,2)stem(k,abs(Xk)

32、;，n=0,10運(yùn)算結(jié)果Xk = Columns 1 through 5 0.0000 -0.0000 - 8.0000i -0.0000 - 8.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i Columns 6 through 10 -0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i Columns 11 through 15 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000

33、 - 0.0000i 0.0000 + 8.0000i Column 16 0.0000 + 8.0000i用MATLAB計(jì)算序列-2 0 1 1 3和序列1 2 0 -1的離散卷積。 解 MATLAB程序如下： a=-2 0 1 -1 3; b=1 2 0 -1; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel(n); ylabel(幅度); 圖1.1給出了卷積結(jié)果的圖形，求得的結(jié)果存放在數(shù)組c中為：-2 -4 1 3 1 5 1 -3。 1-2 用MATLAB計(jì)算差分方程 當(dāng)輸入序列為 時(shí)的輸出結(jié)果 。 解 MATLAB程序如下

34、： N=41;a=0.8 -0.44 0.36 0.22;b=1 0.7 -0.45 -0.6;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y)xlabel(n);ylabel(幅度) 圖 1.2 給出了該差分方程的前41個(gè)樣點(diǎn)的輸出，即該系統(tǒng)的單位脈沖響應(yīng)。 用MATLAB計(jì)算例1-2差分方程 所對(duì)應(yīng)的系統(tǒng)函數(shù)的DTFT。 解 例1-2差分方程所對(duì)應(yīng)的系統(tǒng)函數(shù)為： 其DTFT為 用MATLAB計(jì)算的程序如下：k=256;num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=

35、freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);gridtitle(實(shí)部)xlabel(omega/pi);ylabel(幅度)subplot(2,2,2);plot(w/pi,imag(h);gridtitle(虛部)xlabel(omega/pi);ylabel(Amplitude)subplot(2,2,3);plot(w/pi,abs(h);gridtitle(幅度譜)xlabel(omega/pi);ylabel(幅值)subplot(2,2,4);plot(w/pi,angle(h);gridtitle(相位譜)xlabel(om

36、ega/pi);ylabel(弧度) 2-1 對(duì)連續(xù)的單一頻率周期信號(hào) 按采樣頻率 采樣，截取長度N分別選N =20和N =16，觀察其DFT結(jié)果的幅度譜。 解 此時(shí)離散序列 ，即k=8。用MATLAB計(jì)算并作圖，函數(shù)fft用于計(jì)算離散傅里葉變換DFT，程序如下： k=8; n1=0:1:19; xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1) xlabel(t/T);ylabel(x(n); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1) xlabel(k);ylabel(X(k);

37、n2=0:1:15; xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2) xlabel(t/T);ylabel(x(n); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2) xlabel(k);ylabel(X(k); 計(jì)算結(jié)果示于圖2.1，(a)和(b)分別是N=20時(shí)的截取信號(hào)和DFT結(jié)果，由于截取了兩個(gè)半周期，頻譜出現(xiàn)泄漏；(c) 和(d) 分別是N=16時(shí)的截取信號(hào)和DFT結(jié)果，由于截取了兩個(gè)整周期，得到單一譜線的頻譜。上述頻譜的誤差主要是由于時(shí)域中對(duì)信號(hào)的非整周期截?cái)喈a(chǎn)生的頻譜泄漏

38、。 五、快速傅立葉變換當(dāng)序列長度N很大時(shí)，直接計(jì)算序列的離散傅立葉變換的計(jì)算量非常大，因此引入了快速傅立葉變換(FFT),所有這些高效算法合起來統(tǒng)稱為快速傅立葉變換算法。Matlab中自帶有FFT函數(shù)，試用help函數(shù)察看其具體使用方法。已知帶有測量噪聲的信號(hào)其中f150Hz，f2120Hz，為均值為0的隨機(jī)信號(hào)，數(shù)據(jù)采樣頻率為1000Hz，數(shù)據(jù)點(diǎn)數(shù)N=1024，繪制信號(hào)的頻譜圖和無噪聲信號(hào)的頻譜圖。fs=1000;N=1024;n=0:N-1;t=n/fs;f1=50;f2=120;%信號(hào)沒有噪聲情況x=sin(2*pi*f1*t)+sin(2*pi*f2*t);y=fft(x,N);mag

39、=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:N/2),mag(1:N/2);xlabel(Frequence(Hz);ylabel(Magnitude);title(N=1024 沒有噪聲);grid%信號(hào)有噪聲情況x=x+2*randn(1,length(t);y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*fs/length(y);subplot(211);plot(f(1:N/2),mag(1:N/2);xlabel(Frequence(Hz);ylabel(Magnitude)

40、;title(N=1024 有噪聲);grid1 復(fù)指數(shù)信號(hào)的離散傅里葉變換。其中 ，n=0,10用MATLAB求這一有限時(shí)寬的序列的傅里葉變換。4.3 z正反變換序列的z變換定義為 （3-1）其中，符號(hào)表示取z變換，z是復(fù)變量。相應(yīng)地，單邊z變換定義為 （3-2）MATLAB符號(hào)數(shù)學(xué)工具箱提供了計(jì)算離散時(shí)間信號(hào)單邊z變換的函數(shù)ztrans和z反變換函數(shù)iztrans，其語句格式分別為Z=ztrans(x)x=iztrans(z)上式中的x和Z分別為時(shí)域表達(dá)式和z域表達(dá)式的符號(hào)表示，可通過sym函數(shù)來定義?！緦?shí)例4-4】 試用ztrans函數(shù)求下列函數(shù)的z變換。（1）； （2）。解：（1）z變

41、換MATLAB源程序?yàn)閤=sym(an*cos(pi*n);Z=ztrans(x);simplify(Z)ans=z/(z+a)（2）z變換MATLAB源程序?yàn)閤=sym(2(n-1)-(-2)(n-1);Z=ztrans(x);simplify(Z)ans=z2/(z-2)/(z+2)【實(shí)例4-5】 試用iztrans函數(shù)求下列函數(shù)的z反變換。（1） （2）解：（1）z反變換MATLAB源程序?yàn)閆=sym(8*z-19)/(z2-5*z+6);x=iztrans(Z);simplify(x)ans=-19/6*charfcn0(n)+5*3(n-1)+3*2(n-1)其中，charfcn0(

42、n)是函數(shù)在MATLAB符號(hào)工具箱中的表示，反變換后的函數(shù)形式為。（2）z反變換MATLAB源程序?yàn)閆=sym(z*(2*z2-11*z+12)/(z-1)/(z-2)3);x=iztrans(Z);simplify(x)ans=-3+3*2n-1/4*2n*n-1/4*2n*n2 -3+3*2n-1/4*2n*n-1/4*2n*n2其函數(shù)形式為。如果信號(hào)的z域表示式是有理函數(shù)，進(jìn)行z反變換的另一個(gè)方法是對(duì)進(jìn)行部分分式展開，然后求各簡單分式的z反變換。設(shè)的有理分式表示為 （3-3）MATLAB信號(hào)處理工具箱提供了一個(gè)對(duì)進(jìn)行部分分式展開的函數(shù)residuez，其語句格式為R,P,K=residu

43、ez(B,A)其中，B，A分別表示X(z)的分子與分母多項(xiàng)式的系數(shù)向量；R為部分分式的系數(shù)向量；P為極點(diǎn)向量；K為多項(xiàng)式的系數(shù)。若X(z)為有理真分式，則K為零?！緦?shí)例4-6】 試用MATLAB命令對(duì)函數(shù)進(jìn)行部分分式展開，并求出其z反變換。解：MATLAB源程序?yàn)?B=18; A=18,3,-4,-1; R,P,K=residuez(B,A)R= 0.3600 0.2400 0.4000P= 0.5000 -0.3333 -0.3333K= 從運(yùn)行結(jié)果可知，表示系統(tǒng)有一個(gè)二重極點(diǎn)。所以，X(z)的部分分式展開為因此，其z反變換為（五）、IIR數(shù)字濾波器設(shè)計(jì)1、 基于巴特沃斯法直接設(shè)計(jì)IIR數(shù)字

44、濾波器例.1：設(shè)計(jì)一個(gè)10階的帶通巴特沃斯數(shù)字濾波器，帶通頻率為100Hz到200Hz，采樣頻率為1000Hz，繪出該濾波器的幅頻于相頻特性，以及其沖擊響應(yīng)圖clear all;N=10;Wn=100 200/500;b,a=butter(N,Wn,bandpass);freqz(b,a,128,1000)figure(2)y,t=impz(b,a,101);stem(t,y)2、 基于切比雪夫法直接設(shè)計(jì)IIR數(shù)字濾波器例5.2：設(shè)計(jì)一個(gè)切比雪夫型數(shù)字低通濾波器，要求：Ws=200Hz,Wp=100Hz,Rp=3dB,Rs=30dB,Fs=1000Hzclear all;Wp=100;Rp=3;Ws=200;Rs=30;Fs=1000;N,Wn=cheb1ord(Wp/(Fs/2),Ws/(Fs/2),Rp,Rs);b,a=cheby1(N,Rp,Wn);freqz(b,a,512,1000