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

数字信号处理MATLAB

第1页

习题数字信号处理MATLAB习题

M1-1 已知g1(t)cos(6t)，g2(t)cos(14t)，g3(t)cos(26t)，以抽样频率fsam10Hz对上述三个信号进行抽样。在同一张图上画出g1(t)，

g2(t)和g3(t)及抽样点，对所得结果进行讨论。

解：

第2页

从以上两幅图中均可看出，三个余弦函数的周期虽然不同，但它们抽样后相应抽样点所对应的值都相同。那么这样还原回原先的函数就变成相同的，实际上是不一样的。这是抽样频率太小的原因，我们应该增大抽样频率才能真实还原。如下图：f=50Hz

第3页

程序代码

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(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.',k,h3,'b.')

第4页

xlabel('t'); ylabel('g(t)');

legend('g1(t)','g2(t)','g3(t)');

M2-1 利用DFT的性质，编写一MATLAB程序，计算下列序列的循环卷积。

(1) g[k]={1,-3,4,2,0,-2,},h[k]={3,0,1,-1,2,1}； (2) x[k]=cos(k/2),y[k]=3k,k=0,1,2,3,4,5。 解： (1)

循环卷积结果

6.0000 -3.0000 17.0000 -2.0000 7.0000 -13.0000

程序代码

第5页

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,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); end end y2

stem(0:l-1,y2) xlabel('k') ylabel('y(k)') title('循环卷积')

(2)

循环卷积结果

-71.0000 -213.0000 .0000 267.0000 73.0000 219.0000

第6页

程序代码

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<=L

y2(n)=y1(n)+y1(n+l); else

y2(n)=y1(n); end end y2

stem(0:l-1,y2) xlabel('k') ylabel('y’(k)') title('循环卷积')

第7页

M2-2 已知序列x[k]cos(k/2N),|k|N

0,其他(1)计算序列DTFT的表达式X(ej),并画出N=10时，X(ej)的曲线。

(2)编写一MATLAB程序，利用fft函数，计算N=10时，序列x[k]的DTFT在m2m/N的抽样值。利用hold函数，将抽样点画在X(ej)的曲线上。

解：

(1) X(e)DTFT{x[k]}jkx[k]ejkkNcos(k/2N)eNjk

程序代码

N=10; k=-N:N;

x=cos(k.*pi./(2*N)); W=linspace(-pi,pi,512);

第8页

X=zeros(1,length(W)); for k=-N:N

X1=x(k+N+1).*exp(-j.*W.*k); X=X+X1; end

plot(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; end

第9页

plot(W,abs(X));hold on;

plot(2*pi*L/21,fftshift(abs(X_21)),'o'); xlabel('W'); ylabel('abs(X)');

M2-3 已知一离散序列为x[k]Acos0kBcos[(0)k]。用长度N=的Hamming窗对信号截短后近似计算其频谱。试用不同的A和B的取值，确定用Hamming窗能分辨的最小的谱峰间隔wc的值。

解：f1=100Hz f2=120Hz时

2中cN

f2=140Hz时

第10页

f2=160Hz时

第11页

由以上三幅图可见

f2=140Hz时，各谱峰可分辨。则f又

wc2N

40Hz

且

wT2fT2401 800所以c=3.2（近似值）

程序代码

N=; 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*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);

第12页

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 已知一离散序列为x[k]cos0k0.75cos1k,0k63。其中, 02/15，12.3/15。

(1) 对x[k]做点FFT, 画出此时信号的谱。

(2) 如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的点信号补0而分辨出两个谱峰。通过编程进行证实，并解释其原因。

解： (1)

第13页

程序代码

W0=2*pi/15; W1=2.3*pi/15; N=; k=0:N-1;

x=cos(W0*k)+0.75*cos(W1*k); X=fft(x);

plot(k/N,abs(X)); grid on;

title('点FFT');

(2)

第14页

第15页

由以上三幅图看出：不能对(1)中的点信号补零而分辨出两个谱峰，这样的方法只能改变屏幕分辨率，但可以通过加hamming窗来实现对谱峰的分辨。 程序代码

W0=2*pi/15; W1=2.3*pi/15; N=; 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近似分析

第16页

其频谱。若要求频率分辨率为1Hz，试确定抽样频率fsam、抽样点数N以及持续时间Tp。

解：

本题使用矩形窗，则Nfsamfsam1fsam，Tp1 f1f

第17页

由以上三幅图可以看出当fsam越来越大时，近似值越来越接近

第18页

于实际值。即fsam越大拟合效果越好，造成的混叠也是在可以允许的范围内。 程序代码

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('近似值','实际值);

M2-6 试用DFT近似计算高斯信号g(t)exp(dt2)的频谱抽样值。

π2通过和频谱的理论值G(j)exp()比较，讨论如何根据时域的信

d4d号来恰当地选取截短长度和抽样频率使计算误差能满足精度要求。

解：

第19页

第20页

由以上三幅图可以看出：

当时域截取长度相同时，抽样间隔越小时误差越小，当抽样间隔一定时，时域截取长度越长，误差越小。当取抽样间隔为1S，时域截取长度为2S时，误差较大，绝对误差在0.5左右；当抽样间隔为0,5S，时域截取长度为2S时，误差比间隔为1S时小，绝对误差不大于0.2；当抽样间隔为0.5S时域截取长度为4S时，误差更小，绝对误差不大于0.04。因为时域截取长度越长，保留下来的原信号中的信息越多，抽样间隔越小，频谱越不容易发生混叠，所以所得频谱与理论值相比，误差更小。

程序代码

Ts=0.5; N=4; N0=;

k=(-N/2:(N/2))*Ts;

第21页

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) subplot(2,1,2) plot(w/pi,abs(X)-XT) ylabel('实验误差')

xlabel('\\omega/\\pi');

'N=',num2str(N)]);第22页

' '