信號與系統(tǒng) 第二版 重大 例題程序 第8章

例8.1用Matlab命令繪制連續(xù)時間信號。Matlab程序如下：t=0:0.01:10;%設(shè)定時間向量f=5*exp(-0.5*t).*cos(2*t);%產(chǎn)生信號f(t)plot(t,f)%繪制f(t)xlabel('t'),ylabel('f(t)')%設(shè)定坐標(biāo)軸名程序運(yùn)行結(jié)果如圖：圖8.2.1例8.1信號波形例8.2用Matlab命令繪制沖擊函數(shù)信號。Matlab程序如下：t=-5:0.01:5;%設(shè)定時間向量y=0*(t>=-5&t<0)+1*(t==1)+0*(t>0&t<=5);%根據(jù)沖擊函數(shù)性質(zhì)給出沖擊函數(shù)表達(dá)式plot(t,y)%繪制沖擊函數(shù)axis([-5,5,-1,1.5]);%設(shè)置坐標(biāo)范圍程序運(yùn)行結(jié)果如圖：圖8.2.2例8.2信號波形例8.3用Matlab命令繪制階躍函數(shù)信號。在Matlab中，產(chǎn)生階躍函數(shù)有多種方法：用符號函數(shù)實(shí)現(xiàn)；用stepfun()實(shí)現(xiàn)；用比較表達(dá)式實(shí)現(xiàn)。（1）用符號函數(shù)實(shí)現(xiàn)的Matlab程序如下：t=-200:200;%設(shè)定時間向量f1=sign(t);%產(chǎn)生符號函數(shù)f2=1/2+f1/2;%產(chǎn)生單位階躍信號subplot(1,2,1),plot(t,f1)%繪制符號函數(shù)axis([-200200-1.11.1])%設(shè)置坐標(biāo)范圍xlabel('t'),ylabel('sgn(t)')%設(shè)置坐標(biāo)軸名title('符號函數(shù)')%設(shè)置圖名subplot(1,2,2),plot(t,f2)%繪制單位階躍信號axis([-200200-0.11.1])%設(shè)置坐標(biāo)范圍title('單位階躍信號')%設(shè)置圖名xlabel('t'),ylabel('u(t)')%設(shè)置坐標(biāo)軸名程序運(yùn)行結(jié)果如圖：圖8.2.3例8.3用符號函數(shù)實(shí)現(xiàn)的信號波形（2）stepfun（T,T0）：T為時間向量，T0為階躍時刻，當(dāng)T<T0時返回0，當(dāng)T>T0時返回1，且返回的向量與T具有相同的長度。用stepfun()實(shí)現(xiàn)的Matlab程序如下：t=linspace(-200,200,40000);%設(shè)定時間向量y=stepfun(t,1);%設(shè)定1為階躍時刻plot(t,y);%繪制階躍函數(shù)xlabel('t'),ylabel('x'),%設(shè)置坐標(biāo)軸名axis([-200200-0.11.1])%設(shè)置坐標(biāo)范圍程序運(yùn)行結(jié)果如圖：圖8.2.4例8.3用stepfun()實(shí)現(xiàn)實(shí)現(xiàn)的信號波形（3）用比較表達(dá)式實(shí)現(xiàn)的Matlab程序如下：t=linspace(-200,200,40000);%設(shè)定時間向量x=t>0;%設(shè)置比較時刻plot(t,x)%繪制階躍函數(shù)xlabel('t'),ylabel('x'),%設(shè)置坐標(biāo)軸名axis([-200200-0.11.1])%設(shè)置坐標(biāo)范圍程序運(yùn)行結(jié)果如圖：圖8.2.5例8.3用比較表達(dá)式實(shí)現(xiàn)實(shí)現(xiàn)的信號波形例8.4用Matlab求解例2.6中的LTI系統(tǒng)的全響應(yīng)。解：先求零輸入響應(yīng)，再求零狀態(tài)響應(yīng)，兩部分加起來就是系統(tǒng)的全響應(yīng)。1）Matlab程序如下：eq1='D2y+5*Dy+6*y=0';%設(shè)定零輸入條件下的微分方程ic1='y(0)=2,Dy(0)=-1';%輸入初始狀態(tài)yzi=dsolve(eq1,ic1);%求解微分方程，得到零輸入響應(yīng)yzi=simplify(yzi)%化簡零輸入響應(yīng)eq2='D2y+5*Dy+6*y=2*exp(-t)';%設(shè)定給定輸入條件下的微分方程ic2='y(0)=0,Dy(0)=0';%設(shè)定初始狀態(tài)為0yzs=dsolve(eq2,ic2);%求解微分方程，得到零狀態(tài)響應(yīng)yzs=simplify(yzs)%化簡零狀態(tài)響應(yīng)y=yzi+yzs%全響應(yīng)運(yùn)行程序，得到：yzi=exp(-3*t)*(5*exp(t)-3)yzs=exp(-3*t)*(exp(t)-1)^2y=exp(-3*t)*(5*exp(t)-3)+exp(-3*t)*(exp(t)-1)^2全響應(yīng)結(jié)果經(jīng)化解整理后與例2.6第1）求得的全解一致。2）Matlab程序如下：eq1='D2y+5*Dy+6*y=0';%設(shè)定零輸入條件下的微分方程ic1='y(0)=1,Dy(0)=0';%輸入初始狀態(tài)yzi=dsolve(eq1,ic1);%求解微分方程，得到零輸入響應(yīng)yzi=simplify(yzi)%化簡零輸入響應(yīng)eq2='D2y+5*Dy+6*y=exp(-2*t)';%設(shè)定給定輸入條件下的微分方程ic2='y(0)=0,Dy(0)=0';%設(shè)定初始狀態(tài)為0yzs=dsolve(eq2,ic2);%求解微分方程，得到零狀態(tài)響應(yīng)yzs=simplify(yzs)%化簡零狀態(tài)響應(yīng)y=yzi+yzs%全響應(yīng)運(yùn)行程序，得到：yzi=exp(-3*t)*(3*exp(t)-2)yzs=exp(-3*t)*(t*exp(t)-exp(t)+1)y=exp(-3*t)*(3*exp(t)-2)+exp(-3*t)*(t*exp(t)-exp(t)+1)全響應(yīng)結(jié)果經(jīng)化解整理后與例2.6第2）求得的全解一致。例8.5用Matlab求解第2章例2.9中的LTI系統(tǒng)的零輸入響應(yīng)、零狀態(tài)響應(yīng)和全響應(yīng)。解：此處采用求解零狀態(tài)響應(yīng)函數(shù)initial()、系統(tǒng)的時域響應(yīng)函數(shù)lsim()進(jìn)行編程作圖求解。Matlab程序如下：a=[132];%設(shè)定系數(shù)向量ab=[26];%設(shè)定系數(shù)向量bsys=tf(b,a);%得到系統(tǒng)的傳遞函數(shù)模型[A,B,C,D]=tf2ss(b,a);%由傳遞函數(shù)模型求得狀態(tài)空間模型sys1=ss(A,B,C,D);%狀態(tài)空間模型y0=[20];%輸入系統(tǒng)的初始條件t=0:0.01:5;f=heaviside(t);%系統(tǒng)輸入y1=initial(sys1,y0,t);%零輸入響應(yīng)subplot(2,2,1);plot(t,y1);%繪制多子圖xlabel('t');ylabel('y1')title('零輸入響應(yīng)')y2=lsim(sys1,f,t);%零狀態(tài)響應(yīng)subplot(2,2,2);plot(t,y2);xlabel('t');ylabel('y2')title('零狀態(tài)響應(yīng)')y3=y1+y2;%全響應(yīng)subplot(2,1,2);plot(t,y3);xlabel('t');ylabel('y3')title('全響應(yīng)')程序運(yùn)行結(jié)果如下圖：圖8.2.6例8.5各種響應(yīng)曲線例8.6用Matlab求解第3章例3.1中的周期為T，脈沖幅度為1的矩形波信號的傅里葉級數(shù)展開。解：對矩形波信號進(jìn)行傅里葉級數(shù)展開最關(guān)鍵的是求得傅里葉級數(shù)的系數(shù)、，Matlab編程需要制定周期為固定值，故這里假設(shè)周期為4，脈沖寬度為2。程序如下：T=4;%設(shè)定矩形脈沖信號的周期w=2*pi/T;%基頻F=@(t)abs(t)<=0;a0=quad(F,-2,2)/T;%求a0N=50;%設(shè)定最大諧波數(shù)an=zeros(1,N);bn=zeros(1,N);%初始化和fork=1:Nan(k)=quad(@rectcos,-2,2,[],[],k,w)*2/T;%計算anbn(k)=quad(@rectsin,-2,2,[],[],k,w)*2/T;%計算bnend%繪制周期矩形脈沖信號t=-4:0.01:4;%設(shè)定時間向量x=pulstran(t,-6:4:6,'rectpuls',2);%生成周期矩形脈沖信號subplot(3,2,1);plot(t,x);%繪制周期矩形脈沖信號xlabel('t');ylabel('Amplitude');axis([-4,4,-0.2,1.2]);gridon;title('Rectangularwave');%有限項(xiàng)級數(shù)逼近k_max=[110304050];%最大諧波數(shù)向量form=1:length(k_max)f=a0;fork=1:k_max(m)f=f+an(k)*cos(k*w*t)+bn(k)*sin(k*w*t);endsubplot(3,2,m+1);plot(t,f);gridon;xlabel('t');ylabel('PartialSum');axis([-4,4,-1,1.2]);title(['MaxHarmonics=',num2str(k_max(m))]);end仿真結(jié)果如下圖：圖8.2.7例8.6有限項(xiàng)傅里葉級數(shù)和例8.7用Matlab求解第3章例3.17中均勻沖擊序列的傅里葉變換。解：從圖3.6.1中可得均勻沖激序列的波形圖和頻譜圖。根據(jù)歐拉公式，可將展開成三角函數(shù)形式的傅里葉級數(shù)。級數(shù)中的每一項(xiàng)都是連續(xù)函數(shù)，這里將用傅里葉級數(shù)的有限項(xiàng)來近似單位沖激序列，觀察有限項(xiàng)級數(shù)信號疊加后逐漸演變成離散沖激的過程。Matlab實(shí)現(xiàn)代碼如下：T=1;%設(shè)周期為1w1=2*pi/T;t=-1.2*T:0.01:1.2*T;%設(shè)定時間起始值的間隔N=5;%取傅里葉級數(shù)的前5項(xiàng)近似單位沖激序列xN=ones(1,length(t))/T;%傅里葉級數(shù)首項(xiàng)表達(dá)式plot(t,xN,'m');%作首項(xiàng)圖xlabel('t');%標(biāo)識X軸ylabel('x_N(t)');%標(biāo)識Y軸holdon%疊繪color=['m','r','g','b'];%設(shè)定曲線顏色集forn=1:N-1xn=cos(n*w1*t)/T*2;xN=xn+xN;%求傅里葉級數(shù)的后續(xù)項(xiàng)和plot(t,xN,color(mod(n,length(color))+1));%作前n項(xiàng)和圖endh=plot(t,xN,'k');%作前5項(xiàng)和圖set(h,'linewidth',3);%設(shè)定前5項(xiàng)和的曲線線寬為3axis([-1.2*T,1.2*T,min(xN)*1.1,max(xN)*1.1]);%設(shè)定坐標(biāo)軸范圍(a)N=5(b)N=20圖8.2.8例8.7有限項(xiàng)傅里葉級數(shù)合成單位沖激序列例8.8求的傅里葉變化，并繪制幅度譜。解：傅里葉變換指令較簡單，值得注意的一點(diǎn)的需要先定義符號變量，階躍函數(shù)在Matlab中的表示為Heaviside。Matlab實(shí)現(xiàn)傅里葉變換代碼如下：symst;%定義符號變量tf=sym('heaviside(t+1)')-sym('heaviside(t-1)');%輸入符號函數(shù)fF=fourier(f);%求f的傅立葉變換Fezplot(abs(F));%繪制幅度譜圖8.2.9例8.8幅度譜例8.9用Matlab求解第3章例3.21中LTI系統(tǒng)的頻率響應(yīng)。解：本題編輯程序的重點(diǎn)在于準(zhǔn)確給出LTI系統(tǒng)的常微分方程的系數(shù)，然后使用freqs()指令即可通過圖形展示該系統(tǒng)的幅頻特性和相頻特性。程序如下：b=[143];a=[12];%設(shè)定分子、分母多項(xiàng)式系數(shù)向量w=-20:0.05:20;%設(shè)定頻率向量H=freqs(b,a,w);%求頻率響應(yīng)mag=abs(H);%取頻率響應(yīng)的模phase=angle(H)*180/pi;%取頻率響應(yīng)的相位，并轉(zhuǎn)換成角度subplot(2,1,1),plot(w,mag)%以線性坐標(biāo)繪制幅頻特性曲線xlabel('Frequency(rad/s)'),ylabel('Magnitude'),gridsubplot(2,1,2),plot(w,phase)%以線性坐標(biāo)繪制相頻特性曲線xlabel('Frequency(rad/s)'),ylabel('Phase(degrees)'),grid圖8.2.10例8.9頻率響應(yīng)曲線例8.10用Matlab求例4.4、例4.5信號的單邊拉普拉斯變換。解：例4.4Matlab程序如下：symsta;%定義符號變量t，af=exp(-a*t).*heaviside(t);%輸入符號函數(shù)F=laplace(f);%求f的拉普拉斯變換FF運(yùn)行結(jié)果為:F=1/(a+s)例4.5Matlab程序如下：symsta;%定義符號變量t，af=1+exp(-a*t).*heaviside(t);%輸入符號函數(shù)fF=laplace(f);%求f的拉普拉斯變換FF運(yùn)行結(jié)果為：F=1/(a+s)+1/s注：對比Matlab運(yùn)行結(jié)果與例4.4的計算結(jié)果，Matlab沒有對收斂域進(jìn)行區(qū)分，只是給出了積分收斂時的結(jié)果。例4.5的計算結(jié)果與Matlab運(yùn)行結(jié)果一致。另外，這里沒有給出a的具體值，所以不能使用ezplot()繪圖。例8.11用Matlab求例4.24象函數(shù)的原函數(shù)。解：求原函數(shù)即求拉普拉斯逆變換，采用ilaplace()函數(shù)實(shí)現(xiàn)。Matlab程序如下：symss;%定義符號變量sL=(s+6)/(s*(s+1)*(s+2));%輸入符號函數(shù)Lf=ilaplace(L);%求L的拉普拉斯逆變換ffezplot(f)%原函數(shù)作圖運(yùn)行結(jié)果為：f=2*exp(-2*t)-5*exp(-t)+3圖8.2.11例8.11原函數(shù)曲線注：對比例4.24的計算結(jié)果與Matlab的仿真結(jié)果，Matlab仿真結(jié)果中沒有階躍函數(shù)，因?yàn)閘aplace()和ilaplace()函數(shù)均有單邊變換，默認(rèn)取值大于0。另外，由于該題原函數(shù)中各參數(shù)確定，可以用ezplot()函數(shù)作圖直觀表示原函數(shù)。例8.12用Matlab求例4.29LTI系統(tǒng)的零輸入響應(yīng)、零狀態(tài)響應(yīng)和全響應(yīng)。解：由例4.29可知，對系統(tǒng)的微分方程兩邊取拉普拉斯變換可得：整理得并代入初始狀態(tài)得：則可以用Matlab求其全響應(yīng)、零輸入響應(yīng)、零狀態(tài)響應(yīng)。程序如下：symstsYzs=(s+4)/(s^2+5*s+6);yzs=ilaplace(Yzs)ft=heaviside(t);F=laplace(ft);Yzf=F*2*(s+3)/(s^2+5*s+6);yzf=ilaplace(Yzf)yt=simplify(yzs+yzf)運(yùn)行結(jié)果為：yzs=2*exp(-2*t)-exp(-3*t)yzf=1-exp(-2*t)yt=exp(-2*t)-exp(-3*t)+1例8.13用Matlab求例5.5離散系統(tǒng)差分方程，，，激勵為，求系統(tǒng)的完全響應(yīng)。程序如下：a=[6,-5,1];b=1;n=0:20;x=10*cos(0.5*pi*n);%生成激勵信號x（n）wi=[-10-a(2)*(1/2)-a(1)*(-5/4),0];%定義初始狀態(tài)，wi=f(2)-a(2)*y(1)-a(1)*y(2),從題目中給出%的初始條件可以得知f(2)=-10，y(2)=-5/4[y,wf]=filter(b,a,x,wi);%濾波得到完全響應(yīng)plot(n,y);xlabel('n');ylabel('y')運(yùn)行結(jié)果如下：圖8.3.5例8.13系統(tǒng)完全響應(yīng)圖例8.14用Matlab求例5.8離散系統(tǒng)差分方程的單位序列響應(yīng)，并與例5.8理論值做對比。程序如下：

k=0:10;a=[1,-3,3,-1];b=[1];h=impz(b,a,k);subplot(2,1,1);stem(k,h);title('單位沖擊響應(yīng)近似值');gridon;hk=1+3/2.*k+1/2.*k.^2;%在時為1subplot(2,1,2);stem(k,hk);title('單位沖擊響應(yīng)理論值');gridon;運(yùn)行結(jié)果如下：圖8.3.6例8.14離散時間系統(tǒng)單位序列響應(yīng)對比兩個圖可知，用求得的單位序列響應(yīng)近似值與理論值基本相符。例8.15用Matlab求例5.10離散系統(tǒng)差分方程的單位序列響應(yīng)和單位階躍響應(yīng)，并與理論值、做對比。程序如下：k=0:10;a=[1,-5/6,1/6];b=[1];h=impz(b,a,k);g=stepz(b,a,k);subplot(2,2,1)stem(k,h);title('單位序列響應(yīng)近似值');gridon;subplot(2,2,2)stem(k,g);title('單位階躍響應(yīng)近似值');gridon;hk=3*(1/2).^k-2*(1/3).^k;subplot(2,2,3)stem(k,hk);title('單位序列響應(yīng)理論值')gk=3-3*(1/2).^k+(1/3).^k;subplot(2,2,4)stem(k,gk);title('單位階躍響應(yīng)理論值')圖8.3.7例8.15離散時間系統(tǒng)單位序列響應(yīng)及單位階躍響應(yīng)對比上圖可知，用求得的單位序列響應(yīng)近似值與理論值基本相符，用求得的單位階躍響應(yīng)近似值與理論值基本相符。例8.16用Matlab求例5.11中，的卷積和。程序如下：n=[0:20];h1=exp(-n);h2=(n>=0);y=conv(single(h1),single(h2));%conv指令要求兩個向量是一種精度，故需要進(jìn)行轉(zhuǎn)換；卷積和%的長度length(y)=length(h1)+length(h2)-1。n=0:2*length(n)-2;%stem繪圖要求兩個參數(shù)長度一致，故進(jìn)行n的轉(zhuǎn)換stem(n,y);xlabel('n');ylabel('y');圖8.3.8例8.16卷積和圖形例8.17用Matlab求例6.1中單位序列、例6.3中因果序列的z變換。（1）在Matlab中，kroneckerDelta(n,0)代表單位序列函數(shù)，參考程序如下：f=sym('kroneckerDelta(n,0)');F=ztrans(f)運(yùn)行結(jié)果為：F=1程序運(yùn)行結(jié)果與例6.1理論計算結(jié)果一致。（2）參考程序如下：f=sym('a^n');F=ztrans(f)運(yùn)行結(jié)果為：F=-z/(a-z)程序運(yùn)行結(jié)果與例6.3理論計算結(jié)果一致。例8.18用Matlab求例6.20中象函數(shù)的原序列。程序如下：F=sym('(z^2+z)/(z-1)^2');f=iztrans(F)運(yùn)行結(jié)果如下：f=2*n+1注：利用iztrans指令求得的是收斂域時的原序列，需要在程序運(yùn)行結(jié)果后加單位階躍函數(shù)，結(jié)果為f=（2*n+1），與理論結(jié)果一致。例8.19用Matlab求例6.26中二階離散系統(tǒng)差分方程，，的零輸入響應(yīng)、零狀態(tài)響應(yīng)和全響應(yīng)。解：首先對差分方程左右兩邊取單邊z變換，得整理得：接下來可以利用Matlab求其完全響應(yīng)、零狀態(tài)響應(yīng)和零輸入響應(yīng)，程序如下：symsnzf(z)=2.^n;F(z)=ztrans(f(z))Ys=((5-6*z^-1)-6)/(1-5*z^-1+6*z^-2);ys=iztrans(Ys)Yf=z.^-1*F(z)/(1-5*z.^-1+6*z.^-2);yf=iztrans(Yf)yn=simplify(ys+yf)運(yùn)行結(jié)果為：F(z)=z/(z-2)ys=8*2^n-9*3^nyf=3*3^n-4*2^n-2^n*(n-1)yn=5*2^n-2^n*n-6*3^n即得系統(tǒng)的零輸入響應(yīng)為：零輸入響應(yīng)為：系統(tǒng)的全響應(yīng)為：例8.20用Matlab判斷例6.28中離散系統(tǒng)的穩(wěn)定性，系統(tǒng)的差分方程為程序如下：a=[1];b=[1/3,1/3,1/3];zplane(b,a);運(yùn)行結(jié)果如下：圖8.3.9例8.19離散系統(tǒng)的零極點(diǎn)圖從圖中可以看出，該系統(tǒng)有兩個零點(diǎn)，兩個極點(diǎn)，極點(diǎn)為0，在單位圓內(nèi)，故系統(tǒng)穩(wěn)定，與理論計算結(jié)果一致。例8.21用Matlab求例6.28中離散系統(tǒng)的系統(tǒng)函數(shù)，并判斷穩(wěn)定性，系統(tǒng)的差分方程為程序如下：a=[1,0.1,-0.2];b=[1,1];[r,p,k]=tf2zp(b,a)zplane(b,a);運(yùn)行結(jié)果為：r=-1p=-0.50000.4000k=1圖8.3.10例8.20離散系統(tǒng)的零極點(diǎn)圖由運(yùn)行結(jié)果知：該系統(tǒng)的零點(diǎn)為-1和0，極點(diǎn)為-0.5和0.4，系數(shù)為1，故得該系統(tǒng)的系統(tǒng)函數(shù)為：，由于極點(diǎn)都在單位圓內(nèi)，故系統(tǒng)穩(wěn)定，與例6.29理論計算結(jié)果一致。例8.22分別設(shè)計20階低通模擬chebyshevI型濾波器，通帶內(nèi)的最大衰減為0.3dB；20階低通模擬chebyshevII型濾波器，阻帶內(nèi)的最小衰減為45dB，并給出其頻率特性圖。程序如下：[z1,p1,k1]=cheb1ap(20,0.3);[num1,den1]=zp2tf(z1,p1,k1);%將零極點(diǎn)形式轉(zhuǎn)換為系統(tǒng)函數(shù)形式[z2,p2,k2]=cheb2ap(20,45);[num2,den2]=zp2tf(z2,p2,k2);figure(1)freqs(num1,den1)figure(2)freqs(num2,den2)運(yùn)行結(jié)果如下圖：圖8.4.3低通模擬chebyshevI型濾波器圖8.4.4低通模擬chebyshevII型濾波器例8.23設(shè)計數(shù)字信號處理系統(tǒng)，采樣頻率，希望設(shè)計成butterworth型高通數(shù)字濾波器，通帶內(nèi)允許的最大衰減為0.4dB，阻帶內(nèi)的最小衰減為40dB，通帶上線臨界頻率為400Hz，阻帶下限頻率為300Hz。解：設(shè)計數(shù)字高通濾波器的流程如下圖8.4.5：圖8.4.5數(shù)字高通濾波器設(shè)計流程Matlab程序如下：wp=400*2*pi;ws=300*2*pi;rp=0.4;rs=40;fs=2000;%轉(zhuǎn)換成模擬濾波器[N,Wc]=buttord(wp,ws,rp,rs,'s');%計算階數(shù)和截止頻率[Z,P,K]=buttap(N);%產(chǎn)生模擬低通濾波器[A,B,C,D]=zp2ss(Z,P,K);%零極點(diǎn)形式轉(zhuǎn)換為狀態(tài)方程形式[AT,BT,CT,DT]=lp2hp(A,B,C,D,Wc);%低通向高通的轉(zhuǎn)換[num1,den1]=ss2tf(AT,BT,CT,DT);%狀態(tài)方程形式向傳遞函數(shù)形式轉(zhuǎn)換[num2,den2]=bilinear(num1,den1,fs);%雙線性法得到的濾波器[H,W]=freqz(num2,den2);%離散系統(tǒng)頻率響應(yīng)F=W*fs/2/pi;semilogy(F,abs(H));grid%表示y坐標(biāo)軸是對數(shù)坐標(biāo)系xlabel('頻率（Hz）');ylabel('幅值');運(yùn)行結(jié)果如下：圖8.4.6例8.23數(shù)字高通濾波器仿真結(jié)果例8.24以中國移動的上行頻率為例（下行頻段設(shè)計類似），設(shè)計數(shù)字帶通濾波器，設(shè)抽樣頻率為，使用橢圓型濾波器。要求：①通帶范圍為885-909Hz，在通帶邊緣頻率出的衰減不大于0.5dB；②在780Hz以下和1000Hz以上衰減不小于15dB。Matlab程序如下：fs=2500;ws=[885,909]/(fs/2);wp=[780,1000]/(fs/2);rp=0.5;rs=15;[N,wc]=ellipord(wp,ws,rp,rs,'s');%計算階數(shù)和截止頻率[num,den]=ellip(N,rp,rs,wc);%設(shè)計橢圓濾波器[H,W]=freqz(num,den);F=W*fs/2/pi;plot(F,20*log10(abs(H)));gridxlabel('頻率（Hz）');ylabel('幅值');axis([02000-500])運(yùn)行結(jié)果如下：圖8.4.7例8.24帶通濾波器仿真結(jié)果例8.25用矩形窗設(shè)計線性相位FIR低通濾波器，通帶的截止頻率為0.2，濾波器單位脈沖響應(yīng)的長度為21。Matlab程序如下：M=21;wc=0.2*pi;r=(M-1)/2nr=-r:r;hdn=sin(wc*nr)/pi./nr;ifrem(M,2)~=0,hdn(r+1)=wc/pi;endw=boxcar(M);hn=hdn.*w';%相乘的時域沖激響應(yīng)hw=fft(hn,512);w1=2*[0:511]/512;%頻譜特性subplot(211),stem(1:M,hn,'o');xlabel('n'),ylabel('h(n)');subplot(212),plot(w1,20*log10(abs(hw)));xlabel('w/pi'),ylabel('幅值（dB）');axis([02-9010]);運(yùn)行結(jié)果如下圖：圖8.4.8例8.25FIR低通濾波器仿真結(jié)果例8.26設(shè)計一個多通帶濾波器，歸一化的通帶是：[0,0.2]、[0.4,0.6]、[0.8,1]。注意高頻端為通帶，故濾波器的階數(shù)應(yīng)該為偶數(shù)，取N=30。Matlab程序如下：wn=[0.2,0.4,0.6,0.8];N=30;b=fir1(N,wn,'dc-1');freqz(b)figure(2)stem(b,'.');xlabel('n'),ylabel('h(n)');運(yùn)行結(jié)果如下：圖8.4.9例8.26FIR多通帶濾波器頻譜特性圖8.4.10例8.26單位沖激響應(yīng)例8.27用fir2設(shè)計【例8.26】中的濾波器。Matlab程序如下：f=0:0.01:1m(1:21)=1;m(22:41)=0;m(42:61)=1;m(62:81)=0;m(82:101)=1;plot(f,m,'k');holdonN=30;b=fir2(N,f,m);[h,f1]=freqz(b);f1=f1./pi;plot(f1,abs(h));legend('理想濾波器','設(shè)計濾波器');xlabel('歸一化頻率'),ylabel('幅值');圖8.4.11例8.27濾波器結(jié)果例8.28已知三角脈沖信號為，試用Matlab實(shí)現(xiàn)該信號經(jīng)沖激脈沖抽樣后得到的抽樣信號及其頻譜。解：設(shè)采用的抽樣間隔為時，程序如下：Ts=1;dt=0.01;t1=-4:dt:4;xt=(4-abs(t1));subplot(221)plot(t1,xt),gridonaxis([-4.24.2-0.14.1]);xlabel('t'),title('三角脈沖信號')N=500;k=-N:N;W=pi*k/(N*dt);Xw=dt*xt*exp(-j*t1'*W);subplot(222)plot(W,abs(Xw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('三角脈沖信號頻譜')t2=-4:Ts:4;xst=(4-abs(t2));subplot(223)plot(t1,xt),holdonstem(t2,xst),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('經(jīng)過沖激抽樣后的信號'),holdoffXsw=Ts*xst*exp(-j*t2'*W);subplot(224)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('抽樣后信號的頻譜')運(yùn)行結(jié)果如下：圖8.4.13例8.28三角脈沖信號的沖激抽樣例8.29試用三角脈沖信號來驗(yàn)證抽樣定理。解：例8.27中三角脈沖信號的截止頻率為，故奈奎斯特間隔，在例8.27中可以通過修改的值來得到不同的結(jié)果。取時和取時結(jié)果如下：Ts=1;dt=0.01;t1=-4:dt:4;xt=(4-abs(t1));subplot(421)plot(t1,xt),gridonaxis([-4.24.2-0.14.1]);xlabel('t'),title('三角脈沖信號')N=500;k=-N:N;W=pi*k/(N*dt);Xw=dt*xt*exp(-j*t1'*W);subplot(422)plot(W,abs(Xw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('三角脈沖信號頻譜')t2=-4:Ts:4;xst=(4-abs(t2));subplot(423)plot(t1,xt),holdonstem(t2,xst),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('經(jīng)過沖激抽樣后的信號'),holdoffXsw=Ts*(4-abs(t2))*exp(-j*t2'*W);subplot(424)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('抽樣后信號的頻譜')t3=-4:2:4;subplot(425)plot(t1,xt),holdonstem(t3,(4-abs(t3))),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('t=2經(jīng)過沖激抽樣后的信號'),holdoffXsw=2*(4-abs(t3))*exp(-j*t3'*W);subplot(426)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('t=2抽樣后信號的頻譜')t4=-4:3:4;subplot(427)plot(t1,xt),holdonstem(t4,(4-abs(t4))),gridonaxis([-4.24.2-0.14.1]),xlabel('t'),title('t=3經(jīng)過沖激抽樣后的信號'),holdoffXsw=3*(4-abs(t4))*exp(-j*t4'*W);subplot(428)plot(W,abs(Xsw)),gridonaxis([-1010-0.118]),xlabel('\omega'),title('t=3抽樣后信號的頻譜')運(yùn)行結(jié)果如下：圖8.4.14例8.29三角脈沖信號的沖激抽樣（Ts取不同的值）觀察圖8.4.14可以發(fā)現(xiàn)，抽樣時間間隔大于奈奎斯特間隔時，信號的頻譜會產(chǎn)生十分嚴(yán)重的混疊現(xiàn)象。例8.30對于三角脈沖信號，假設(shè)截止頻率為，抽樣間隔為，試用Matlab恢復(fù)抽樣信號，并計算恢復(fù)后的信號與原信號的絕對誤差。解：去低通濾波器的截止頻率為，則Matlab程序如下：wm=pi/2;wc=1.1*wm;Ts=1;n=-50:50;nTs=n*Ts;xs=(4-abs(nTs));t=-4:0.01:4;xt=xs*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));t1=-4:0.01:4;x1=4-abs(t1);subplot(311)stem(nTs,xs),holdonplot(t1,x1),gridonaxis([-4.14.1-0.14.1])title('Ts=1時的抽樣信號')holdoffsubplot(312)plot(t,xt),gridonaxis([-4.14.1-0.14.1])title('恢復(fù)后所得的三角脈沖信號')error=abs(xt-x1);subplot(313)plot(t,error),gridontitle('恢復(fù)信號與原信號的絕對誤差')運(yùn)行結(jié)果如下：圖8.4.15例8.30抽樣信號的恢復(fù)例8.31對信號、進(jìn)行幅度調(diào)制。程序如下：fs=100;fc=15;t=0:0.01:2;x=sin([pi*t',2*pi*t']);y=ammod(x,fc,fs);z=amdemod(y,fc,fs);subplot(211)plot(t,x(:,1),'-',t,z(:,1),'o',t,y(:,1),'r')title('sin(pi*t)')subplot(212)plot(t,x(:,2),'--',t,z(:,2),'*',t,y(:,2),'r')title('sin(2*pi*t)')運(yùn)行結(jié)果如下：圖8.4.16例8.31幅度調(diào)制信號例8.32對信號、進(jìn)行頻率調(diào)制和解調(diào)。程序如下：fs=2000;fc=750;t=[0:fs]'/fs;s1=sin(2*pi*300*t)+2*sin(2*pi*600*t);s2=sin(2*pi*150*t)+2*sin(2*pi*900*t);x=[s1,s2];dev=50;y=fmmod(x,fc,fs,dev);z=fmdemod(y,fc,fs,dev);subplot(211)plot(t,x(:,1),'-',t,z(:,1),'o',t,y(:,1),'r')title('s1')subplot(212)plot(t,x(:,2),'--',t,z(:,2),'*',t,y(:,2),'r')title('s2')運(yùn)行結(jié)果如下：圖8.4.17例8.32頻率調(diào)制信號例8.33對模擬信號進(jìn)行相位調(diào)制，經(jīng)信道疊加白噪聲，解調(diào)并繪制原始信號與解調(diào)后信號波形。程序如下：fs=100;fc=10;t=[0:2*fs+1]'/fs;x=sin(2*pi*t)+sin(4*pi*t);phasedev=pi/2;y=pmmod(x,fc,fs,phasedev);y=awgn(y,10,'measured',103)%疊加噪聲z=pmdemod(y,fc,fs,phasedev);figure;plot(t,x,'k-',t,z,'r-');legend('originalsignal','recoveredsignal')運(yùn)行結(jié)果如下：圖8.4.18例8.33相位調(diào)制解調(diào)信號例8.34給一鋸齒波信號添加白噪聲，并繪制出原信號和含噪聲信號的波形。程序如下：t=0:0.1:10;x=sawtooth(t);y=awgn(x,10,'measured');plot(t,x,'*',t,y,'r')legend('originalsignail','signalwithawgn')運(yùn)行結(jié)果如下：圖8.4.18例8.34添加高斯白噪聲的信號與原信號波形例8.35利用Simulink仿真實(shí)現(xiàn)對隨機(jī)二進(jìn)制信號（即隨機(jī)出現(xiàn)的矩形脈沖函數(shù)）x(t)用cos200t進(jìn)行調(diào)制，疊加隨機(jī)噪聲，再進(jìn)行相干解調(diào)。觀察各部分信號的情況及噪聲對解調(diào)后的信號產(chǎn)生的影響。系統(tǒng)仿真模型圖如下圖8.4.19所示。圖8.4.19例8.35隨機(jī)矩形脈沖信號調(diào)制-解調(diào)Simulink仿真模型例8.36讀取電腦中一段已經(jīng)存在一段wav格式的語音信號語音一。程序如下：clc%清空當(dāng)前工作區(qū)域[data,fs]=audioread('C:\Users\G\Documents\錄音\ly1.wav');%讀取數(shù)據(jù)left=data(:,1);right=data(:,2);%分離出左/右聲道數(shù)據(jù)time=(1/fs)*length(left);%劃分時間間隔t=linspace(0,time,length(left));%給定時間范圍及步進(jìn)plot(t,data)%做出頻譜圖像;soundsc(data,fs)%音頻播放xlabel('time(sec)');ylabel('relativesignalstrength');%定義坐標(biāo)軸表示運(yùn)行結(jié)果可播放音頻，頻譜如下圖8.4.22：圖8.4.22音頻信號頻譜圖例8.37對語音信號一和語音信號二進(jìn)行預(yù)處理。程序如下：[y1,fs]=audioread('C:\Users\G\Documents\錄音\ly1.wav');%讀取語音一信號[y2,fs]=audioread('C:\Users\G\Documents\錄音\ly2.wav');%讀取語音二信號L1=length(y1);%測定語音一信號長度L2=length(y2);%測定語音二信號長度a1=y1(:,1).*hamming(L1);%加窗預(yù)處理a2=y2(:,1).*hamming(L2);%加窗預(yù)處理L1=length(a1);%測定語音一信號長度L2=length(a2);%測定語音二信號長度%采樣信號的時域顯示figure(1);subplot(211);plot(a1);title('語音一載波信號時域波形');subplot(212);plot(a2);title('語音二調(diào)幅信號時域波形');%傅里葉頻譜繪制F1=fft(a1,L1);F2=fft(a2,L2);AF1=abs(F1);AF2=abs(F2);figure(2);subplot(211);plot(AF1);title('語音一載波信號幅頻特性顯示');subplot(212);plot(AF2);title('語音二調(diào)幅信號幅頻特性顯示');figure(3);freqz(F1);title('語音一載波信號FFT頻譜顯示');figure(4);freqz(F2);title('語音二載波信號FFT頻譜顯示');圖8.4.23信號預(yù)處理之后時域波形圖8.4.24信號預(yù)處理之后頻域波形圖8.4.25載波信號FFT頻譜顯示例8.38對例8.37中完成預(yù)處理的信號獲取初始位置并進(jìn)行信號合成。程序如下：%獲取語音一信號的開始位置fori=1:L1-4g(i)=a1(i).*a1(i+1).*a1(i+2).*a1(i+3).*a1(i+4);%認(rèn)為連續(xù)4個幅值不為0的信號即為開始ifg(i)~=0break;elsei=i+1;endendI=i;%獲取語音二信號開始位置forj=1:L2-4m(j)=a2(j).*a2(j+1).*a2(j+2).*a2(j+3).*a2(j+4);ifm(j)~=0break;elsej=j+1;endendJ=j;%語音二信號hilbert變換H=hilbert(a2);figure(5);plot(abs(H));title('語音二信號包絡(luò)顯示');%信號對齊,語音二包絡(luò)調(diào)制語音一振幅max1=max(I,J);fork=1:5*max1N(k)=a1(i).*H(j);i=i+1;j=j+1;end%N=N';N=N/(max(abs(N))*1.05);audiowrite('HC.wav',N,fs);%存儲合成語音soundsc(imag(N),fs)%播放合成信號figure(6);plot(imag(N));title('合成信號時域顯示');pause(1);FN=fft(N);figure(7);freqz(FN);title('合成聲音信號FFT顯示');figure(8);plot(abs(FN));title('合成聲音信號的幅頻特性');圖8.4.26語音二信號包絡(luò)圖合成信號的時域顯示結(jié)果、幅頻特性、快速傅里葉變換結(jié)果分別如圖8.4.27-8.4.29所示，該合成信號是以語音一信號的特性和語音二信號的幅度變化的，由其快速傅里葉變換的結(jié)果更證實(shí)了這一點(diǎn)，其幅頻特性與語音一信號的幅頻特性更接近。圖8.4.27合成語音信號的時域波形圖8.4.28合成語音信號的幅頻特性圖8.4.29合成語音信號快速傅里葉變換結(jié)果4、重建的趣味語音functionct2%定義常數(shù)FL=80;%幀長WL=240;%窗長P=10;%預(yù)測系數(shù)個數(shù)hw=hamming(WL);%漢明窗[s,Fs]=audioread('C:\Users\G\Documents\錄音\ly1.wav');%載入語音sload('mtlb.mat');s=mtlb;s=buffer(s,FL);%分割信號x成不重疊的長度為FL的數(shù)據(jù)幀[~,FN]=size(s);%計算幀數(shù)%預(yù)測和重建濾波器exc=zeros(FL,FN);%激勵信號（預(yù)測誤差）s_rec=exc;%重建語音zi_pre=zeros(P,1);%預(yù)測濾波器的狀態(tài)zi_rec=zeros(P,1);%重建濾波器的狀態(tài)%合成濾波器exc_syn=exc;%合成的激勵信號（脈沖串）s_syn=exc;%合成語音last_syn=0;%存儲上一個幀的最后一個脈沖在下一幀結(jié)束的位置zi_syn=zeros(P,1);