利用DFT近似计算信号频谱专题研讨

1、数字信号处理课程研究性学习报告DFT近似计算信号频谱专题研讨姓名 学号 同组成员 指导教师 薛 健 时间 2013年5月17日 利用DFT近似计算信号频谱专题研讨【目的】(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法；(2) 理解误差产生的原因及减小误差的方法；(3)研究用DFT近似计算连续周期信号的方法；(4) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨内容】 基本题问题一已知某离散序列为 (1)用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；(2)对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱

2、峰的频率；(3)讨论所获得的结果，从中你能得到了什么结论？该结论对序列的频谱计算有何指导意义？【题目分析】本题讨论补零对离散序列频谱计算的影响。 补零可以使DFT计算得出的频谱更加细致，但是不能改变序列的DTFT【温磬提示】在计算离散非周期序列频谱时常用W/p作为横坐标，称W/p为归一化频率(normalized frequency)。在画频谱时需给出横坐标。每幅图下都需给出简要的文字说明。由于离散非周期序列频谱是周期的，所以在计算时不必用fftshift 函数对fft计算的结果进行重新排列。【序列频谱计算的基本方法】在MATLAB中，用函数fft（x，N）可以计算Xk序列的N点DFT【仿真结

3、果】L=32:L=64:L=128:L=256:L=512L=32:L=128:【结果分析】对长度为32点的序列xk进行32点的DFT,只可以获得其频谱函数X很少的细节，谱峰的频率为0.1Hz和0.9Hz。我们知道，理论上在0.2处和-0.2处有谱峰，但是由于周期性，所以所以在0,2出的频谱如图所示。对长度为32点的序列xk进行64点的DFT,相比32点的DFT，其频谱函数获得了X很多的细节；对长度为32点的序列xk进行128点的DFT,相比64点的DFT，其频谱函数获得了X更多的细节。 通过对该序列不断的补零，我们从图中可以看出补零越多，则频谱间隔越减小，频谱更清晰，频谱函数越详细，细节越多

4、，越精确，频谱的显示分辨率越大。而且不管序列补零到多少点，频谱的峰值处的对应点都是0.2和1.8，所以，序列补零并不能改变频谱的值，只能提高频谱分辨率。【自主学习内容】利用DFT分析连续非周期信号的频谱 栏栅现象【阅读文献】数字信号处理（第三版） 陈后金主编【发现问题】(专题研讨或相关知识点学习中发现的问题)：在序列后面补零的作用【问题探究】在序列后面补零不能提高频谱的分辨率。DFT是对信号fourier变换的离散化处理，从后两个图中可以明显地看出：对序列后面补零，只是增加了信号fourier变换后的离散抽样点，并不能改变信号本身的采样点，故不能提高频谱的分辨率。在连续信号在离散化或时域加窗过

5、程中，由于混叠或泄露等过程已经造成信号频谱中信息的失真，则无论怎么补零也无法再恢复已损失的信息。【仿真程序】k=0:31x=sin(0.2*pi*k)L=0:31X=fft(x,32)plot(2*L/32,abs(X)title('L=32') k=0:31x=sin(0.2*pi*k)L=0:63X=fft(x,64)plot(2*L/64,abs(X)title('L=64') k=0:31x=sin(0.2*pi*k)L=0:127X=fft(x,128)plot(2*L/128,abs(X)title('L=128') k=0:31x=

6、sin(0.2*pi*k)L=0:255X=fft(x,256)plot(2*L/256,abs(X)title('L=256') k=0:31x=sin(0.2*pi*k)L=0:511X=fft(x,512)plot(2*L/512,abs(X)title('L=512') k=0:31;x=sin(0.2*pi*k)L=0:1023X=fft(x,1024)plot(2*L/1024,abs(X),'b')hold on;k1=0:31X1=fft(x,32)stem(2*k1/32,abs(X1),'x','r&#

7、39;)title('L=32') k=0:31;x=sin(0.2*pi*k)L=0:1023X=fft(x,1024)plot(2*L/1024,abs(X),'b')hold on;k1=0:127X1=fft(x,128)stem(2*k1/128,abs(X1),'x','r')title('L=128')问题二 某离散序列为 xk=AcosW0k+Bcos ( (W0+DW)k)。用长度N=64的哈明窗对信号截短后近似计算其频谱。试用不同的A和B的值(如 A和B近似相等，A和B近差距较大)，确定用哈明窗

8、能分辩的最小的谱峰间隔中c的值。【题目分析】本题讨论用哈明窗计算序列频谱时的频率分辨率问题。就是针对不同的AB值，调整频率的差值来比较哈明窗所分辨的最小的谱峰。利用控制变量的思想，首先固定AB值，比较不同频率所造成的影响。然后固定频率，比较不同AB值所造成的影响。解题思路： Xk加窗截断，得到Xnk，对其进行DFT计算，然后画出频谱。【仿真结果】（1） AB近似的情况(A=B=1) c=0.5 c=1 C=1.5 c=2 C=2.5 c=3 C=4（2） AB相差较大的情况(A=5;B=1)【结果分析】将实验结果与教材中定义的哈明窗有效宽度相比较，发表你的看法。（1） 当AB近似时，N=3时就

9、能有效的区分能分辨出谱峰间隔。C=2是教材中定义的哈明窗的有效宽度，在实际中有时候是区分不开的。从图像可以看出c>=3的时候就能完全分辨出来,c值越大，显示的频谱越接近真实值。c=4是哈明窗的主瓣宽度，能够完全区别谱峰间隔。（2） 当AB相差较大的时候，采用相同的分析方法。我们可以确定大约c>=3能够有效的分辨出谱峰间隔。【自主学习内容】几种常见的窗函数的特点与应用。有关DFT频谱分析，对连续离散信号的分析处理，频谱分辨率，谱峰间隔的知识。【阅读文献】数字信号处理高等教育出版社 陈后金主编，薛健胡键编著【发现问题】 (专题研讨或相关知识点学习中发现的问题)：通过仿真我们发现一个很奇

10、怪的问题，那就是c=1的时候，我们看到的仿真结果好像也能分开谱峰间隔。这有悖了c值越小，分辨越差的理论特点。我们认为这种情况是由于DFT计算频谱时候误差，计算机把一些离散的点相连得到连续的点，可是离散的点的位置不同，可能会影响连续的点的取值，从而造成较大的误差，形成了这种的现象，看上去能够分辨出来谱峰间隔实际上不能分辨出来。【问题探究】在离散序列频谱计算中为何要用窗函数？用不同的窗函数对计算结果有何影响？与矩形窗相比哈明窗有何特点?如何选择窗函数？我们要用窗函数队对无限长的信号截断，不同的船函数有不同的影响，比如哈明窗与矩形窗相比，主瓣宽度加大加上旁瓣泄露，就是牺牲频谱的分辨率减少频谱的泄漏。

11、因此，我们在加窗截断时候一定要采用合理的窗函数。【仿真程序】N=64;k=0:63;L=0:511;A=1;B=1;dw1=pi/32;x1=A*cos(pi*k./4)+B*cos(pi*k./4+dw1*k);wh=(hamming(N)'x1=x1.*wh;y1=fft(x1,512);figure;plot(L/512,abs(y1);xlabel('归一化频率');ylabel('幅值');title('dw1=pi/32；C=1');dw2=pi/16;x2=A*cos(pi*k./4)+B*cos(pi*k./4+dw2*k

12、);wh=(hamming(N)'x2=x2.*wh;y2=fft(x2,512);figure;plot(L/512,abs(y2);xlabel('归一化频率');ylabel('幅值');title('dw2=pi/16;C=2');dw3=3*pi/32;x3=A*cos(pi*k./4)+B*cos(pi*k./4+dw3*k);wh=(hamming(N)'x3=x3.*wh;y3=fft(x3,512);figure;plot(L/512,abs(y3);xlabel('归一化频率');ylabel(

13、'幅值');title('dw3=3*pi/32;C=3');dw4=pi/8;x4=A*cos(pi*k./4)+B*cos(pi*k./4+dw4*k);wh=(hamming(N)'x4=x4.*wh;y4=fft(x4,512);figure;plot(L/512,abs(y4);xlabel('归一化频率');ylabel('幅值');title('dw4=pi/8;C=4');问题三已知一离散序列为 xk=cos(W0k)+0.75cos(W1k), 0£ k £ 63 其中

14、W0=0.4p, W1=W0+p/64 (1) 对xk做64点FFT, 画出此时信号的频谱。(2) 如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的64点信号补零而分辨出两个谱峰。通过编程进行证实，并解释其原因 。(3) 给出一种能分辨出信号中两个谱峰的计算方案，并进行仿真实验。【题目分析】分析影响谱峰分辨率的主要因数，进一步认识补零在在频谱计算中的作用。 影响谱峰分辨率的主要因数是序列的频率，即连续时间的抽样点数。【仿真结果】（1）（2）L=128L=256L=512L=1024（3）64点信号改为128点信号,即能分辨出来L=128L=256L=51264点信号改为256点信号，即

15、能分辨出来L=256L=512L=1024【结果分析】由结果可知，对信号增加零点不能提高频谱的分辨率；为了提高频谱的分辨率，可以增加序列的抽样点数，比如开始的时候，k的范围是0 64，后来k的范围是0 128,然后是0 256随着抽样点数的增加，频谱的分辨率逐渐提高。【自主学习内容】如何增加频谱分辨率【阅读文献】 数字信号处理（第二版）陈后金 主编。【发现问题】 (专题研讨或相关知识点学习中发现的问题)：无【问题探究】问题一、二、三讨论的是离散信号频谱的计算问题。与连续信号频谱计算问题相比较，其计算误差有何不同？ 离散信号的频谱具有周期性，会因抽样频率的不适而产生混叠，导致失真。【仿真程序】(

16、1) k=0:63;L=64;f1=0.4*pi;f2=f1+pi/64;x=cos(f1*k)+0.75*cos(f2*k);x=x zeros(1,L-length(x);X=fftshift(fft(x);f=-pi:2*pi/L:pi-2*pi/L;plot(f/pi,abs(X)grid(2)L=128;k=0:63;f1=0.4*pi;f2=f1+pi/64;x=cos(f1*k)+0.75*cos(f2*k);x=x zeros(1,L-length(x);X=fftshift(fft(x);f=-pi:2*pi/L:pi-2*pi/L;plot(f/pi,abs(X)grid(

17、3)k=0:127;L=128; f1=0.4*pi;f2=f1+pi/64;x=cos(f1*k)+0.75*cos(f2*k);x=x zeros(1,L-length(x);X=fft(x);f=-pi:2*pi/L:pi-2*pi/L;plot(f/pi,abs(X);title('128点DFT');gridk=0:255;L=256; f1=0.4*pi;f2=f1+pi/64;x=cos(f1*k)+0.75*cos(f2*k);x=x zeros(1,L-length(x);X=fft(x);f=-pi:2*pi/L:pi-2*pi/L;plot(f/pi,ab

18、s(X);title('128点DFT');grid问题四 试用DFT近似计算高斯信号的频谱抽样值。高斯信号频谱的理论值为通过与理论值比较，讨论信号的时域截取长度和抽样频率对计算误差的影响。【题目分析】学会连续非周期信号频谱计算的基本方法。分析计算中出现误差的主要原因及减小误差的方法。连续非周期信号频谱与离散Fourier变换的关系：在计算中会出现3个问题：1、混叠现象，当连续信号不是带限信号时，在连续信号离散化时，就会出现信号的频谱混叠 。在用DFT分析连续信号频谱时，信号抽样频率f对分析时的精度影响较大，因为其直接影响频谱的混叠程度，选用不同的抽样频率即可观察。2、泄漏现象

19、，由于矩形窗在其两个端点突然截断，使得频谱中存在许多高频分量。而增加信号的长度不能减少频率泄漏，而且信号的谱峰宽度与信号的长度成反比。3、栅栏现象，由于DFT得出的是离散数列，无法反映连续信号频谱中抽样点之间的细节。为了改善这种现象，需要增加DET的点数。本题以d=0.05为例，讨论fs、N变化带来的影响；再讨论d不同带来的影响。【温磬提示】信号关于t=0偶对称，在计算过程中一定要注意这个问题。【仿真结果】N变化：fs=1,N=4,d=0.05fs=1,N=6,d=0.05fs=1,N=8,d=0.05fs=1,N=16,d=0.05fs=1,N=32,d=0.05fs变化 ：fs=0.1;N

20、=128; d=0.05fs=0.15;N=128; d=0.05fs=0.2;N=128; d=0.05fs=0.3;N=128; d=0.05fs=0.4;N=128; d=0.05d变化时：fs=1;N=128; d=0.05fs=1;N=128; d=0.5fs=1;N=128; d=1fs=1;N=128; d=2fs=1;N=128; d=3【结果分析】由于信号在时域和频谱都有理论表达式，在进行误差分析时希望给出一些定量的结果。例如，可参考概率论中的对高斯概率密度函数的讨论，定义信号时域的有效宽度，频域的有效宽度，分析时域、频域有效宽度对计算误差的影响。当N逐渐变大时，谱线逐渐变为

21、谱峰，而且谱峰宽度随N在一定范围内变化时逐渐变窄；当f变化时，对DFT分析信号频谱的精度影响较大，因为其直接影响频谱混叠的精度；当d逐渐变大时，N合适大小时，需要改变抽样频率来使DFT的近似值逼近理论曲线。【自主学习内容】利用DFT分析连续飞周期信号的频谱的相关概念及分析过程中出现的现象【阅读文献】MATLAB数字信号处理数字信号处理【发现问题】 (专题研讨或相关知识点学习中发现的问题)：信号关于t=0偶对称，在计算过程中如何处理这个问题。【问题探究】N、fs、d的取值问题【仿真程序】fs=1;N=4;T=1/fs;ws=2*pi*fs;d=0.05;L=512;g=exp(-d*(-N+1:

22、N-1)*T.*(-N+1:N-1)*T);X=fftshift(fft(g,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);wt=linspace(-2*pi,2*pi,2048);G=1/T*sqrt(pi/d).*exp(-wt.*wt)/(4*d);plot(w,abs(X),wt/(2*pi),abs(G);legend('近似值','理论值');扩展题问题五本题研究连续周期信号频谱的近似计算问题。 周期为T0的连续时间周期信号x(t)可用Fourier级数表示为其中X(nw0)称为连续时间周期信号x(t)的频谱函数。称为信号的基频。

23、如果信号x(t)函数表达式已知，则可由积分得出信号的频谱。如果信号x(t)函数表达式未知，或者x(t)函数表达式非常复杂，则很难由积分得信号的频谱。本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。(1)若在信号x(t)的一个周期T0内抽样N个点，即， T为抽样周期(间隔)，可获得序列xk试推导序列xk的DFT与连续时间周期信号x(t)的频谱X(nw0)的关系；(2)由(1)的结论，给出由DFT近似计算周期信号频谱X(nw0)的方案；(3)周期信号x(t)的周期T0=1，x(t)在区间0，1的表达式为x(t)=20t2(1-t)4cos(12pt)(a)试画出信号x(t)在区间0

24、，1的波形；(b)若要用6次以内的谐波近似表示x(t)，试给出计算方案，并计算出近似表示的误差。讨论出现误差的原因及减小误差的方法。 【题目分析】根据“时域抽样，频域周期化”可知，xk的频谱xm是X(nw0)以为周期周期化后的结果。为了利用DFT近似计算连续周期信号的频谱，可先计算出序列DFT与连续周期信号的频谱X(nw0)之间的关系。然后进一步利用DFT近似计算连续周期信号的频谱。【理论推导】DFT计算所得结果Xm与连续周期信号频谱X(nw0)的关系。xk相当于是对x(t)进行了抽样，由时域离散化频域周期化可知，Xm为X(nw0)的周期化。理论推导：DFT计算所得结果Xm与连续周期信号频谱X

25、(nw0)的关系。 由于于是得 令n=m + r*N; m=0,1,2,N-1，r为整数，上式可化为化简可得对比IDFT 可得 【计算方案】把连续周期序列的一个周期进行N点抽样，并求该N点的DFT，所得DFT结果的前N点即为所求连续周期序列的频谱。主要误差原因：混叠误差从上面的关系式可以看出Xm是连续周期信号频谱X(nw0)发生周期化，再乘上比例因子N所得，由于对于一般的连续周期信号的频谱，当n±时，X(nw0)0，所以用Xm来近似为较低的谐波分量，当然这一过程中的误差很显然主要在于低次谐波分量与高次谐波 分量发生的混叠。如果在一个周期内取较多 的抽样点N，可以提高与低频分量发生混叠

26、的高频分量的频率值，高频分量如果足够高，也就是幅度值可以趋近于零，就可以进行有效的近似，其混叠误差也会有效地减少。同时，在进行信号的重建时，谐波分量选取的次数的多少对误差也有影响。对于一般的周期信号，通常时域是无限的，而我们只能选择次数较低的谐波分量（占有原信号的主要能量）进行信号重建，这是不可避免的会发生泄漏，产生误差。【扩展分析】如果周期信号x(t)是带限信号，即信号的最高频率分量为Mw0(是正整数)，试确定在一个周期内的最少抽样点N，使得在频谱的计算过程当中不存在混叠误差。与抽样定理给出的结论比较，发表你的看法。 若不发生混叠，则抽样频率fs>=Mw0，又T=T0/N,则N>=2 T0Mw0，与抽样定理基本相同，当取到等号时，若边界值不为0，则会产生混叠，故不能取到等号，若边界值为0，则可取到等号。【仿真结果】x(t)的波形：N=32，取11次谐波： N=16，取11次谐波： N=32，取4次谐波：N=32，取5次谐波：N=32，取6次谐波： N=32，取8次谐波： 【结果分析】讨论DFT点数对近似计算的影响，讨论所取谐波项的多少对近似计算的影响。误差分析要给出定量的结果，如平均误差，最大误差等。与连续非周期信号频谱计算过程中存在的误差相比较，连续周期信号频谱的计算计算误