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

1、 实 验 报 告课程名称实验项目名称实验类型实验学时班级学号姓名指导教师实验室名称实验时间实验成绩预习部分实验过程表现实验报告部分总成绩教师签字日期实验二 离散时间傅立叶变换 实验报告1、 实验原理 1、经由正、逆离散时间傅里叶变换表达的信号傅里叶表达式是信号分析的一个关键部分，下面方程分别是分析方程和综合方程。X()=类似地，当LTI系统用于滤波时，作为冲击响应离散时间傅里叶变换的频率响应，提供了LTI系统简介的描述。离散时间傅里叶变换X()是的周期复值函数，周期总是2，并且基周期通常选在区间-，）上。对离散时间傅里叶变换DTFT来说有两个问题：第一个问题是：DTFT的定义对无限长信号是有效

2、的第二个问题是：DTFT是连续变量的函数在MATLAB中，任何信号必须是有限长度的，这就使第一点成为问题。因此，不可能使用MATLAB计算无限长信号的DTFT。有一个值得注意的例外情形，当能从变换定义式推导出解析式并只是计算它时，可以使用MATLAB计算无限长信号的DTFT。例如在xn=un，xn具有有理的DTFT的情形下。 2、第二个问题是频率抽样问题。MATLAB擅长在有限网格点上计算DTFT。通常选择足够多的频率以使绘出的图平滑，逼近真实的DTFT。对计算有利的最好选择是在（-，）区间上一组均匀地隔开的频率，或者对共轭对称变换选择0,区间。采用上述抽样办法，DTFT式变成X()= （3.

3、11）DTFT的周期性意味着在-<0区间上的数值是那些对k>N/2的数值。因为上式是在有限数量的频率点=2k/N处计算，并在有限范围内求和，因此它是可计算的。由于信号长度必须是有限的（0n<L），这个求和式不适用于xn=un的情形。在对DTFT进行抽样时，并不要求N=L，尽管通常经由DFT进行计算（使用FFT算法）时，如果N=L很便利。的确，在使用N点DFT时，如果N>L，只须想像xn是增补零的即可；如果N<L，会较复杂。这种情形下，在正确应用FFT计算N点DFT前，需要对xn进行时间混叠，因此此时应确保总是比远处是时间序列（即N>=L）点数多得多的频率处来

4、计算DTFT。3、 计算DTFT需要两个函数，MATLAB的freqz函数计算无限长信号，dtft（h，H）函数计算有限长信号的DTFT。二：实验内容 实验内容包括脉冲信号的DTFT，asinc的M文件，无限长信号的DTFT，指数信号，复指数信号这五个部分。1. 脉冲信号的DTFTa、已知矩形脉冲rn的定义式为在0n<L上rn=1，n为其余值时，rn均为0。将rn代入DTFT的分析方程X()=中，计算化简后可得：。b、建立DTFT函数（.M文件）实验程序清单1：functionA,B=dtft(h,N)%DTFT calculate DTFT at N equally spaced fr

5、equencies%Usage:%A,B=dtft(h,N)%h:finite-length input vector,whose length is L%N:number of frequencies for evaluation over -pi,pi)%=>constraint:N>=L%H:DTFT values(complex)%W:(2nd output)vector of freqs where DTFT is computed%N=fix(N);L=length(h); h=h(:); %<-for vectors ONLY!if(N<L) error(

6、'DTFT: # data samples cannot exceed # freq samples')endB=(2*pi/N)*0:(N-1)'mid=ceil(N/2)+1;B(mid:N)=B(mid:N)-2*pi;%movepi,2pi)to-pi,0)B=fftshift(B);A=fftshift(fft(h,N);%move negative freq components实验程序清单2:（当w=80时）h=ones(1,12);A,B=dtft(h,80); %w=80时subplot(311),plot(B,real(A); %把绘图窗口分成三行一

7、列三个子图，然后在第一块区域作图xlabel('B'),ylabel('A,real part'); %以B为横轴，A的实部为纵轴subplot(312),plot(B,imag(A); %把绘图窗口分成三行一列三个子图，然后在第二块区域作图xlabel('B'),ylabel('A,imag part'); %以B为横轴，A的虚部为纵轴subplot(313),plot(B,abs(A); %把绘图窗口分成三行一列三个子图，然后在第三块区域作图xlabel('B'),ylabel('A amplitude

8、'); %以B为横轴，A的幅值为纵轴实验仿真图：实验程序清单3:（当w=100时）h=ones(1,12);A,B=dtft(h,100); %w=100时subplot(311),plot(B,real(A);xlabel('B'),ylabel('A,real part');subplot(312),plot(B,imag(A);xlabel('B'),ylabel('A,imag part');subplot(313),plot(B,abs(A);xlabel('B'),ylabel('A a

9、mplitude');实验仿真图：实验程序清单4:（当w=120时）h=ones(1,12);A,B=dtft(h,120); %w=120时subplot(311),plot(B,real(A);xlabel('B'),ylabel('A,real part');subplot(312),plot(B,imag(A);xlabel('B'),ylabel('A,imag part');subplot(313),plot(B,abs(A);xlabel('B'),ylabel('A amplitud

10、e');实验仿真图：实验结果分析：由以上三组不同频率（w=80,w=100,w=120）的仿真图形可知，当w越大，取的点数越多时，绘制的图形越来越平滑。c、注意asinc函数零点的位置是规则分布的。对奇数长脉冲，比如L=15的脉冲重复进行DTFT计算并绘出幅度；同样再次检验零点位置，注意峰值高度。实验程序清单：n=0:14;h=ones(1,15);X,W=dtft(h,150);subplot(211),plot(W,abs(X); %把绘图窗口分成二行一列二个子图，然后在第一块区域作图grid; %图层为网格xlabel('W'),ylabel('|X(w)

11、|'); %X轴为w，Y轴为X(w)幅值H,X=dtft(X,150);subplot(212),plot(X,abs(H); %把绘图窗口分成二行一列二个子图，然后在第二块区域作图grid;xlabel('X'),ylabel('|H(w)|');%X轴为X，Y轴为H(w)幅值实验仿真图： 实验结果分析：由仿真结果可以看出函数零点关于x=0对称。在对称轴处幅值最大，离对称轴越远，幅值逐渐减小。d、对asinc函数零点的间距与asinc函数的直流值，确定出通用规则。 零点间距为/6，由数学推导：在a部分中，当L=12时，由R()=0可得sin(wL/2)

12、=0，即wL/2=k* 则w=k*/6。故得出以上结论。在实验a部分的仿真图中，可观察到幅值高度为12。则零点间距*幅值=(/6)*12=2*。2、 asinc的M文件编写一个MATLAB的函数如asinc(w,L),直接从混叠函数表达式中计算在频率格上的asinc(w,L)。该函数应该有两个输入：唱的L和频率w的向量。函数必须检查被零除的情况，如w=0时。直接计算混叠sinc函数得到脉冲信号的DTFT。绘出幅度，保存该图以便将其与用dtft得到的结果进行比较。function H,W=asinc(w,L) %定义asinc函数N=fix(L);if(N<0) error('Pl

13、ease input a positive integer.(L>0)');end%if(sin(w/2)=0)%endW=(2*pi/N)*0:(N-1)'H=sin(w*N/2)/sin(w/2); 3、 指数信号对于信号xn=(0.9)n*un，使用freqz函数计算其DTFT X()。a、对w在区间-，）上绘出幅值与相位特性。这需要从freqz返回的X,W向量的移位。解释为什么幅度特性是w的偶函数，而相位特性是w的奇函数。实验程序清单1：N=100;a=1,-0.9;b=1;H,W=freqz(b,a,N);W=-pi:0.1:pi;H=freqz(b,a,W);

14、subplot(211);%把绘图窗口分成二行一列二个子图，然后在第一块区域作图plot(W,abs(H);grid,title('MAGNITUDE RESPONSE')xlabel('FREQUENCY W'),ylabel('|H(w)|')subplot(212);%把绘图窗口分成二行一列二个子图，然后在第二块区域作图plot(W,angle(H);grid,title('PHASE RESPONSE')xlabel('FREQUENCY W'),ylabel('DEGREES')实验仿真图

15、：实验结果分析：将该指数信号代入分析方程X()=中，计算化简后可从结果X()看出：该结果的幅值关于w是偶函数，相位关于w是奇函数。b、推算一阶系统的幅度特性与相位特性的表达式。设一阶差分方程：yn=xn+ayn-1,其系统方程H()=其幅度相位c、直接以这些表达式来计算幅度特性与相位特性，并与用freqz函数计算出的结果相对比。 xn=un的幅度特性为 相位特性为实验程序清单：N=1000;a=0.9;w=(2*pi/N)*-(N-1)/2:(N-1)/2; m=1./(1-a*cos(w)+(a*sin(w).2);n=1./(1-a*cos(w).2+(a*sin(w).2);H=sqrt

16、(m.2+n.2);%H的表达式subplot(2,1,1);%把绘图窗口分成二行一列二个子图，然后在第一块区域作图plot(w,sqrt(H);grid;xlabel(' W'),ylabel('|H(w)|');phase=(a*sin(w)./(1-a*cos(w);subplot(2,1,2);plot(w,atan(-phase);%绘制连续图形grid,title('PHASE ');xlabel(' W'),ylabel('DEGREES');实验仿真图：实验结果分析：以上两个图对比，可以发现用表达式

17、来计算幅度特性和相位特性这种方式与用freqz函数计算这种方式，它们最终得出的结果是一样的。3、复指数信号：a、取z=0.95ej3pi/11,对0,30，绘出xn=zn*un。用subplot指令将实部和虚部对n绘出双子图。实验程序清单：z0=0.95*exp(j*3*pi/11);n=0:30;x=z0.n;subplot(311);%用subplot指令绘出子图stem(n,x); %以x为横轴，n为纵轴绘制离散信号波形title('复指数信号');xlabel('n'),ylabel('x');subplot(3,1,2); stem(n

18、,real(x);title('REAL PART');xlabel('n'),ylabel('real');subplot(313);stem(n,imag(x); %以n为横轴，x的虚部为纵轴绘制离散信号波形title('IMAG PART');xlabel('n'),ylabel('imag');实验仿真图：b、再取一次z=0.95ej3pi/11，计算其DTFT，并对w绘出幅度。注意作为w的函数，幅度响应的尖峰位于何处。把尖峰位置与z的极坐标形式联系起来。实验程序清单：b=1;a=1,-0.

19、95*exp(j*3*pi/11);w=-(8/11)*pi:0.05*pi:(14/11)*pi;H=freqz(b,a,w);stem(w,abs(H); %以w为横轴，H的幅度为纵轴绘制离散图形grid;xlabel('w'),ylabel('|H(w)|');实验仿真图：c、如果z的角度变为3/5，粗略绘出预期的DTFT，并对w绘出幅度，并通过对freqz的计算进行绘图来验证。实验程序清单：w=-0.4*pi:0.05*pi:1.6*pi; %W的表达式z0=0.95*exp(j*3*pi/5);%Z0的表达式H=1./(1-z0*exp(-j.*w);

20、%H的表达式subplot(211);stem(w,abs(H);grid;xlabel('w'),ylabel('|H(w)|');b=1;a=1,-0.95*exp(j*3*pi/5);c=-0.4*pi:0.05*pi:1.6*pi;H=freqz(b,a,c); %调用freqz函数subplot(212); %绘制子图stem(c,abs(H); %绘制离散图形grid;xlabel('c'),ylabel('|H(c)|');实验仿真图：d、试建立一个联系带宽与r的简单表达式。实验程序清单：a1=1,-0.975*exp(j*3*pi/5);a2=1,-0.95*exp(j*3*pi/5);a3=1,-0.9*exp(j*3*pi/5);a4=1,-0.8*exp(j*3*pi/