例2.4參考代碼：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%a.利用conv函數(shù)T=0.1;%%%時間步長t1=0:T:10;%%%時間序列f1=cos(t1);%%%信號f1t2=t1;f2=exp(t2)+exp(-2*t2);%%%信號f2f=T*conv(f1,f2);%%%計(jì)算卷積k0=t1(1)+t2(1);%%%卷積輸出序列的起始k3=length(f1)+length(f2)-2;t=k0:T:(k0+T*k3);%%%卷積結(jié)果對應(yīng)的時間向量%subplot(3,1,1);%%%繪制信號f1plot(t1,f1,'linewidth',2);title('f1(t)');subplot(3,1,2);%%%繪制信號f2plot(t2,f2,'linewidth',2);title('f2(t)');subplot(3,1,3);%%%繪制卷積結(jié)果plot(t,f,'linewidth',2);title('convolutionoff1(t)andf2(t)');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%b.編程實(shí)現(xiàn)T=0.1;%%%時間步長t1=0:T:10;%%%時間序列f1=cos(t1);%%%信號f1t2=t1;f2=exp(t2)+exp(-2*t2);%%%信號f2lf1=length(f1);%%%信號f1長度lf2=length(f2);%%%信號f2長度fork=1:lf1+lf2-1y(k)=0;%%%y賦初始值forii=max(1,k-(lf2-1)):min(k,lf1)y(k)=y(k)+f1(ii)*f2(k-ii+1);%%%信號相乘和求和endyzsappr(k)=T*y(k);%%%用乘和加運(yùn)算來近似積分運(yùn)算endt0=t1(1)+t2(1);%%%卷積輸出序列起點(diǎn)t3=lf1+lf2-2;t=t0:T:(t0+t3*T);%%%卷積輸出對應(yīng)的時間序列subplot(3,1,1);%%%繪制信號f1的波形plot(t1,f1,'linewidth',2);title('f1(t)');subplot(3,1,2);%%%繪制信號f2的波形plot(t2,f2,'linewidth',2);title('f2(t)');subplot(3,1,3);%%%繪制卷積結(jié)果的波形plot(t,f,'linewidth',2);title('convolutionoff1(t)andf2(t)');例3.1%例3.1求解周期信號傅里葉級數(shù)的MATLAB程序%先求出信號傅里葉級數(shù)的系數(shù)symsTEtnx;pi=sym('pi');%創(chuàng)建符號對象a0=2/T*int(-E,t,-T/2,0)+2/T*int(E,t,0,T/2),A0=a0/2;%計(jì)算系數(shù)an=2/T*int(-E*cos(2*pi*n*t/T),t,-T/2,0)+2/T*int(E*cos(2*pi*n*t/T),t,0,T/2)bn=2/T*int(-E*sin(2*pi*n*t/T),t,-T/2,0)+2/T*int(E*sin(2*pi*n*t/T),t,0,T/2)%求出E=1時各諧波分量的幅度symsEAncbn=2*E*c;E=1;c=-(cos(n*pi)-1)/(n*pi);bn=subs(bn,{sym('E'),sym('c')},{E,c})%符號變量置換%畫出頻譜圖A=[00000000000];%存儲11個諧波分量幅度的數(shù)組n=1:1:11;forn=1:11bn=2*(-cos(n*pi)+1)/(n*pi);A(n)=double(vpa(bn));%對任意精度的符號類數(shù)據(jù)進(jìn)行規(guī)范An=A(n)endx1=1;x2=11;x=x1:x2;stem(x,A,'r','filled')%畫出11個諧波分量的幅度譜例3.3%產(chǎn)生周期為2π的方波信號t=0:1/10:10;y=square(2*pi*(0.05*pi)*t);subplot(2,2,1);plot(t,y);xlabel('(a)方波波形')axis([010-1.51.5])%基波信號t=0:1/10:10;y=sin(t);subplot(2,2,2);plot(t,y);xlabel('(b)基波波形')%基波+3次諧波合成的波形t=0:1/10:10;y=sin(t)+sin(3*t)/3;subplot(2,2,3);plot(t,y);xlabel('(c)基波+3次諧波')%基波+3次+5次+7次+9次諧波合成的波形t=0:1/10:10;y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+sin(9*t)/9;subplot(2,2,4);plot(t,y);xlabel('(d)基波+3次+5次+7次+9次諧波') %最高諧波階次為19時的吉伯斯現(xiàn)象t=0:31/1000:5;y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:19,x=x+sin(k*t)/k;y((k+1)/2,:)=x;end;pause,plot(y(1:9,:)),pause,mesh(y),pauseclc%最高諧波階次為45時的吉伯斯現(xiàn)象y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:45,x=x+sin(k*t)/k;y((k+1)/2,:)=x;end;pause,plot(y(1:9,:)),pause,mesh(y),pauseclc例3.10%例3.10中傅里葉變換的MATLAB程序symsEtaw;f=E*Heaviside(t+a)-E*Heaviside(t-a)F=fourier(f);F1=simplify(F)%得到F的簡化形式E=1;a=1;F2=subs(F1,{sym('E'),sym('a')},{Ea})%將F1中的變量替換為新變量ezplot(abs(F2),[-3*pi,3*pi])title('F(Ω)')xlabel('Ω')f2=ifourier(F2)運(yùn)行結(jié)果為F=(E*(sin(a*w)+cos(a*w)*1i))/w-(E*(cos(a*w)*1i-sin(a*w)))/w例3.14%例3.14高斯信號的頻譜symstwsigma=sym('sigma','positive');f=exp(-t^2/(2*sigma^2))/(sqrt(2*pi)*sigma);F=fourier(f,t,w);F1=simplify(F);digits(4);F2=vpa(F1)sigma=1;F3=subs(F2,'sigma',sigma)ezplot(F3,[-2*pi,2*pi]),gridontitle('F(Ω)');xlabel('Ω')例3.15%信號最高頻率wm＝2ws（ws:采樣頻率）臨界采樣的情況wm=1;%信號最高頻率（帶寬）wc=wm;%預(yù)處理用低通濾波器截止頻率Ts=pi/wm;%采樣間隔ws=2*pi/Ts;%變?yōu)椴蓸咏穷l率n=-100:100;%采樣點(diǎn)數(shù)nTs=n*Ts;f=sinc(nTs/pi);%對信號Sa(t)的抽樣信號Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));%信號的重構(gòu)t1=-10:0.5:10;f1=sinc(t1/pi);%Sa(t)的連續(xù)波形表示%畫出抽樣信號與抽樣信號的重構(gòu)波形subplot(2,1,1);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的臨界抽樣');subplot(2,1,2);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由抽樣信號重構(gòu)sa(t)');grid%欠抽樣情況，抽樣頻率ws與信號最高頻率wm之間，wm>2ws.wm=1;wc=1.1*wm;Ts=2*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Ts;f=sinc(nTs/pi);Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));error=abs(fa-sinc(t/pi));%求解重構(gòu)信號與原信號的誤差t1=-10:0.5:10;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的欠抽樣信號');subplot(312);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由欠抽樣信號重構(gòu)的sa(t)');grid;subplot(313);plot(t,error);xlabel('t');ylabel('error(t)');title('欠抽樣信號與原信號的誤差error(t)')%過抽樣情況，wm<2ws.wm=1;wc=1.1*wm;Ts=0.5*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Ts;f=sinc(nTs/pi);Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));error=abs(fa-sinc(t/pi));t1=-10:0.5:10;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的采樣信號');subplot(312);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由過抽樣信號重構(gòu)的sa(t)');grid;subplot(313);plot(t,error);xlabel('t');ylabel('error(t)');title('過抽樣信號與原信號的誤差error(t)')例4.17%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%symst;%%%定義符號變量f=cos(2*t);%%%信號表達(dá)式F=laplace(f);%%%調(diào)用MATLAB函數(shù)對x做拉普拉斯變換fl=cos(2*(t-1));%%%信號xl表達(dá)式Fl=laplace(fl);%%%對xl做拉普拉斯變換例4.18%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%symss;%%%定義符號變量F=5*(s+2)*(s+5)/(s*(s+1)*(s+3));%%%信號表達(dá)式f=ilaplace(F);%%%逆拉普拉斯變換例4.19%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%A=[83-21];%%%分子系數(shù)B=[10-7-6];%%%分母系數(shù)z=roots(A);%%%解得零點(diǎn)p=roots(B);%%%解得極點(diǎn)x=max(abs([z;p]));%%%取零極點(diǎn)的最大值x=x+0.1;y=x;figure;%%%繪制零極點(diǎn)圖holdon;%%%使多次繪圖同時保留plot([-xx],[00],'--');%%%繪制實(shí)軸（或x軸）plot([00],[-yy],'--');%%%繪制虛軸（或y軸）plot(real(z),imag(z),'bo',real(p),imag(p),'kx');xlabel('RealPart');ylabel('ImaginaryPart');axis([-xx-yy]);例4.20%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%symst;%%%定義符號變量f=exp(-2*t)*heaviside(t);%%%信號表達(dá)式F=laplace(f);%%%對信號做拉普拉斯變換%繪制零極點(diǎn)圖%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%A=[1];%%%多項(xiàng)式分子B=[12];%%%多項(xiàng)式分母z=roots(A);%%%解得零點(diǎn)p=roots(B);%%%解得極點(diǎn)x=max(abs([z;p]));%%%取零極點(diǎn)的最大值x=x+0.1;y=x;figure;%%%繪制零極點(diǎn)圖holdon;plot([-xx],[00],'--');%%%繪制實(shí)軸plot([00],[-yy],'--');%%%繪制虛軸plot(real(z),imag(z),'bo',real(p),imag(p),'kx');xlabel('RealPart');ylabel('ImaginaryPart');axis([-xx-yy]);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%length=100;%%%定義仿真的長度xs=zeros(length,length);%%%仿真區(qū)域初始化delta_t=0.1;%%%時間步長formm=1:length%%%對仿真區(qū)域每一個點(diǎn)fornn=1:lengthdelta_xs=0;%%%保存積分運(yùn)算的中間結(jié)果s=(mm/10-5)+1j*(nn/10-5)%%%s的取值范圍為[-5-5i,5+5i]xs(mm.nn)=1./(s+2);%%%信號的拉普拉斯變換ifmm<30%%%相當(dāng)于sigma<-2fortt=0:delta_t:50delta_xs=delta_xs+exp(-(2+s)*tt);%%%對delta_xs累加endxs(mm,nn)=delta_xs*delta_t;%%%對xt進(jìn)行積分endendendxs=log10(abs(xs));figure;mesh(1:length,1:length,xs);set(gca,'xtick',[0102030405060708090100]);set(gca,'xticklabel',{'-5','-4','-3','-2','-1','0','1','2','3','4','5'});set(gca,'ytick',[0102030405060708090100]);set(gca,'yticklabel',{'-5','-4','-3','-2','-1','0','1','2','3','4','5'});xlabel('j\omega');ylabel('sigma');zlabel('log10(X(s))');例5.3%MATLABPROGRAMExample5.3b=input('差分方程左端系數(shù)b=[b(1),…,b(nb)]=');a=input('差分方程右端系數(shù)a=[a(1),...,a(na)]=');x=input('輸入信號序列x=');nb=length(b);na=length(a);nx=length(x);s=['現(xiàn)時刻前',int2str(nb-1),'點(diǎn)ym的值=[ym(1),...,ym(na-1)='];ym=zeros(1,nb+nx-1);%先預(yù)設(shè)ym序列長度為nb+nx,并初始化為0ym(1:nb-1)=input(s);%對ym序列的前nb個點(diǎn)賦初值xm=[zeros(1,nb-1),x];%給序列xm賦初值forn=nb:nb+nx-1%n以ym的起點(diǎn)為初值ys=ym(n-1:-1:n-nb+1);%得出ysxs=xm(n-1:-1:n-na);%得出xsym(n)=(a*xs'-b(2:nb)*ys')/a(1);%遞推求ymend%把ym左移nb-1位,求得yy=ym(nb:nb+nx-1);stem(y),gridline([0,60],[0,0])%給出起點(diǎn)、終點(diǎn)坐標(biāo)，畫橫坐標(biāo)例6.2%x(n)的圖形表示n=-5:5;x=(-0.95).^n;subplot(3,1,1);stem(n,x);title('序列x(n)的圖形');xlabel('n');ylabel('x(n)')%x(n)的傅里葉變換k=-200:200;w=(pi/100)*k;%橫坐標(biāo)(頻率軸)的示值縮小比例系數(shù)X=x*(exp(-j*pi/100)).^(n'*k);magX=abs(X);angX=angle(X);subplot(3,1,2);plot(w/pi,magX);gridon;title('x(n)的幅度譜');axis([-2,2,0,15]);xlabel('頻率ω(單位:π)');ylabel('|X|')subplot(3,1,3);plot(w/pi,angX/pi);gridon;title('x(n)的相位譜');axis([-2,2,-1,1]);xlabel('頻率ω(單位:π)');ylabel('φ/π')例6.4%構(gòu)成并畫出受噪聲污染的信號clear,randn('state',0);%正態(tài)隨機(jī)數(shù)發(fā)生器置0t=linspace(0,10,512);%生成均勻采樣數(shù)組y=6*cos(6*t)-4*sin(12*t)+6*randn(size(t));%生成波形和%標(biāo)準(zhǔn)差為6，標(biāo)準(zhǔn)偏差為0的正態(tài)隨機(jī)噪聲subplot(311);plot(t,y)xlabel('t');%標(biāo)注坐標(biāo)軸ylabel('信號值f(t)');%計(jì)算并畫出受噪聲污染信號的幅度譜Y=fft(y);%計(jì)算信號的DFTTs=t(2)-t(1)%計(jì)算信號的采樣周期Ws=2*pi/Ts;%計(jì)算信號的采樣頻率Wn=Ws/2w=linspace(0,Wn,length(t)/2);%頻率軸的頻率間隔Ya=abs(Y(1:length(t)/2));%頻譜幅度subplot(312);plot(w,Ya)xlabel('\omega')%標(biāo)注坐標(biāo)軸ylabel('頻譜幅值F(\omega)')ii=find(w<=20);%在有效信號頻率段尋找數(shù)組中非零元素序號subplot(313);plot(w(ii),Ya(ii))xlabel('\omega')%標(biāo)注坐標(biāo)軸ylabel('局部放大譜幅F(\omega)')grid例6.5clear%序列的表示a=ones(1,15);a([1,2,3,4,5])=0;%序列a(n)b=ones(1,12);b([1,2,3])=0;%序列b(n)%序列直接卷積c=conv(a,b);%直接卷積運(yùn)算%應(yīng)用FFT及DFT性質(zhì)，進(jìn)行序列的快速卷積M=32;%加長序列符合快速卷積要求AF=fft(a,M);%序列a(n)FFT得AFBF=fft(b,M);%序列b(n)FFT得BFCF=AF.*BF;%AF和BF相乘得CFcc=real(ifft(CF));%CF進(jìn)行IFFT%求直接卷積和快速卷積兩種計(jì)算結(jié)果之間的誤差nn=0:(M-1);c(M)=0;error=c-cc;%圖形表示出序列、直接卷積和快速卷積以及兩種卷積結(jié)果的誤差subplot(411);stem([1:15],a),ylabel('a(n)');subplot(412);stem([1:12],b),ylabel('b(n)');subplot(413),stem(nn,c,'fill'),grid,axis([0,31,0,9])ylabel('a(n)*b(n)')subplot(414),stem(nn,error,'fill'),axis([0,31,-1,1])xlabel('n'),ylabel('誤差')例6.6N=64;T=5/N;t=0:T:5;x=exp(-1.5*t);%連續(xù)時間信號x(t)Xx=T*fft(x,N);%x(t)的近似頻譜n=0:1:63;x0=1/T.*ifft(Xx,N);%重構(gòu)x(n)%分別畫出x(t)波形、幅譜X(k)和x(n)subplot(311);plot(t,x);xlabel('t');ylabel('x(t)')subplot(312);stem(n,Xx);xlabel('k'),ylabel('X(k)')subplot(313)stem(n,x0);xlabel('n'),ylabel('x(n)')例6.7%example6.7MATLABPROGRAMN=256;%采樣點(diǎn)數(shù)dt=0.0004%采樣時間間隔wn=500;%二階系統(tǒng)的固有頻率seta=0.47;%系統(tǒng)阻尼比dw=2*pi/(N*dt);%頻率間隔a=wn^2;b=[1,2*seta*wn,a];t=[0:dt:(N-1)*dt];c=step(a,b,t);%求二階系統(tǒng)的階躍響應(yīng)w=[0:dw:(N-1)*dw];[mag,phase]=bode(a,b,w);%計(jì)算二階系統(tǒng)理論頻率特性ycw=fft(c);%求系統(tǒng)階躍響應(yīng)瞬態(tài)分量的傅里葉變換Re=real(ycw);%ycw的實(shí)部Im=imag(ycw);%ycw的虛部fori=1:NRw(i,1)=1-Im(i,1)*(i-1)*dw*dt;Iw(i,1)=Re(i,1)*(i-1)*dw*dt;end%計(jì)算頻率特性的實(shí)部和虛部分量ffw=Rw+Iw*sqrt(-1);Aw=abs(ffw)%系統(tǒng)幅頻特性semilogx(w,20*log10(mag),'r-')%理論幅頻特性曲線axis([100,10000,-30,10])text(600,12,'對數(shù)幅頻特性')ylabel('A(ω)(dB)');holdonsemilogx(w,20*log10(Aw))%由階躍響應(yīng)求得的幅頻特性曲線axis([100,10000,-30,10])xlabel('ω(1/s)');gridon例7.11%MATLABPROGRAMofExample7.11b=[0,1];a=[3,-4,1];[R,p,k]=residuez(b,a)例7.12%例7.12中求解z變換的MATLAB程序symsaznTf=a^(n*T)F=factor(ztrans(f))%做因式分解處理%若T=1,a=1求其z變換symsknwzF=ztrans(1^n)例7.13symsznabF=k*z^2/((z-a)*(z-b))f=iztrans(F)%若k=1,a=1,b=0.5時,求其z反變換symsnzF=z^2/((z-1)*(z-0.5))f=iztrans(F)例7.18%Exampe7.18MATLABProgram%畫出單位響應(yīng)h(n)clear;n=0:7;hn=(0.9).^n;subplot(3,1,1);stem(n,hn,'k')xlabel('n');ylabel('h(n)');%畫出頻率響應(yīng)（幅頻和相頻）k=0:200;w=(pi/100)*k;Hew=hn*(exp(-j*pi/100).^(n'*k));mag_H=abs(Hew);ang_H=angle(Hew);subplot(3,1,2);plot(w/pi,mag_H,'k');xlabel('w');ylabel('H(w)');subplot(3,1,3);plot(w/pi,ang_H/pi,'k');xlabel('w');ylabel('\Phi');例8.4%例8.4的MATLAB程序n=0:0.01:2;fori=1:4switchicase1N=2;case2N=5;case3N=10;case4N=20;end[z,p,k]=buttap(N);%巴特沃思濾波器原型設(shè)計(jì)函數(shù)，n:階數(shù)；%z,p,k：濾波器零點(diǎn)、極點(diǎn)和增益[b,a]=zp2tf(z,p,k);%零極點(diǎn)增益轉(zhuǎn)換為傳遞函數(shù)[H,w]=freqs(b,a,n);%模擬濾波器頻率響應(yīng)函數(shù)magH2=(abs(H)).^2;%幅度平方函數(shù)holdonplot(w,magH2);axis([0201]);endxlabel('w/wc');ylabel('|H(jw)|^2');grid例8.5%繪制切比雪夫低通濾波器幅度平方函數(shù)的NATLAB程序n=0:0.01:2;fori=1:4switchicase1N=2;case2N=4;case3N=6;case4N=8;endRp=1;%濾波器通帶波紋系數(shù)[z,p,k]=cheb1ap(N,Rp);%設(shè)計(jì)切比雪夫模擬低通濾波器原型函數(shù),%z,p,k分別為濾波器的零點(diǎn)、極點(diǎn)和增益[b,a]=zp2tf(z,p,k);%零點(diǎn)、極點(diǎn)和增益轉(zhuǎn)換為傳遞函數(shù)[H,w]=freqs(b,a,n);%模擬濾波器的頻率響應(yīng)magH2=(abs(H)).^2;%幅度平方函數(shù)%Outputposplot=['22',num2str(i)];%定義字符串變量subplot(posplot)plot(w,magH2,'k');ylabel('|H(jw)^2');title(['N=',num2str(N)]);gridend例8.7%設(shè)計(jì)切比雪夫模擬低通濾波器的MATLAB程序wp=200*pi;ws=300*pi;Rp=1;%通帶波紋Rs=16;%阻帶衰減%計(jì)算濾波器階數(shù)ebs=sqrt(10^(Rp/10)-1);%波紋系數(shù)A=10^(Rs/20);%A為參變量Wc=wp;%截止頻率Wr=ws/wp;g=sqrt(A*A-1)/ebs;%g為參變量N1=log10(g+sqrt(g*g-1))/log10(Wr+sqrt(Wr*Wr-1));%濾波器階數(shù)計(jì)算N=ceil(N1);%N應(yīng)取整例8.8%設(shè)計(jì)巴特沃思帶通濾波器MATLAB程序%濾波器設(shè)計(jì)指標(biāo)wp=[20003000]*2*pi;ws=[15003500]*2*pi;Rp=1;Rs=100;%計(jì)算階數(shù)和截止頻率[N,Wn]=buttord(wp,ws,Rp,Rs,'s');Fc=Wn/(2*pi);%計(jì)算濾波器傳遞函數(shù)多項(xiàng)式系數(shù)[b,a]=butter(N,Wn,'s');%畫出濾波器幅頻特性w=linspace(1,4000,1000)*2*pi;%生成線性等間隔的向量H=freqs(b,a,w);magH=abs(H);phaH=unwrap(angle(H));plot(w/(2*pi),20*log10(magH),'k');xlabel('頻率(Hz)');ylabel('幅度(dB)');title('巴特沃思模擬濾波器')gridon例8.9%設(shè)計(jì)切比雪夫帶阻濾波器MATLAB程序%濾波器設(shè)計(jì)指標(biāo)wp=[20003000]*2*pi;ws=[15003500]*2*pi;Rp=1;Rs=60;%計(jì)算階數(shù)和截止頻率[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');%計(jì)算濾波器傳遞函數(shù)多項(xiàng)式系數(shù)[b,a]=cheby1(N,Rp,Wn,'stop','s');%畫出濾波器幅頻特性w=linspace(1,4000,1000)*2*pi;%生成線性等間隔的向量H=freqs(b,a,w);%相位展開magH=abs(H);phaH=unwrap(angle(H));plot(w/(2*pi),20*log10(magH),'k');xlabel('頻率(Hz)');ylabel('幅度(dB)');title('切比雪夫模擬濾波器')gridon例8.15%例8.15中的MATLAB程序%定義濾波器的性能指標(biāo)wp=2000*2*pi;ws=4000*2*pi;Rp=3;Rs=15;Fs=20000;Nn=128;%給出濾波器的序列點(diǎn)數(shù)%計(jì)算濾波器階數(shù)和截止頻率[N,Wn]=buttord(wp,ws,Rp,Rs,'s');%設(shè)計(jì)模擬濾波器原型[z,p,k]=buttap(N);[Bap,Aap]=zp2tf(z,p,k);%將系統(tǒng)函數(shù)的零極點(diǎn)形式變?yōu)榉肿臃帜付囗?xiàng)式形式[b,a]=lp2lp(Bap,Aap,Wn);%低通至低通模擬濾波器的變換%沖激響應(yīng)不變法實(shí)現(xiàn)模擬數(shù)字濾波器的變換[bz,az]=impinvar(b,a,Fs);freqz(bz,az,Nn,Fs)例8.19%例8.19雙線性變換法設(shè)計(jì)數(shù)字低通濾波器MATLAB程序%給定濾波器指標(biāo)wp=0.2*pi;ws=0.3*pi;Rp=3;Rs=15;Ts=0.001;Nn=128;%數(shù)字頻率模擬頻率非線性變換Wp=(2/Ts)*tan(wp/2);Ws=(2/Ts)*tan(ws/2);%計(jì)算濾波器階次和截止頻率[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s');%設(shè)計(jì)模擬原型[z,p,k]=buttap(N);[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2lp(Bap,Aap,Wn);%雙線性變換法設(shè)計(jì)數(shù)字濾波器[bz,az]=bilinear(b,a,1/Ts);freqz(bz,az,Nn,1/Ts)例8.22%給出濾波器設(shè)計(jì)要求Fs=10000;wp=[20003000]*2/Fs;ws=[15003500]*2/Fs;Rp=1;Rs=20;Nn=128;%計(jì)算模擬原型的階次和截止頻率[N,Wn]=cheb1ord(wp,ws,Rp,Rs);[b,a]=cheby1(N,Rp,Wn);%得到頻響特性freqz(b,a,Nn,Fs)例8.23%例8.23中設(shè)計(jì)巴特沃思高通數(shù)字濾波器MATLAB程序%給定濾波器設(shè)計(jì)指標(biāo)Fs=5000;wp=2000*2/Fs;ws=1500*2/Fs;Rp=1;Rs=20;Nn=128;%計(jì)算濾波器階次和截止頻率[N,Wn]=buttord(wp,ws,Rp,Rs);%設(shè)計(jì)數(shù)字濾波器[b,a]=butter(N,Wn,'high');%數(shù)字濾波器頻響特性freqz(b,a,Nn,Fs)例8.25%濾波器設(shè)計(jì)要求wn=[0.280.58];N=50;%使用FIR1類函數(shù)計(jì)算并求出濾波器特性b=fir1(2*N,wn);freqz(b,1,512)例8.27%例8.27中設(shè)計(jì)FIR濾波器MATLAB程序N=20;alpha=(N-1)/2;l=0:N-1;wl=(2*pi/N)*l;Hrs=[1,1,1,zeros(1,15),1,1];Hdr=[1,1,0,0];wdl=[0,0.25,0.25,1];k1=0:floor((N-1)/2);k2=floor((N-1)/2)+1:N-1;angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];H=Hrs.*exp(j*angH);h=real(ifft(H,N));[db,mag,pha,grd,w]=freqz_m(h,1);[Hr,ww,a,L]=Hr_Type2(h);subplot(1,1,1)subplot(2,2,1);plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);axis([0,1,-0.1,1.1]);title('頻率樣本:N=20')xlabel('頻率(單位:pi)');ylabel('Hr(k)')set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])set(gca,'YTickMode','manual','YTick',[0,1]);gridsubplot(2,2,2);stem(l,h);axis([-1,N,-0.1,0.3])title('單位抽樣響應(yīng)');ylabel('h(n)');text(N+1,-0.1,'n')subplot(2,2,3);plot(ww/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');axis([0,1,-0.2,1.2]);title('振幅響應(yīng)')xlabel('頻率(單位:pi)');ylabel('Hr(w)')set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])set(gca,'YTickMode','manual','YTick',[0,1]);gridsubplot(2,2,4);plot(w/pi,db);axis([0,1,-200,10]);gridtitle('幅度響應(yīng)');xlabel('頻率(單位:pi)');ylabel('分貝數(shù)');set(gca,'XTickMode','Manual','XTick',[0;0.2;0.3;1]);set(gca,'YTickMode','Manual','YTick',[-16;0]);set(gca,'YTickLabelMode','manual','YTickLabels',['16';'0'])%function[db,mag,pha,grd,w]=freqz_m(b,a);%------------------------------------%[db,mag,pha,grd,w]=freqz_m(b,a);%db：0-π區(qū)間內(nèi)的相對振幅(dB)%mag：0-π區(qū)間內(nèi)的絕對振幅%pha：0-π區(qū)間內(nèi)的相移%grd：0-π區(qū)間內(nèi)的群延遲%w：0-π區(qū)間內(nèi)的501個抽樣點(diǎn)頻率%b：系統(tǒng)函數(shù)H(z)中分子多項(xiàng)式系數(shù)(對FIRb=h)%a：系統(tǒng)函數(shù)H(z)中分母多項(xiàng)式系數(shù)(對FIR:a=[1])%[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);grd=grpdelay(b,a,w);擴(kuò)展函數(shù)Hr_Type2的MATLAB程序?yàn)閒unction[Hr,w,b,L]=Hr_Type2(h);%-----------------------------------------------------------%[Hr,w,b,L]=Hr_Type2(h)%Hr:頻響振幅%w：0-π區(qū)間內(nèi)計(jì)算Hr的頻率點(diǎn)%b：低通濾波器的系數(shù)%L：Hr的階次%h：低通濾波器的單位抽樣響應(yīng)%M=length(h);L=M/2;b=2*[h(L:-1:1)];n=[1:1:L];n=n-0.5;w=[0:1:500]'*pi/500;Hr=cos(w*n)*b';例8.28%例8.28求解三種結(jié)構(gòu)形式及轉(zhuǎn)換的濾波系數(shù)MATLAB程序%b=[1,-3,11,-27,18];a=[16,12,2,-4,-1];[b0,B,A]=dir2cas(b,a);[C,B1,A1]=dir2par(b,a);N=24;n=1:N+1;formatlong;delta=impseq(0,0,N);h1=filter(b,a,delta);h2=casfiltr(b0,B,A,delta);h3=parfiltr(C,B1,A1,delta);figure(1)subplot(3,1,1)plot(n,h1)title('直接型')subplot(3,1,2)plot(n,h2)title('級聯(lián)型')subplot(3,1,3)plot(n,h3)title('并聯(lián)型')function[b,a]=dir2cas(b0,B,A);%將直接型結(jié)構(gòu)變成級聯(lián)型結(jié)構(gòu),求系統(tǒng)函數(shù)的系數(shù)%b0為常數(shù);B、A分別為級聯(lián)型系統(tǒng)函數(shù)的分子、分母%b、a分別為直接型系統(tǒng)函數(shù)的分子、分母%計(jì)算增益系數(shù)b0b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;%補(bǔ)零,使b和a等長M=length(b);N=length(a);ifN>Mb=[b,zeros(1,N-M)];elseifM>Na=[a,zeros(1,M-N)];N=M;elseNM=0;end%計(jì)算多項(xiàng)式的根K=floor(N/2);B=zeros(K,3);A=zeros(K,3);ifK*2==Nb=[b0];a=[a,0];end%用函數(shù)cplxpair把根以共軛復(fù)根對的次序排列,并用poly函數(shù)轉(zhuǎn)換成二階多項(xiàng)式broots=cplxpair(roots(b));aroots=cplxpair(roots(a));fori=1:2:2*KBrow=broots(i:1:i+1,:);Brow=real(poly(Brow));B(fix((i+1)/2),:)=Brow;Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow));A(fix((i+1)/2),:)=Arow;end%直接型到并聯(lián)型的轉(zhuǎn)換%--------------------------------------%function[C,B,A]=dir2par(b,a)%C=當(dāng)length(b)>=length(a)時的多項(xiàng)式部分%B=包含各系數(shù)bk的實(shí)系數(shù)矩陣%A=包含各ar系數(shù)的實(shí)系數(shù)矩陣%b=直接型的分子多項(xiàng)式系數(shù)%a=直接型的分母多項(xiàng)式系數(shù)function[C,B,A]=dir2par(b,a)M=length(b);N=length(a);[r1,p1,C]=residuez(b,a);p=cplxpair(p1,10000000*eps);I=cplxcomp(p1,p);r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);ifK*2==N;%Neven,orderofA(z)odd,onefactorisfirstorderfori=1:2:N-2Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix((i+1)/2),:)=real(Brow);A(fix((i+1)/2),:)=real(Arow);end[Brow,Arow]=residuez(r(N-1),p(N-1),[]);B(K,:)=[real(Brow)0];A(K,:)=[real(Arow)0];elsefori=1:2:N-1Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix((i+1)/2),:)=real(Brow);A(fix((i+1)/2),:)=real(Arow);endendfunctiony=casfiltr(b0,B,A,x);%IIR%-----------------------------------------------%y=casfiltr(b0,B,A,x);%y=輸出序列%b0=級聯(lián)型的增益系數(shù)%B=包含各系數(shù)bk的實(shí)系數(shù)矩陣%A=包含各ar系數(shù)的實(shí)系數(shù)矩陣%x=輸入序列%[K,L]=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=x;fori=1:1:Kw(i+1,:)=filter(B(i,:),A(i,:),w(i,:));endy=b0*w(K+1,:);functiony=parfiltr(C,B,A,x);%IIR濾波器的并聯(lián)型實(shí)現(xiàn)%----------------------------------------%[y]=parfiltr(C,B,A,x);%y=輸出序列%C=當(dāng)M>=N時(FIR)的多項(xiàng)式部分%B=包含各bk系數(shù)的實(shí)系數(shù)矩陣%A=包含各ar系數(shù)的實(shí)系數(shù)矩陣%x=輸入序列%[K,L]=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=filter(C,1,x);fori=1:1:Kw(i+1,:)=filter(B(i,:),A(i,:),x);endy=sum(w);function[x,n]=impseq(n0,n1,n2)%產(chǎn)生單位抽樣序列n=[n1:n2];x=[(n-n0)==0];functionI=cplxcomp(p1,p2)%比較兩個包含同樣標(biāo)量元素,但可能有不同下標(biāo)的復(fù)數(shù)對,返回其數(shù)組下標(biāo)%這一函數(shù)必須用在cplxpair函數(shù)之后,以便對頻率極點(diǎn)矢量及其響應(yīng)的留數(shù)矢量重新排序p2=cplxpair(p1)I=[];forj=1:1:length(p2)fori=1:1:length(p1)if(abs((p1(i))-p2(j))<0.0001)I=[I,i];endendendI=I';例8.29%例8.29中計(jì)算FIR濾波器級聯(lián)型結(jié)構(gòu)的MATLAB程序b=[1,0,0,0,16+1/16,0,0,0,1];[b0,B,A]=dir2cas(b,1)例9.6%例9.6中周期圖法計(jì)算信號功率譜的Matlab程序clfFs=2000;%情況1:數(shù)據(jù)長度N1=256N1=256;N1fft=256;n1=0:N1-1;t1=n1/Fs;f1=50;f2=100;xn1=sin(2*pi*f1*t1)+2*sin(2*pi*f2*t1)+randn(1,N1);Pxx1=10*log10(abs(fft(xn1,N1fft).^2)/N1);f1=(0:length(Pxx1)-1)*Fs/length(Pxx1);subplot(2,1,1)plot(f1,Pxx1)ylabel(′功率譜(dB)′);title(′數(shù)據(jù)長度N1=256′)grid%情況2:數(shù)據(jù)長度N2=1024N2=1024;N2fft=1024;n2=0:N2-1;t2=n2/Fs;f1=50;f2=100;xn2=sin(2*pi*f1*t2)+2*sin(2*pi*f2*t2)+randn(1,N2);Pxx2=10*log10(abs(fft(xn2,N2fft).^2)/N2);f2=(0:length(Pxx2)-1)*Fs/length(Pxx2);subplot(2,1,2)plot(f2,Pxx2)xlabel(‘頻率(Hz)’);ylabel(‘功率譜(dB)’);title(‘?dāng)?shù)據(jù)長度N2=1024’)例9.8%基于功率譜估計(jì)的系統(tǒng)辨識%產(chǎn)生輸入信號N=1024;nFFT=256;window=hanning(256);noverlap=128;dflag=′none′;Fs=1000;n=0:N-1;t=n/Fs;randn(′state′,0);xn=sin(2*pi*50*t)+randn(1,N);%構(gòu)建系統(tǒng)：一個濾波器特性h=ones(1,10)/10;%計(jì)算系統(tǒng)輸出yn=filter(h,1,xn);%進(jìn)行系統(tǒng)函數(shù)估計(jì)[Txy,f]=tfe(xn,yn,nFFT,Fs,window,noverlap,dflag);%計(jì)算系統(tǒng)頻率響應(yīng)H=freqz(h,1,f,Fs);%比較響應(yīng)結(jié)果subplot(2,1,1);plot(f,abs(H));xlabel(′頻率（Hz）′);ylabel(′幅度′);title(′頻率響應(yīng)′);axis([050001]);gridsubplot(2,1,2);plot(f,abs(Txy));xlabel(′頻率（Hz）′);ylabel(′幅度′);title(′頻率響應(yīng)′);axis([050001])grid例9.9%example9.9MatlabProgramN=1000;n=0:N-1;Fs=500;t=n/Fs;Lag=100;%相關(guān)信號的最大延遲量x1=sin(2*pi*10*t)+0.6*randn(1,length(t));%含白噪聲的正弦信號x1[c,lags]=xcorr(x1,Lag,′unbiased′);%無偏自相關(guān)函數(shù)的計(jì)算subplot(2,2,1)%畫出x1曲線plot(t,x1);xlabel(′t′);ylabel(′x1(t)′);title(′含白噪聲的正弦信號x1′);gridsubplot(2,2,2);%畫x1自相關(guān)曲線plot(lags/Fs,c);xlabel(′t′);ylabel(′Rxx1(t)′);title(′x1的自相關(guān)函數(shù)′);gridx2=randn(1,length(t));%發(fā)生白噪聲x2\[c,lags\]=xcorr(x2,Lag,′unbiased′);%白噪聲x2的無偏自相關(guān)函數(shù)subplot(2,2,3)%畫x2曲線plot(t,x2);xlabel(′t′);ylabel(′x2(t)′);title(′白噪聲x2′);gridsubplot(2,2,4)%畫x2自相關(guān)曲線plot(lags/Fs,c)xlabel(′t′);ylabel(′Rxx2(t)′);title(′x2的自相關(guān)曲線′);grid例9.10%采用白噪聲輸入，進(jìn)行系統(tǒng)辨識h=fir1(30,0.2,boxcar(31));x=randn(16384,1);y=filter(h,1,x);tfe(x,y,1024,[],[],512)例10.3t=0:0.01:2.10;%抽樣頻率為100Hzs1=sin(2*pi*30*t);%產(chǎn)生頻率為30Hz的正弦信號s2=s1+0.5*[zeros(1,50)s1(1:161)];%加上回聲，其幅度是原信號的一半%時域上延遲50個抽樣周期即0.5秒c=cceps(s2);%用函數(shù)cceps求出倒譜subplot(121),plot(t,s2)%畫出倒譜圖xlabel('\fontsize{14}t')ylabel('\fontsize{14}\fontname{courier}信號+噪聲的幅值')subplot(122),plot(t,c)xlabel('\fontsize{14}t')ylabel('\fontsize{14}\fontname{courier}倒譜')例10.4%%例10.4，利用FFT對信號進(jìn)行STFT分析A=2;w0=10*pi*0.00001;%產(chǎn)生一個被分析的信號M=20000;n=0:M-1;x=A*sin(n.^2*w0);%%進(jìn)行短時傅里葉變換L=200;P=L/2;N=256;Fs=10000;%%漢寧窗的窗口長度L，做FFT運(yùn)算w=0.5*(1-cos(2*pi*(0:L-1)/(L-1)));Q=fix((M-P)/(L-P));forq=0:Q-1x0=x(q*(L-P)+1:q*(L-P)+L).*w;X(q+1,:)=fft(x0,N);endtn=((0:Q-1)*(L-P))/Fs;fk=(0:N/2)*Fs/N;zhX=X