数字信号处理课后习题Matlab作业.docx

1、数字信号处理课后习题Matlab作业数字信号处理MATLAB习题数字信号处理MATLAB习题M1-1 已知，以抽样频率对上述三个信号进行抽样。在同一张图上画出，和及抽样点，对所得结果进行讨论。解：从以上两幅图中均可看出，三个余弦函数的周期虽然不同，但它们抽样后相应抽样点所对应的值都相同。那么这样还原回原先的函数就变成相同的，实际上是不一样的。这是抽样频率太小的原因，我们应该增大抽样频率才能真实还原。如下图：f=50Hz程序代码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

2、:0.2;h1=cos(6.*pi.*k);h2=cos(14.*pi.*k);h3=cos(26.*pi.*k);% subplot(3,1,1);% plot(k,h1,r.,t,g1,r);% xlabel(t);% ylabel(g1(t);% subplot(3,1,2);% plot(k,h2,g.,t,g2,g);% xlabel(t);% ylabel(g2(t);% subplot(3,1,3);% plot(k,h3,b.,t,g3,b);% xlabel(t);% ylabel(g3(t);plot(t,g1,r,t,g2,g,t,g3,b,k,h1,r.,k,h2,g.

3、,k,h3,b.)xlabel(t);ylabel(g(t);legend(g1(t),g2(t),g3(t);M2-1 利用DFT的性质，编写一MATLAB程序，计算下列序列的循环卷积。(1) gk=1,-3,4,2,0,-2,hk=3,0,1,-1,2,1；(2) xk=cos(k/2),yk=3k,k=0,1,2,3,4,5。解：(1)循环卷积结果6.0000 -3.0000 17.0000 -2.0000 7.0000 -13.0000程序代码g=1 -3 4 2 0 -2;h=3 0 1 -1 2 1;l=length(g);L=2*l-1;GE=fft(g,L);HE=fft(h,

4、L);y1=ifft(GE.*HE);for n=1:l if n+l=L y2(n)=y1(n)+y1(n+l); else y2(n)=y1(n); endendy2stem(0:l-1,y2)xlabel(k)ylabel(y(k)title(循环卷积)(2)循环卷积结果-71.0000 -213.0000 89.0000 267.0000 73.0000 219.0000程序代码k=0:5;x=cos(pi.*k./2);y=3.k;l=length(x);L=2*l-1;GE=fft(x,L);HE=fft(y,L);y1=ifft(GE.*HE);for n=1:l if n+l=

5、L y2(n)=y1(n)+y1(n+l); else y2(n)=y1(n); endendy2stem(0:l-1,y2)xlabel(k)ylabel(y(k)title(循环卷积)M2-2 已知序列(1)计算序列DTFT的表达式,并画出N=10时，的曲线。(2)编写一MATLAB程序，利用fft函数，计算N=10时，序列xk的DTFT在的抽样值。利用hold函数，将抽样点画在的曲线上。解：(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:N X1=x(k+

6、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:N X1=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(a

7、bs(X); M2-3 已知一离散序列为。用长度N=64的Hamming窗对信号截短后近似计算其频谱。试用不同的A和B的取值，确定用Hamming窗能分辨的最小的谱峰间隔中c的值。解：f1=100Hzf2=120Hz时f2=140Hz时f2=160Hz时由以上三幅图可见f2=140Hz时，各谱峰可分辨。则又且所以c=3.2（近似值）程序代码N=64;L=1024;f1=100;f2=160;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

8、*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

9、*pi);subplot(2,2,1);plot(W,abs(X1);title(A=1,B=1);xlabel(W);ylabel(X1);subplot(2,2,2);plot(W,abs(X2);title(A=1,B=0.5);xlabel(W);ylabel(X2);subplot(2,2,3);plot(W,abs(X3);title(A=1,B=0.25);xlabel(W);ylabel(X3);subplot(2,2,4);plot(W,abs(X4);title(A=1,B=0.05);xlabel(W);ylabel(X4);M2-4 已知一离散序列为,0k63。其中,

10、，。(1) 对xk做64点FFT, 画出此时信号的谱。(2) 如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的64点信号补0而分辨出两个谱峰。通过编程进行证实，并解释其原因。解：(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)由以上三幅图看出：不能对(1)中的64点信号补零而分辨出两个谱峰，这样的方法只能改变屏幕分辨率，但可以通过加hamming窗来实现对谱峰的分辨。程序代码W0=2*pi/1

11、5;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 已知一连续信号为x(t)=exp(-3t)u(t)，试利用DFT近似分析其频谱。若要求频率分辨率为1Hz，试确定抽样频率fsam、抽样点数N以及持续时间Tp。解：本题使用矩形窗，则， 由以上三幅图可以看出当fsam越来越大时，近似值越来越接近于实际值。即fsam越大拟合效果越好，造成的混叠也是在可以允许的范围内。程序代码fs=100;ws=2*pi*f

12、s;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(近似值,实际值);M2-6 试用DFT近似计算高斯信号的频谱抽样值。通过和频谱的理论值比较，讨论如何根据时域的信号来恰当地选取截短长度和抽样频率使计算误差能满足精度要求。解：由以上三幅图可

13、以看出：当时域截取长度相同时，抽样间隔越小时误差越小，当抽样间隔一定时，时域截取长度越长，误差越小。当取抽样间隔为1S，时域截取长度为2S时，误差较大，绝对误差在0.5左右；当抽样间隔为0,5S，时域截取长度为2S时，误差比间隔为1S时小，绝对误差不大于0.2；当抽样间隔为0.5S时域截取长度为4S时，误差更小，绝对误差不大于0.04。因为时域截取长度越长，保留下来的原信号中的信息越多，抽样间隔越小，频谱越不容易发生混叠，所以所得频谱与理论值相比，误差更小。程序代码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(jomega);legend(试验值,理论值);title(Ts=,num2str(Ts) N=,num2str(N);subplot(2,1,2)plot(w/pi,abs(X)-XT)ylabel(实验误差)xlabel(omega/pi);