哈尔滨工程大学数字信号处理实验二

-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50510MAGNITUDERESPONSENORMALIZEDFREQUENCY|H(w)|-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5-100-50050100NORMALIZEDFREQUENCYDEGREESPHASERESPONSE实验原理对离散时间傅里叶变换DTFT来说有两个问题：（1）DTFT的定义对无限长信号时有效的；（2）DTFT是连续变量w的函数。在MATLAB中，任何信号必须是有长度的，使（1）成为问题，当能从变换定义式推导出解析式并只是计算它，可以使用MATLAB计算无限长信号的DTFT。对于频率抽样问题，通常在有限网格点上选择足够多的频率一试绘出逼近真实的DTFT。最好的选择是在（-pi，pi）区间上一组均匀隔开的频率，如式3.11所示。在数量有限的频率点wk=2*pi*k/N出计算，并在有限范围内求和。要求N=L。MATLAB的freqz函数足以应对无限长信号，但计算有限长信号的DTFT需用一个新的函数，dtft函数。对于无限长信号中的指数信号，其DTFT是exp(j*w)的有理函数。可以推导出|a|1时h[n]=a^n*u[n]的DTFT表达式H(exp(j*w))=1/(1-a*exp(-j*w))，使用这一有理式容易计算这个DTFT的一组样本，分子和分母的计算都可以使用FFT，即计算有理式等同于做两个dtft函数计算，定义freqz函数[HH,WW]=freqz(b,a,N,'whole')。建立dtft函数如下所示：function[H,W]=dtft(h,N)%UNTITLED2Summaryofthisfunctiongoeshere%DetailedexplanationgoeshereN=fix(N);L=length(h);h=h(:);if(NL)error('DTFT:#datasamplescannotexceed#freqsamples')endW=(2*pi/N)*[0:(N-1)]';mid=ceil(N/2)+1;W(mid:N)=W(mid:N)-2*pi;W=fftshift(W);H=fftshift(fft(h,N));在变换域上使用subplot绘制双子图，程序如下：formatcompact,subplot(111)a=0.88*exp(sqrt(-1)*2*pi/5);nn=0:40;xn=a.^nn;[X,W]=dtft(xn,128);subplot(211),plot(W/2/pi,abs(X));grid,title('MAGNITUDERESPONSE')xlabel('NORMALIZEDFREQUENCY'),ylabel('|H(w)|')subplot(212),plot(W/2/pi,180/pi*angle(X));gridxlabel('NORMALIZEDFREQUENCY'),ylabel('DEGREES')title('PHASERESPONSE')一、脉冲信号的DTFTa．证明3.12的矩形脉冲DTFT可由3.13数学表达式得出，改变换的第一项成为混叠函数。程序如下：formatcompact,subplot(111)%利用dtft函数求矩形序列的傅里叶变换，取N=200a=1;%并利用subplot将傅里叶变换的幅频特性和相频特性nn=0:19;xn=a.^nn;%在双子图上显示出来。[X,W]=dtft(xn,200);subplot(211),plot(W/2/pi,abs(X));grid,title('MAGNITUDERESPONSE')xlabel('NORMALIZEDFREQUENCY'),ylabel('|H(w)|')subplot(212),plot(W/2/pi,180/pi*angle(X));gridtitle('PHASERESPONSE')xlabel('NORMALIZEDFREQUENCY'),ylabel('DEGREES')图形显示如下：w=W;R=sin(w*20/2)./sin(w/2).*exp(-j.*w*19/2);%验证矩形序列的DTFT可由式3.13得出subplot(211),plot(w,abs(R));%并利用subplot将傅里叶变换的幅频特性xlabel('w'),ylabel('abs')%和相频特性在双子图上显示出来。subplot(212),plot(w,180/pi*angle(R));xlabel('w'),ylabel('phase')-0.5-0.4-0.3-0.2-0.100.10.20.30.40.505101520MAGNITUDERESPONSENORMALIZEDFREQUENCY|H(w)|-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5-200-1000100200PHASERESPONSENORMALIZEDFREQUENCYDEGREES-4-3-2-101234-2024681012wreal-4-3-2-101234-10-8-6-4-20246810wimagb．使用dtft函数计算12点脉冲信号的dtft。绘出在区间pi=w=pi上对w的DTFT。把实部和虚部分开绘出。另绘出DTFT的幅度。选择样本的数量是脉冲长度的5到10倍，以使绘出的图看上去平滑。程序及图形显示如下：1．绘出实部formatcompact,subplot(111)a=1;nn=0:11;xn=a.^nn;[X,W]=dtft(xn,120);plot(W,real(X));xlabel('w'),ylabel('real')2.绘出虚部formatcompact,subplot(111)a=1;nn=0:11;xn=a.^nn;[X,W]=dtft(xn,120);-4-3-2-10123405101520wabs-4-3-2-101234-200-1000100200wphase-4-3-2-101234024681012wabsplot(W,imag(X));xlabel('w'),ylabel('imag')3.绘制DTFT的幅度formatcompact,subplot(111)a=1;nn=0:11;xn=a.^nn;[X,W]=dtft(xn,120);plot(W,abs(X));xlabel('w'),ylabel('abs')C．注意asinc函数零点的位置是规则分布的，对奇数脉冲，L=15的脉冲重复进行DTFT计算并绘出幅度；同样再次检验零点位置，注意峰值高度。程序及图形显示如下：-4-3-2-101234051015wabs-4-3-2-101234051015wabs1.利用dtft函数formatcompact,subplot(111)a=1;nn=0:14;xn=a.^nn;[X,W]=dtft(xn,200);plot(W,abs(X));xlabel('w'),ylabel('abs')2.利用式3.13w=W;R=sin(w*15/2)./sin(w/2).*exp(-j.*w*14/2);plot(w,abs(R));xlabel('w'),ylabel('abs')-4-3-2-101234051015wabs结论：在w=0处，峰值高度为脉冲长度。d．对asinc函数零点的间距与asinc函数的直流值确定出通用规则。二、asinc的M文件实验内容：编写一个MATLAB的函数如asinc（w，L），直接从式3.14计算在频率格上的asinc（w，L）。应该有w频率和L长度两个输入量。函数必须检查w=0de情形。直接计算混叠sinc函数式3.13得到的脉冲信号DTFT，绘出幅度。建立asinc函数如下所示：function[x]=asinc(w,L)%UNTITLED2Summaryofthisfunctiongoeshere%Detailedexplanationgoeshereifw==0x=L;%当w=0时，幅值为脉冲长度elsex=sin(w*L/2)./sin(w/2);end利用asinc函数得到矩形脉冲的傅里叶变换幅值的程序：w=W;y=asinc(w,15);x=y.*exp(-j.*w*14/2);plot(w,abs(x));xlabel('w'),ylabel('abs')%利用asinc函数得到矩形脉冲%的傅里叶变换幅值四、指数信号对于信号x[n]=(0.9)^n*u[n],使用freqz函数计算其DTFT。a.对w在区间-pi=w=pi上绘出幅度与相位特性，这需要从freqz返回的[X,W]向量的移位。a=[1,-0.9];%H(exp(j*w))=1/(1-a*exp(-j*w))分母系数对应的向量b=1;%上式分子N=100;%0到2pi共有100个抽样点[HH,WW]=freqz(b,a,N,'whole');%利用freqz函数求无限长指数信号的傅里叶变换mid=ceil(N/2)+1;%ceil为从0朝正无穷大方向取整WW(mid:N)=WW(mid:N)-2*pi;%此句为移位操作WW=fftshift(WW);%频率格点WW和变换数值HH调整，让正半轴部分和HH=fftshift(HH);%负半轴部分的图像分别关于各自的中心对称。subplot(211),plot(WW/2/pi,abs(HH));%将幅频特性和相频特性在双子图上显示subplot(212),plot(WW/2/pi,180/pi*angle(HH));-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5051015-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5-100-50050100051015202530-1-0.500.51nreal051015202530-1-0.500.51nimag-0.5-0.4-0.3-0.2-0.100.10.20.30.40.502468101214161820五、复指数信号a.取z0=0.95*exp(j*3*pi/11),0=n=30,绘出x[n]=z0^n*u[n]。用subplot指令将实部和虚部对n绘制双子图。程序如下：nn=0:30;z0=0.95*exp(j*3*pi/11);x=z0.^nn;subplot(211),plot(nn,real(x));xlabel('n'),ylabel('real')subplot(212),plot(nn,imag(x));xlabel('n'),ylabel('imag')%用subplot指令将实部和%虚部在双子图上绘制b.再一次取z0=0.95*exp(j*3*pi/11),计算其DTFT，并对w画出幅度。注意幅度响应的尖峰位于何处，吧尖峰位置与z0的极坐标形式联系起来。程序如下：formatcompact,subplot(111)a=[1,-0.95*exp(j*3*pi/11)];b=1;N=200;[HH,WW]=freqz(b,a,N,'whole');mid=ceil(N/2)+1;WW(mid:N)=WW(mid:N)-2*pi;WW=fftshift(WW);HH=fftshift(HH);plot(WW/2/pi,abs(HH));%利用subplot指令绘制%z0=0.95*exp(j*3*pi/11)的幅频特性-0.5-0.4-0.3-0.2-0.10