信号处理引论第五章及matlab作业答案

第五章习题与上机题5.1 已知序列12()(),0 1 , ()()()nx n a u n a x n u n u n N =<<=--，分别求它们的自相关函数，并证明二者都是偶对称的实序列。

解：111()()()()()nn mx n n r m x n x n m a u n au n m ∞∞-=-∞=-∞=-=-∑∑当0m ≥时，122()1mmnx n ma r m aaa∞-===-∑ 当0m <时，122()1m mnx n a r m aaa -∞-===-∑ 所以，12()1mx ar m a =-2 ()()()()N x n u n u n N R n =--=22210121()()()()()1,0 =1,00, =()(1)x NN n n N mn N n m N r m x n x n m Rn R n m N m N m N m m Nm N m R m N ∞∞=-∞=-∞--=-=-=-=-⎧=--<<⎪⎪⎪⎪=-≤<⎨⎪⎪⎪⎪⎩-+-∑∑∑∑其他从1()x r m 和2()x r m 的表达式可以看出二者都是偶对称的实序列。

5.2 设()e()nTx n u n -=，T 为采样间隔。

求()x n 的自相关函数()x r m 。

解：解：()()()()e()e ()nTn m T x n n r m x n x n m u n u n m ∞∞---=-∞=-∞=-=-∑∑用5.1题计算1()x r m 的相同方法可得2e()1e m Tx Tr m --=-5.3 已知12()sin(2)sin(2)s s x n A f nT B f nT ππ=+，其中12,,,A B f f 均为常数。

求()x n 的自相关函数()x r m 。

解：解：()x n 可表为)()()(n v n u n x +=的形式，其中)2sin()(11s nT f A n u π=，=)(n v 22sin(2)s A f nT π，)(),(n v n u 的周期分别为 s T f N 111=，sT f N 221=，()x n 的周期N 则是21,N N 的最小公倍数。

MATLAB作业5参考答案1、 试求出下面线性微分方程的通解。

543225432()()()()()136415217680()[sin(2)cos(3)]3t d y t d y t d y t d y t dy t y t e t t dt dt dt dt dt π-+++++=++假设上述微分方程满足已知条件(0)1,(1)3,()2,(0)1,(1)2y y y y y π=====，试求出满足该条件的微分方程的解析解。

【求解】先定义t 为符号变量，求出等号右侧的函数，则可以由下面命令求出方程的解析 解，解的规模较大，经常能占数页。

>> syms texp(-2*t)*(sin(2*t+sym(pi)/3)+cos(3*t))ans =exp(-2*t)*(sin(2*t+1/3*pi)+cos(3*t))>> y=dsolve(['D5y+13*D4y+64*D3y+152*D2y+176*Dy+80*y=',...'exp(-2*t)*(sin(2*t+1/3*pi)+cos(3*t))'],'y(0)=1','y(1)=3','y(pi)=2',...'Dy(0)=1','Dy(1)=2')略：事实上，仔细阅读求出的解析解就会发现，其中大部分表达式是关于系数的，所以如果能对 系数进行近似则将大大减小解的复杂度。

>> vpa(y)ans =.20576131687242798353909465020576e-2*exp(-2.*t)*cos(3.*t)+.15538705805619602372728107411086e-1*exp(-2.*t)*sin(2.*t)+.76830587084294035590921611166287e-2*exp(-2.*t)*cos(2.*t)-106.24422608844727797303237726774*exp(-2.*t)*t^2+98.159206062620455331994871615083*exp(-2.*t)*t+59.405044899367325888329709780356*exp(-2.*t)*t^3-30.741892776456442808809983330755*exp(-2.*t)+.20576131687242798353909465020576e-2*exp(-2.*t)*sin(3.*t)+31.732152104579289125415500223136*exp(-5.*t)2、 试求解下面微分方程的通解以及满足(0)1,()2,(0)0x x y π===条件下的解析解。

light('position',[2,1,1])
colormap([jet; flipud(jet)])
disp('按任意键，观察色图变幻。退出按Ctrl+C')
pause
spinmap(80,9)
%习题5_10p
function f=prob5_10(K,ki)
%prob5_10函数产生动态衰减正弦函数，K控制动态曲线动态变化的循环次数，ki控制曲线动态变化的快慢
% size(find(isnan(Z)))
% sum(sum(isnan(Z)))
%习题5_8
ezplot('y/(1+x^2+y^2)-0.1',[-2*pi,2*pi,-pi/3,3.5*pi])
hold on
ezplot('sin(x+cos(y))',[-2*pi,2*pi,-pi/3,3.5*pi])
mesh(X,Y,Z)
hidden off
% colormap(cool),
% shading interp,
syms x y z
% z=4.*x.*exp(-x.^2-y.^2);
z=4*x*exp(-x^2-y^2);
ezmesh(z,[-3,3])
hidden off
%习题5_7
hold on
for k=1:L
zk=zeta(k);
beta=sqrt(abs(1-zk^2));
if zk<1 缺陷在此，由于计算机的精度，zeta(5)<1 ,可改为zk-1<-2*eps %满足此条件，绘蓝色线

数字信号处理课后Matlab习题解答学院：电子信息工程学院专业：通信工程指导老师：钱满义2015年6月26日目录课后Matlab习题解答 (1)1、第一章 (4)M1-1 (4)M1-2 (5)M1-3 (7)M1-4 (9)M1-5 (10)M1-6 (11)2、第二章 (12)M2-1 (12)M2-2 (15)M2-3 (16)M2-4 (18)M2-5: (20)M2-6 (21)3、第四章 (23)M4-1 (23)M4-2 (24)M4-3 (25)M4-4 (27)M4-5 (28)M4-6 (29)M4-7 (31)M4-8 (32)M4-9 (34)4、第五章 (35)M5-1 (35)M5-2 (35)M5-3 (37)M5-4 (39)M5-5 (39)M5-6 (39)M5-7 (40)M5-8 (40)M5-9 (40)M5-10 (40)5、第六章 (40)M6-1 (40)M6-2 (40)6、第七章 (44)M7-1 (44)M7-2 (45)M7-3 (47)M7-4 (48)M7-5 (48)1、第一章M1-1代码：f=10;t=-0.2:0.001:0.2;g1=cos(6.*pi.*t);g2=cos(14.*pi.*t);g3=cos(26.*pi.*t);k=-0.2:1/f:0.2;h1=cos(6.*pi.*k);h2=cos(14.*pi.*k);h3=cos(26.*pi.*k);subplot(4,1,1);plot(k,h1,'r.',t,g1,'r');xlabel('t');ylabel('g1(t)');subplot(4,1,2);plot(k,h2,'g.',t,g2,'g');xlabel('t');ylabel('g2(t)');subplot(4,1,3);plot(k,h3,'b.',t,g3,'b');xlabel('t');ylabel('g3(t)');subplot(4,1,4);plot(t,g1,'r',t,g2,'g',t,g3,'b',k,h1,'r.',k,h2,'g.',k,h3,'b.') xlabel('t');ylabel('g(t)');legend('g1(t)','g2(t)','g3(t)');图像：M1-2(1)代码：b=[1];a=[1 -1.845 0.850586];w=linspace(0,2*pi,200);y=filter(b,a,w);plot(w,y);ylabel('单位脉冲响应'); xlabel('Normalized frequency');%b=[1];%a=[1 -1.845 0.850586];h=impz(b,a,21);%figure(1);stem(0:20,h);ylabel('单位脉冲响应'); xlabel('Normalized frequency');图像：(2)代码：b=[1];a=[1 -1.85 0.85];w=linspace(0,2*pi,200);y=filter(b,a,w);plot(w,y);ylabel('单位脉冲响应'); xlabel('Normalized frequency'); b=[1];a=[1 -1.85 0.85];h=impz(b,a,21);figure(1);stem(0:20,h);ylabel('单位脉冲响应');图像：M1-3(1)代码：x=firls(511,[0 0.4 0.404 1],[1 1 0 0]); k=linspace(0,1,512);stem(k,x);图像：(2)代码：x=firls(511,[0 0.4 0.404 1],[1 1 0 0]); k=linspace(0,1,512);stem(k,x);k1=0.4*pi;y1=x.*cos(k1*k); A=abs(y1); subplot(4,1,1); plot(A);k2=0.8*pi;y2=x.*cos(k2*k); B=abs(y2); subplot(4,1,2); plot(B);k3=0.9*pi;y3=x.*cos(k3*k); C=abs(y3); subplot(4,1,3); plot(C);k4=pi;y4=x.*cos(k4*k); D=abs(y4); subplot(4,1,4); plot(D);图像：M1-4代码：b=[1 1];a=2;w=linspace(0,pi,512);H1=freqz(b,a,w);plot(w/pi,unwrap(abs(H1)),'r'); hold onb1=[1.8];a1=[1 -0.8];w=linspace(0,pi,512);H2=freqz(b1,a1,w);plot(w/pi,unwrap(abs(H2)),'g'); hold onb2=conv(b,b1);a2=conv(a,a1);w=linspace(0,pi,512);H3=freqz(b2,a2,w);plot(w/pi,unwrap(abs(H3),'b')); ylabel('幅度');xlabel('Normalized frequency'); 图像：M1-5(1)代码：b1=[0.0534 0.0534];b2=[1 2.0332 2.0169 -1.0166 1]; a1=[1 -0.683];a2=[1 -1.4661 0.7957];b=conv(b1,b2);a=conv(a1,a2);w=linspace(0,pi,200);H=freqz(b,a,w);subplot(2,1,1);plot(w/pi,abs(H));ylabel('幅度');xlabel('Normalized frequency'); subplot(2,1,2);plot(w,angle(H));ylabel('相位');xlabel('Normalized frequency'); 图像：(2)代码：b1=[1 -2 1];b2=[1 -2 1];a1=[1 -1.499 0.8482];a2=[1 -1.5548 0.6493];b=conv(b1,b2);a=conv(a1,a2);w=linspace(0,pi,200);H=freqz(b,a,w);subplot(2,1,1);plot(w/pi,abs(H));ylabel('幅度');xlabel('Normalized frequency'); subplot(2,1,2);plot(w,angle(H));ylabel('相位');xlabel('Normalized frequency'); 图像：M1-6代码：b=[1 2 0.99];a=[1 1.55 0.6];w=linspace(0,pi,512);H=freqz(b,a,w);plot(w,unwrap(angle(H)),'r'); hold on;b1=[1.1 1.99 0.9];a1=[1 1.55 0.6];w=linspace(0,pi,512);Hmin=freqz(b1,a1,w);plot(w,unwrap(angle(Hmin)),'g'); hold on;b2=[0.9 1.99 1.1];a2=[1 1.55 0.6];w=linspace(0,pi,512);Hmax=freqz(b2,a2,w);plot(w,unwrap(angle(Hmax)),'m'); ylabel('相位');xlabel('omega');图像：2、第二章M2-1（1）代码：g=[1 -3 4 2 0 -2];h=[3 0 1 -1 2 1];l=length(g);GE=fft(g,L);HE=fft(h,L);y1=ifft(GE.*HE);for n=1:lifn+l<=Ly2(n)=y1(n)+y1(n+l);elsey2(n)=y1(n);endendy2stem(0:l-1,y2)xlabel('k')ylabel('y(k)')title('循环卷积')结果：y2 =6.0000 -3.0000 17.0000 -2.0000 7.0000 -13.0000 图像：（2）代码：k=0:5;x=cos(pi.*k./2);y=3.^k;L=2*l-1;GE=fft(x,L);HE=fft(y,L);y1=ifft(GE.*HE);for n=1:lifn+l<=Ly2(n)=y1(n)+y1(n+l);elsey2(n)=y1(n);endendy2stem(0:l-1,y2)xlabel('k')ylabel('y’(k)')title('循环卷积')结果：y2 =-71.0000 -213.0000 89.0000 267.0000 73.0000 219.0000 图像：M2-2（1）代码：N=10;k=-N:N;x=cos(k.*pi./(2*N));W=linspace(-pi,pi,512);X=zeros(1,length(W));for k=-N:NX1=x(k+N+1).*exp(-j.*W.*k);X=X+X1;endplot(W,abs(X))xlabel('W');ylabel('abs(X)');图像：（2）代码：N=10;k=-N:N;x=cos(k.*pi./(2*N));X_21=fft(x,21);L=-10:10;W=linspace(-pi,pi,1024);X=zeros(1,length(W));for k=-N:NX1=x(k+N+1).*exp(-j.*W.*k);X=X+X1;endplot(W,abs(X));hold on;plot(2*pi*L/21,fftshift(abs(X_21)),'o');xlabel('W');ylabel('abs(X)');图像：M2-3代码：N=64;L=1024;f1=100;f2=120;fs=800;A=1;B1=1;B2=0.5;B3=0.25;B4=0.05;T=1/fs;ws=2*pi*fs;k=0:N-1;x1=A*cos(2*pi*f1*T*k)+B1*cos(2*pi*f2*T*k); x2=A*cos(2*pi*f1*T*k)+B2*cos(2*pi*f2*T*k); x3=A*cos(2*pi*f1*T*k)+B3*cos(2*pi*f2*T*k); x4=A*cos(2*pi*f1*T*k)+B4*cos(2*pi*f2*T*k); hf=(hamming(N))';x1=x1.*hf;x2=x2.*hf;x3=x3.*hf;x4=x4.*hf;X1=fftshift(fft(x1,L));X2=fftshift(fft(x2,L));X3=fftshift(fft(x3,L));X4=fftshift(fft(x4,L));W=T*(-ws/2+(0:L-1)*ws/L)/(2*pi); subplot(2,2,1);plot(W,abs(X1));title('A=1,B=1');xlabel('W');ylabel('X1');subplot(2,2,2);图像：F2=120HZ时：F2=140HZ时：F2=160HZ时：M2-4（1）代码：W0=2*pi/15;W1=2.3*pi/15;N=64;k=0:N-1;x=cos(W0*k)+0.75*cos(W1*k); X=fft(x);plot(k/N,abs(X));grid on;title('64点FFT');图像：（2）代码：W0=2*pi/15;W1=2.3*pi/15;N=64;L=1024;k=0:N-1;x=cos(W0*k)+0.75*cos(W1*k); X=fft(x,L);plot((0:L-1)/N,abs(X));grid on;title('1024点FFT');图像：M2-5:代码：fs=100;ws=2*pi*fs;Ts=1/fs;N=fs;x=exp(-3*Ts*(0:N-1));y=fft(x,N);l=length(y);k=linspace(-ws/2,ws/2,l); plot(k,Ts*fftshift(abs(y)),'b:'); hold on;w=linspace(-ws/2,ws/2,1024); y1=sqrt(1./(9+w.^2));plot(w,y1,'r')title('fs=100Hz时的频谱') legend('近似值','实际值');图像：fs=10HZ时：fs=50Hz时：fs=100HZ时:M2-6 代码：Ts=0.5;N=4;N0=64;k=(-N/2:(N/2))*Ts;x=exp(-pi*(k).^2);X=Ts*fftshift(fft(x,N0));w=-pi/Ts:2*pi/N0/Ts:(pi-2*pi/N0)/Ts;XT=(pi/pi)^0.5*exp(-w.^2/4/pi);subplot(2,1,1)plot(w/pi,abs(X),'-o',w/pi,XT);xlabel('\omega/\pi');ylabel('X(j\omega)');legend('试验值','理论值');title(['Ts=',num2str(Ts) ' ''N=',num2str(N)]); subplot(2,1,2)plot(w/pi,abs(X)-XT)ylabel('实验误差')xlabel('\omega/\pi');图像：Ts=1,N=2:Ts=0.5 ,N=2:Ts=0.5 ,N=4:3、第四章M4-1代码：WP=10;WS=2;Ap=1;As=40;wp=1/WP;ws=1/WS;w0=1;B=2;[N,Wc]=cheb1ord(wp,ws,Ap,As,'s'); [num,den] = cheby1(N,Ap,Wc,'s'); [numt,dent] = lp2hp(num,den,1);w=linspace(1,12,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h))) ; grid ; xlabel('Frequency in rad/s');ylabel('Gain in dB')图像：M4-2代码：wp=1;ws=3.3182;Ap=1;As=32;w0=sqrt(48);B=2;[N,Wc]=buttord(wp,ws,Ap,As,'s'); [num,den] = butter(N,Wc,'s'); [numt,dent] = lp2bp(num,den,w0,B); w=linspace(2,12,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h))) ; grid ; xlabel('Frequency in rad/s');ylabel('Gain in dB')图像：M4-3（1）代码：wp=1;ws=3.3182;Ap=1;As=32;w0=sqrt(48);B=2;[N,Wc]=cheb1ord(wp,ws,Ap,As,'s'); [num,den] = cheby1(N,Ap,Wc,'s'); [numt,dent] = lp2bp(num,den,w0,B); w=linspace(2,12,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h))) ; grid ; xlabel('Frequency in rad/s');ylabel('Gain in dB')图像：（2）代码：wp=1;ws=3.3182;Ap=1;As=32;w0=sqrt(48);B=2;[N,Wc]=cheb2ord(wp,ws,Ap,As,'s'); [num,den] = cheby2(N,Ap,Wc,'s'); [numt,dent] = lp2bp(num,den,w0,B); w=linspace(2,12,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h))) ; grid ; xlabel('Frequency in rad/s');ylabel('Gain in dB')图像：（3）代码：wp=1;ws=3.3182;Ap=1;As=32;w0=sqrt(48);B=2;[N,Wc]=ellipord(wp,ws,Ap,As,'s'); [num,den] = ellip(N,Ap,As,Wc,'s'); [numt,dent] = lp2bp(num,den,w0,B); w=linspace(2,12,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h))) ; grid ; xlabel('Frequency in rad/s');ylabel('Gain in dB')图像：M4-4代码：Ap=1;As=10;wp1=6;wp2=13;ws1=9;ws2=11;B=ws2-ws1;w0=sqrt(ws1*ws2);wLp1=B*wp1/(w0*w0-wp1*wp1); wLp2=B*wp2/(w0*w0-wp2*wp2);[N,Wc]=buttord(wLp,1,Ap,As,'s'); [num,den] = butter(N,Wc,'s' ); [numt,dent]=lp2bs(num,den,w0,B); w=linspace(5,14,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h)));w=[wp1 ws1 ws2 wp2];set(gca,'xtick',w);grid;h=freqs(numt,dent,w);A=-20*log10(abs(h));图像：M4-5代码：Ap=1;As=10;wp1=6;wp2=13;ws1=9;ws2=11;B=ws2-ws1;w0=sqrt(ws1*ws2);wLp1=B*wp1/(w0*w0-wp1*wp1); wLp2=B*wp2/(w0*w0-wp2*wp2);[N,Wc]=ellipord(wLp,1,Ap,As,'s'); [num,den] = ellip(N,Ap,As,Wc,'s'); [numt,dent] = lp2bp(num,den,w0,B); w=linspace(5,14,1000);h=freqs(numt,dent,w);plot(w,20*log10(abs(h)));w=[wp1 ws1 ws2 wp2];set(gca,'xtick',w);grid;h=freqs(numt,dent,w);A=-20*log10(abs(h));图像：M4-6代码：Wp=0.1*pi;Ws=0.4*pi;Ap=1;As=25;Fs=1wp=Wp*Fs;ws=Ws*Fs;N=buttord(wp,ws,Ap,As,'s');wc=wp/(10^(0.1*Ap)-1)^(1/2/N); [numa,dena]=butter(N,wc,'s'); [numd,dend]=impinvar(numa,dena,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h));numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm));w=[WpWs];h=freqz(numd,dend,w);fprintf('Ap= %.4f\n',-20*log10( abs(h(1)))); fprintf('As= %.4f\n',-20*log10( abs(h(2))));图像：代码：WP=0.1*pi;WS=0.4*pi;Ap=1;As=25;Fs=1;wp=WP*Fs;ws=WS*Fs;N=cheb2ord(wp,ws,Ap,As,'s');wc=wp/(10^(0.1*Ap)-1)^(1/2/N); [numa,dena]=cheby2(N,As,wc,'s'); [numd,dend]=impinvar(numa,dena,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h));numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm))w=[wpws];fprintf('Ap= %.4f\n',-20*log10( abs(h(1)))); fprintf('As= %.4f\n',-20*log10( abs(h(2))));图像：M4-7代码:Wp=0.2*pi;Ws=0.4*pi;Ap=1;As=15;Fs=1;wp=Wp*Fs;ws=Ws*Fs;[N,wc]=buttord(wp,ws,Ap,As,'s'); [numa,dena]=butter(N,wc,'s');[num,dend]=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(num,dend,w);norm=max(abs(h));plot(w/pi,20*log10(abs(h)/norm)); xlabel('Normalized frequency');ylabel('Gain,dB');grid on;w=[wpws];numd=num/norm;图像：M4-8(1)脉冲响应不变法代码：fp=2000;fs=10000;Ap=0.5;As=50;Fs=44100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;[N,wc]=buttord(wp,ws,Ap,As,'s'); [numa,dena]=butter(N,wc,'s'); [num,dend]=impinvar(numa,dena,Fs); w=linspace(0,2*pi,1024);h=freqz(num,dend,w);norm=max(abs(h));plot(w/pi,20*log10(abs(h)/norm)); xlabel('Normalized frequency'); ylabel('Gain,dB');grid on;图像：（2）双线性变换法代码：fp=2000;fs=10000;Ap=0.5;As=50;Fs=20000;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;[N,wc]=buttord(wp,ws,Ap,As,'s'); [num,den]=butter(N,wc,'high'); [numc,dend]=bilinear(num,den,Fs); w=linspace(0,2*pi,1024);h=freqz(numc,dend,w);norm=max(abs(h));plot(w/pi,20*log10(abs(h)/norm)); xlabel('Normalized frequency'); ylabel('Gain,dB');grid on;图像：M4-9代码：wp=0.8*pi;ws=0.6*pi;Ap=0.5;As=30;[N,wc]=buttord(wp,ws,Ap,As,'s'); [num,den]=butter(N,wc,'high'); [numc,dend]=bilinear(num,den,Fs); w=linspace(0,pi,2048);h=freqz(numc,dend,w);norm=max(abs(h));plot(w/pi,20*log10(abs(h)/norm)); xlabel('Normalized frequency'); ylabel('Gain,dB');grid on;图像：4、第五章M5-1M5-2(1)求增益响应代码：Wp=0.6*pi;Ws=0.4*pi;As=45;Ap=1;N=ceil(7*pi/(Wp-Ws));N=mod(N+1,2)+N;M=N-1;w=hamming(N)';Wc=(Wp+Ws)/2;k=0:M;hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi); hd(0.5*M+1)=hd(0.5*M+1)+1;h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,[1],omega);magdb=20*log10(abs(mag));plot(omega/pi,magdb);xlabel('Normalized frequency'); ylabel('Gain, db');grid;图像：(2)代码：Wp=0.6*pi;Ws=0.4*pi;As=45;Ap=1;N=ceil(7*pi/(Wp-Ws));N=mod(N+1,2)+N;M=N-1;w=hamming(N)';Wc=(Wp+Ws)/2;k=0:M;hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);hd(0.5*M+1)=hd(0.5*M+1)+1;h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag));%plot(omega/pi,magdb);xlabel('Normalized frequency');ylabel('Gain, db');grid;w=[WpWs];H=freqz(h,[1],w);fprintf('Ap=%.4f\n',-20*log10(abs(H(1)))); fprintf('As=%.4f\n',-20*log10(abs(H(2))));结果：Ap=0.0182As=53.5634M5-3(1)blackman窗：代码：Wp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=45;N=ceil(11.4*pi/(Ws-Wp));N=mod(N+1,2)+N;M=N-1;w=blackman(N)';wc=(Ws+Wp)/2;k=0:M;hd=-(wc/pi)*sinc(wc*(k-0.5*M)/pi); h=hd.*w;n=linspace(0,pi,512);mag=freqz(h,[1],n);magdb=20*log10(abs(mag));plot(n/pi,magdb);xlabel('normalized frequency'); ylabel('Gain(dB)');title('blackman');grid图像：(2)kaizer窗（a为通，阻带。

为得到 H ( z ) ，
(1) 由极点构成 H a ( s ) 的分母多项式，分子为分母多项式的常数。 (2) H a ( s ) 展成部分分式。 (3) 据有理分式变换得到对应的 H ( z ) 各分式，整理得到最后的 H ( z ) 。 22、 取 T=1， 预畸， 由已知列出对模拟滤波器的衰减要求， 解出 N=6.04， 取 N=7， 得到
−0.5
Z −1
−1
0.9
−0.81
4、 H ( z ) = −4.9383 +
2.1572 4.7811 − 1.5959 z −1 + 1 + 0.5 z −1 1 − 0.9 z −1 + 0.81z −2
−4.9383
x ( n) y ( n)
2.1572 −0.5
Z −1
4.7811
Z
0.9 −0.81
= H 2 ( z)
α 02 + α12 z -1 -3.1986 + 0.2591z -1 = 1 +z 2 1 + 1.618 z - 4π 2 2 1 + r z 1 - 2rz -cos 5
频率取样型实现流程图：
−10.125
Z −1
18.3236
x ( n)
Z −1
x ( n)
Z −1
Z −1
+
Z −1
− 7 4
+
Z −1
− 69 8
+
y ( n) 4） 频率取样型：取 r=1，N=5，得到 DFT{h(n)}为：
{-10.1250 9.1618 + 6.6564i -1.5993 - 4.9221i -1.5993 + 4.9221i 9.1618 - 6.6564i}

FIR 数字滤波器设计本章知识点：对于一个离散时间系统∑∑=-=--=M 1n nn 1-N 0n nnz a 1z bz H )(,若分母多项式中系数0a a a M 21==== ，则此系统就变成一个FIR 系统∑-=-=1N 0n n nz bz H )(，其中系数1-N 10b ,.b ,b 即为该系统的单位取样响应h ( 0 ) , h ( 1 ) ,… h ( N-1 ),且当n > N-1时，h ( n ) = 0。

FIR 系统函数H(z) 在Z 平面上有N-1个零点，在原点z=0处有N-1个重极点。

这类系统不容易取得较好的通带和阻带特性，要想得到与IIR 系统类似的衰减特性，则要求较高的H(z)阶次。

相比于IIR 系统来说，FIR 系统主要有三大突出优点：1）系统永远稳定；2）易于实现线性相位系统；3）易于实现多通带（或多组带）系统。

线性相位FIR 滤波器实现的充要条件是：对于任意给定的数值N （奇数或偶数），冲激响应h[n] 相对其中心轴21-N 必须成偶对称或奇对称，此时滤波器的相位特性是线性的，且群延时均为常数 21-=N τ。

由于h(n) 有奇对称和偶对称两种情况，h(n)的点数N 有奇数、偶数之分。

因此，h （n ）可以有4种不同的类型，分别对应于4种线性相位FIR 数字滤波器：h[n] 偶对称N 为奇数、h[n] 偶对称N 为偶数、h[n] 奇对称N 为奇数、h[n] 奇对称N 为偶数。

四种线性相位FIR 滤波器的特性归纳对比于表5.１中。

一．FIR DF 设计方法FIR DF 的设计实现不能像IIR DF 设计那样借助于模拟滤波器的设计方法来实现，其设计方法主要是建立在对理想滤波器频率特性进行不同程度逼近的基础上，主要的逼近方法有三种：窗函数法；频率抽样法；最佳一致逼近法。

1． 窗函数法窗函数法是设计FIR 滤波器的最直接方法，它通过采用不同时宽的窗函数，对理想滤波器的无限长冲激响应h d (n)进行截短，从而得到系统的有限长冲激响应 h (n)，这一过程可用式5－1来描述：,021-N ||,(n)h )()()(d ⎪⎩⎪⎨⎧≤=其它＝ nn w n h n h R d （5.1）其中W R (n)是时宽为N 的窗函数。

最小， 而既非通带波纹最小， 又非阻带波动最小。 所以， 用这种优化程序设计的滤波器的阻带最小衰减和通带波纹可能 不满足要求。
第５章 无限脉冲响应(IIR)数字滤波器的设
计
特别是以理想滤波器特性作为Hd(ejω)时， 为了使ε2最小，
优化过程尽可能逼近Hd(ejω)的间断特性（即使过渡带最窄）， 而使通带出现较大过冲、 阻带最小衰减过小， 不能满足工
H(ejω)=|H(ejω)|ejθ(ω) 其中， |H(ejω)|称为幅频特性函数， θ(ω)称为相频特性函数。
常用的典型滤波器|H(ejω)|是归一化的， 即|H(ejω)|max=1， 下 的讨论一般就是针对归一化情况的。 对IIR数字滤波器， 通
常用幅频响应函数|H(ejω)|来描述设计指标， 而对线性相位特 性的滤波器， 一般用FIR数字滤波器设计实现。
计
图5.1.6
第５章 无限脉冲响应(IIR)数字滤波器的设
计
5.1.4 IIR-DF的直接设计法
所谓直接设计法， 就是直接在数字域设计IIR[CD*2]DF 的方法。 相对而言， 因为从AF入手设计DF是先设计相应的 AF， 然后再通过s-z平面映射， 将Ha(s)转换成H(z)， 所以 这属于间接设计法。 该设计法只能设计与几种典型AF相对 应的幅频特性的DF。 而需要设计任意形状幅频特性的DF时， 只能用直接设计法。 直接设计法一般都要借助于计算机进行 设计， 即计算机辅助设计（CAD）。 现在已有多种DF优化 设计程序。 优化准则不同， 所设计的滤波器特点亦不同。所 以最主要的是建立优化设计的概念， 了解各种优化准则的 特点， 并根据设计要求， 选择合适的优化程序设计DF。
≤≤
(5.1.1)
≤
(5.1.2)

张卫华 MATLAB课堂
例13：
分析以下程序， 并运行观察。
clf;x=3*pi*(-1:0.05:1);y=x; [X,Y]=meshgrid(x,y); R=sqrt(X.^2+Y.^2)+eps; Z=sin(R)./R; h=surf(X,Y,Z);colormap(jet); axis off n=12;mmm=moviein(n); for i=1:n rotate(h,[0 0 1],25); mmm(:,i)=getframe; end movie(mmm,5,10)
张卫华 MATLAB课堂
例5程序：
fplot('cos(tan(pi*x))',[ 0,1],1e-4)
张卫华 MATLAB课堂
例6：
绘制r=sin(t)cos(t)的极坐标图，并标记 数据点。
张卫华 MATLAB课堂例6 Nhomakorabea序：t=0:pi/50:2*pi; r=sin(t).*cos(t); polar(t,r,'-*');
title('y1=0.2e^{-0.5x}cos(4\pix) 和y2=2e^{0.5x}cos(4\pix)比较') text(2,2*exp(-1),'\fontsize{20}\bf \leftarrow x_{1}=2,y_{1}=2e^{-1}')
张卫华 MATLAB课堂
例5：
绘制f(x)=cos(tan(πx))的曲线
张卫华 MATLAB课堂
例11：
绘制柱形，剪切掉x、y小于零 的部分
张卫华 MATLAB课堂
例11程序：
t=linspace(0,2*pi,100); r=1-exp(-t/2).*cos(4*t); [X,Y,Z]=cylinder(r,60); ii=find(X<0&Y<0); Z(ii)=NaN; surf(X,Y,Z);colormap(spring), shading interp light('position',[-3,-1,3],'style','local')

（1）如果当前时间为 7-8 点则提醒用户，该吃早饭了； （2）如果当前时间为 9-11 点则提醒用户，该学习了； （3）如果当前时间为 12-14 点则提醒用户，该午休了； （4）如果当前时间为 14-17 点则提醒用户，该锻炼了； 答： %%程序为：remind.fig、remind.m

6.求解六元线性方程组： （1）方程组的系数矩阵由用户通过键盘输入； （2）得到系数矩阵后给出方程的解； （3）程序要具有友好性 答：
%% %该模块实现的功能是：求解六元线性方程组，方程组的系数矩阵由用户通过键盘输入；得到系数矩阵后给 出方程的解； clc; clear all; close all; %% %方程输入模块 x=inputdlg({'第一个方程系数','第二个','第三个','第四个','第五个','第六个'}); A=cell2mat(x); %% y=inputdlg({'第一个方程等号右边数','第二个','第三个','第四个','第五个','第六个'});
you=input('请做出你的选择，石头（1） ，剪刀（2） ，布（3） ： '); end disp('您的选择是： '); disp(b(2*you-1:2*you)); compute=ceil(3*rand(1,1)); disp('电脑的选择是：'); disp(b(2*compute-1:2*compute)); %% %%输赢判断模块 %如果电脑与选手出的一样，则显示平手，否则根据石头剪刀布的规则来判断输赢 if you==compute disp('平手'); end switch(you-compute) case{1,-2} disp('您输了'); case{-1,2} disp('您赢了'); end %% %是否继续判别模块 n=input('是否继续玩该游戏？否（0） ，是（1）'); while n~=0&n~=1 disp('您输入的不是正确数字，请正确输入'); n=input('是否继续玩该游戏？否（0） ，是（1）'); end end >> 请选择，石头（1） ，剪刀（2） ，布（3） ： 1 您的选择是： 石头 电脑的选择是： 、剪 您赢了 是否继续玩该游戏？否（0） ，是（1）1 请选择，石头（1） ，剪刀（2） ，布（3） ： 2 您的选择是： 、剪 电脑的选择是： 石头 您输了 是否继续玩该游戏？否（0） ，是（1）0 >> 5.编写一个日程提醒程序实现如下功能：

第五章习题讲解
1
1、用直接I型及典范结构实现以下系统函数：
H
z
=
3 2
4.2 0.6
z z
1 1
0.8 z 2 0.4 z 2
解：根据IIR滤波器的系统函数标准式
M
H
bm zm
z =
m0 N
1 an zn
Y z X z
n 1
将系统函数整理为：
H
z
=
1.5 2.1z1 0.4z2 1 0.3z1 0.2z2
H
z
=
4
z
z 1z 0.5 z2
2 1.4z 1 0.9z 0.8
试问一共能构成几种级联型网络。
解：H
z
A
k
1 1k z1 2k z2 1 1k z1 2k z2
4 1 z1 11.4z1 z2
1 0.5z1 1 0.9z1 0.8z2
4
考虑分子分母的组合及级联的次序，共有以下 四种级联型网络：
=
5
2z3 3z 1 z1
6
抽样点数 N 6，修正半径 r 0.9。
解：由N = 6，得频率抽样型结构：
H z= 1 6
1 r6z6
H0
z
H3
z
2 k 1
H
k
z
又
5 3z3 1 z3
H z=
1 z1
1 z1
1 z1 z2
5 3z3 1 z1 z2
5
h2 1
h1 h3 3 0.6
5
即 hn是偶对称，对称中心在n N 1 2处，
N为奇数 N 5 。
2
得线性相位结构：
17
1.5 2.1z1 0.4z2 1 0.3z1 0.2z2

第一题（1）a=[1 9 8;7 2 5;3 -2 7] %产生矩阵det(a) %检验是否可逆ans=-442,非0，可逆div(a) %求逆矩阵（2）b=[1 0 -7 5;0 -26 7 2;7 4 3 5;8 -3 2 15]det(b)div(b)第二题(1)A=[1 2 3;2 2 5;3 5 1];B=[11 12 31];X=B/A(2) A=[3 1 0 5;0 6 7 3;0 4 3 0;2 -1 2 6];B=[2 4 7 8];X=B/A第三题（1）t=[1 2 3 4 5 6 7 8 9 10]';y=[4.842 4.362 3.754 3.368 3.169 3.083 3.034 3.016 3.012 3.005]';A=[ones(size(t)) exp(-t)];C=A\yT=[0:.1:10]';Y=[ones(size(T)) exp(-T)]*C;plot(T,Y,'-',t,y,'o')title( '采用y(t)≈c1+c2e–t的拟合' )xlabel('\itt'), ylabel('\ity')（2）t=[0 .3 .8 1.1 1.6 2.3]';y=[.82 .72 .63 .60 .55 .50]';A=[ones(size(t)) t.*exp(-t)];C=A\y;T=[0:.1:2.5]';Y=[ones(size(T)) T.*exp(-T)]*C;plot(T,Y,'-',t,y,'o')title('采用y(t)≈d1+d2te–t拟合')xlabel('\itt'), ylabel('\ity')第四题A=[11.59 12.81 15.66; 15.2 4.18 13.61; 10.597.59 9.22];[L,U]=lu(A)[Q R]=qr(A)B=[16.00 4.41 -10.37 -21.61; 0.88 -20.04 12.86 8.56; -1.43 10.71 18.81 -5.99; -12.48 24.35-23.9 10.34];[C,D]=lu(B)[E F]=qr(B)第五题（1）A=[5 -5 -6;3 -2 5;2 -1 -4];x0=[1;-4;5];X=[];for t=0:.01:1X=[X expm(t*A)*x0];endplot3(X(1,:),X(2,:),X(3,:),'-o')grid on(2)A=[1 2 -3 1;3 0 1 -2;1 -2 0 5;2 3 0 1];x0=[1;-1;2;1];X=[];for t=0:.01:1X=[X expm(t*A)*x0];endplot3(X(1,:),X(2,:),X(3,:),'-o')grid on第六题（1）A=[11.59 12.81 15.66; 15.2 4.18 13.61;10.59 7.59 9.22];lambda=eig(A)[V,D]=eig(A)（2）B=[16.00 4.41 -10.37 -21.61; 0.88 -20.04 12.86 8.56; -1.43 10.71 18.81 -5.99; -12.48 24.35 -23.9 10.34];lambda=eig(B)[V,D]=eig(B)第七题（1）x=[1 2 3 4 5 6 7 8 9 10];y=[15.0 39.5 66.0 85.5 89.0 67.5 12.0 -86.4 -236.9 -448.4];p=polyfit(x,y,2);x2=1:.1:10;y2=polyval(p,x2);figure(1)plot(x,y,'o',x2,y2)grid ontitle('二阶多项式曲线拟合')(2)x=[1 2 3 4 5 6 7 8 9 10];y=[15.0 39.5 66.0 85.5 89.0 67.5 12.0 -86.4 -236.9 -448.4];p=polyfit(x,y,3);x2=1:.1:10;y2=polyval(p,x2);figure(1)plot(x,y,'o',x2,y2)grid ontitle('三阶多项式曲线拟合')第八题p1=[1,-2-3,4,2];p2=[1,-7,5,31,-30];p3=[1,-1,-25,25];p4=[-2,3,1,5,8,0];[L1,U1]=lu(p1)r1=roots(p1)[L2,U2]=lu(p2)r2=roots(p2)[L3,U3]=lu(p3)r3=roots(p3)[L4,U4]=lu(p4)r4=roots(p4)第九题p1=[1,-2-3,4,2];p2=[1,-7,5,31,-30];p3=[1,-1,-25,25];p4=[-2,3,1,5,8];p1_x=polyval(p1,[-1.5,2.1,3.5]) p2_x=polyval(p2,[-1.5,2.1,3.5]) p3_x=polyval(p3,[-1.5,2.1,3.5]) p4_x=polyval(p4,[-1.5,2.1,3.5])第十题a=[2,3,-4];b=[4,-2,5];c=[3,0,-2,5,6];d1=conv(a,b)[d2,r2]=deconv(c,a)[d3,r3]=deconv(c,b)第十一题a=[2,3,-4];b=[4,-2,5];c=[3,0,-2,5,6];dao1=polyder(a,b)[dao2,r2]=polyder(c,a)[dao3,r3]=polyder(c,b)第十二题x=-5:.25:5;y=10*exp(-x);xi=-5:5;y1=interp1(x,y,xi,'nearest');y2=interp1(x,y,xi,'linear');y3=interp1(x,y,xi,'spline');y4=interp1(x,y,xi,'cubic'); figure(1);subplot(2,2,1)plot(x,y,'-',xi,y1,'o');title('最邻近内插');grid on;xlabel('x');ylabel('y');subplot(2,2,2)plot(x,y,'-',xi,y2,'o');title('线性内插');grid on;xlabel('x');ylabel('y');subplot(2,2,3)plot(x,y,'-',xi,y3,'o');title('三次样条内插');grid on;xlabel('x');ylabel('y');subplot(2,2,4)plot(x,y,'-',xi,y4,'o');title('三次曲线内插');grid on;xlabel('x');ylabel('y');第十三题x=rand(1,50);y=randn(1,50);minx=min(x)miny=min(y)maxx=max(x)maxy=max(y)avx=mean(x)avy=mean(y)Ex=(std(x)).^2Ey=(std(y)).^2第十四题t=[0 .2 .4 .6 .8 1.0 2.0 5.0 ]';y=[1.0 1.51 1.88 2.13 2.29 2.40 2.60 24.00]'; X1=[ones(size(t)) t t.^2];a=X1\y;X2=[ones(size(t)) exp(-t) t.*exp(-t)];b=X2\y;T=[0:.1:6]';Y1=[ones(size(T)) T T.^2]*a;Y2=[ones(size(T)) exp(-T) T.*exp(-T)]*b; figure(1)subplot(1,2,1)plot(T,Y1,'-',t,y,'o'),grid ontitle('多项式回归')subplot(1,2,2)plot(T,Y2,'-',t,y,'o'),grid ontitle('指数函数回归')第十五题t=0:1/119:1;x=3*sin(2*pi*20*t)+10*sin(2*pi*200*t+pi/4)+10*randn(size(t)); y=fft(x);m=abs(y);f=(0:length(y) -1)'*119/length(y);figure(1)subplot(2,1,1),plot(t,x),grid ontitle('被噪声污染的信号')ylabel('Input \itx'),xlabel('Time ')subplot(2,1,2),plot(f,m)ylabel('Abs. Magnitude'),grid onxlabel('Frequency (Hertz)')第十六题w=input('w=');t=0:1/119:1;x1=sin(w.*t)+randn(size(t));x2=cos(w.*t)+randn(size(t));x3=sin(w.*t)+randn(size(t));a=corrcoef(x1,x2)b=corrcoef(x1,x3)若没有正弦分量w=input('w=');t=0:1/119:1;x1=randn(size(t));x2=randn(size(t));x3=randn(size(t));a=corrcoef(x1,x2)b=corrcoef(x1,x3)第十七题z1=quad('exp(-2*t)',0,2)z2=quad('exp(2*t)',0,2)z3=quad('exp(t.^2-3*t+.5)',-1,1)第十八题function y=five(x)y=exp(-x)-1.5*exp(2*cos(2*x));%主函数x0=input('x0='); %执行时，按要求输入[-1,1]z=fzero('five',x0)第十九题function f=five(x,y)f=exp(-x.*y)-2*x.*y;%主函数z=dblquad('five',0,1,-1,1)第二十题function dy=five(t,y)dy=[0.5-y(1);y(1)-4*y(2)];%主函数X0=[1; -0.5];tspan=[0,25];[T,X]=ode45('five',tspan,X0);figure(1)subplot(2,1,1),plot(T,X(:,1),'r'),title('x_{1}'),grid onsubplot(2,1,2),plot(T,X(:,2),'k'),title('x_{2}'),grid onfigure(2)plot(X(:,1),X(:,2)),title('系统轨迹'),grid onxlabel('x_{1}'),ylabel('x_{2}')。

数字信号处理第五章习题解答————第五章————数字滤波网络5.1 学习要点本章主要介绍数字滤波器的系统函数()z H 与其网络结构流图之间的相互转换方法，二者之间的转换关系用Masson 公式描述。

由于信号流图的基本概念及Masson 公式已在信号与系统分析课程中讲过，所以下面归纳IIR 系统和FIR 系统的各种网络结构及其特点。

5.1.1 IIR 系统的基本网络结构1. 直接型结构如果将系统函数()z H 化为标准形式(5.1)式：()∑∑=-=--=Nk kkMk kkz az bz H 11 (5.1) 则可根据Masson 公式直接画出()z H 的直接II 型网络结构流图如图5.1所示（取N=4，M=3）。

二阶直接II 型网络结构最有用，它是级联型和并联型网络结构的基本网络单元。

优点：可直接由标准形式(5.1)或差分方程()()()∑∑==-+-=Mk kN k kk n x b k n y a n y 01画出网络结构流图，简单直观。

缺点：对于高阶系统：（1）调整零、极点困难；（2）对系数量化效应敏感度高；（3）乘法运算量化误差在系统输出端的噪声功率最大。

2. 级联型结构将(5.1)式描述的系统函数()z H 分解成多个二阶子系统函数的乘积形式()()()()z H z H z H z H m 21?= (5.2) (),1221122110------++=zzzzz H i i i i i i ααβββ m i ,,2,1 = (5.3)画出的级联型方框图如图5.2所示。

图中每一个子系统均为二阶直接型结构，根据()z H 的具体表达式确定()z H i 的系数i i i i 1210,,,αβββ和i 2α后，可画出()z H i 的网络结构流图如图5.3所示。

优点：（1）系统结构组成灵活；（2）调整零、极点容易，因为每一级二阶子系统()z H i 独立地确定一对共轭零点和一对共轭极点；（3）对系数量化效应敏感度低。

2011/4/10
3
Notice: Skipped Sections

5.7 5.8 (all) 5.9 (all) 5.10 5.11
2011/4/10
4
5.1
Introduction
2011/4/10
5
5.1 Introduction
Digital processing of a real-world continuous-time signal involves the following basic steps:
The multiplication operation yields an impulse train
Note: Analog Filters are very important for digital systems. Since both the anti-aliasing filter and the reconstruction filter are analog lowpass filters, we review first the theory behind the design of such filters Also, the most widely used IIR digital filter design method is based on the conversion of an analog lowpass prototype
g[n ] = ga (nT ), -¥ <n < ¥
(5.1)
With T being the sampling period. The reciprocal of T is called the sampling frequency , i.e., 1⁄ . It is known that the frequency-domain of the analog signal is given by its FT:

第5章MＡTLAB绘图习题5一、选择题1.如果x、ｙ均为４×3矩阵，则执行plot（ｘ,ｙ）命令后在图形窗口中绘制( ）条曲线。

DA。

12 B。

７ C.4 Ｄ.32.下列程序得运行结果就是（）。

Aｘ=0:pi/１00:２＊ｐi;for n=1:2：１０ｐlot（n＊ｓiｎ（x），n＊cos（x)）hold oneｎｄaｘｉｓsｑuarｅA。

5个同心圆Ｂ.5根平行线Ｃ．一根正弦曲线与一根余弦曲线 D.5根正弦曲线与5根余弦曲线3．命令ｔｅxt（1,1，'{\ａｌpｈa｝+｛\ｂeta｝’）执行后,得到得标注效果就是( )。

CＡ.{\ａlｐha}+｛\bｅta｝B。

{\α}＋{\β｝ C.α+βＤ.＼α＋\β4。

subplot(２，2，3）就是指( )得子图。

AA.两行两列得左下图Ｂ.两行两列得右下图Ｃ.两行两列得左上图Ｄ．两行两列得右上图5。

要使函数ｙ=2ｅx得曲线绘制成直线，应采用得绘图函数就是(）.ＣA.pｏｌａｒB。

ｓeｍiloｇｘC。

semiｌogy D。

loｇlｏg６.下列程序得运行结果就是( )。

B[x，ｙ］＝meshｇｒｉd(1：５);surf（x,ｙ,5＊ｏnes(sｉze（x)）)；A.z=x+y平面B.与xy平面平行得平面Ｃ。

与xy平面垂直得平面 D.z=5x平面７.下列函数中不能用于隐函数绘图得就是（)。

ＤA.ezｍeｓhB.eｚsurｆC.ｅzｐloｔD．pｌoｔ38.下列程序运行后，瞧到得图形（）.Ct=0：pi/20：2*pi；［ｘ，ｙ］＝ｍeｓhgrid（-8：0、5:8）;ｚ=siｎ（sqrt（ｘ、＾2+y、＾2）)、／sｑrt（x、^2＋ｙ、^2+eps)；surf（ｘ，y,z)viｅw（0,9０)；axiｓ eｑualA．像墨西哥帽子 B.就是空心得圆C。

边界就是正方形D．就是实心得圆９。

下列程序运行后得到得图形就是( ）.Ａ[x，y］=meshｇrid（－2：2）；z=ｘ＋y;i=ｆind（aｂs（ｘ）<1 ＆ abs（y）<1);z（i）=ＮaＮ;ｓurf（x,ｙ，z)；ｓhａdinｇｉｎterpA。