DSP研究性学习报告频谱计算.

1、站索X.逼乂爭数字信号处理课程研究性学习报告指导教师薛健时间2014年4月DFT近似计算信号频谱专题研讨【目的】(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法。(2) 理解误差产生的原因及减小误差的方法。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】基本题1. 已知一离散序列为xk sin(0.2 水),k 0,1,31(1) 用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱峰的频率；(3)讨论所获得的结果，给出你的结论。该结论对序列的频谱计

2、算有何指导意义？【题目分析】本题讨论补零对离散序列频谱计算的影响。【温磬提示】在计算离散非周期序列频谱时常用/作为横坐标，称/为归一化频率 normalized frequency)。在画频谱时需给出横坐标。每幅图下都需给出简要的文字说明。由于离散非周期序列频谱是周期的，所以在计算时不必用fftshift函数对fft计算的结果进行重新排列。【序列频谱计算的基本方法】(1) 周期为N的离散周期信号 k的频谱为N 1jJmkXm DFSxk ke NK 0(2) 离散非周期信号xk的频谱为X(ej ) DTFT xkxk e j k,k(3) 有限长N的序列xk的频谱为N 1j 乙 mkXm xk

3、e N ,m=0,1,N-1k 0【仿真结杲】(1)15L=321050.80.9100.10.20.30.40.50.60.7Normalized frequency(2 )2015L=128eungM10A j I*J000.10.20.30.40.50.60.70.80.91Normalized freque ncyeungM01 11kK142011511PA-15Ar：-l: :e1:u10! 1f-i-fg-=M.| 1 *5A !_/aA jt y y wrL*0A/r,Ay j v L=256L=512201510500.10.20.30.40.50.60.70.80.91No

4、rmalized freque ncy00.10.20.30.40.50.60.70.80.91Normalized freque ncy【结果分析】在有限长序列后补 0，不会增加任何信息，补0前后的两序列对应的 DTFT完全一致，但补0 后的DFT存在明显差别。从信号表示的角度来看，对于长度为N的时域序列xk,可由N点的DFT对应频域序列 Xm唯一表示，Xm是序列xk的离散时间Fourier变换X ( e)在一个周期2的等间隔抽样，由于抽样间隔不同所以Xm不同。但从信号频谱分析的角度，在序列xk后补可以在X（ej ）的一个周期2内获得更多的抽样值，获得更多信息。【自主学习内容】【阅读文献】【

5、发现问题】（专题研讨或相关知识点学习中发现的问题）：【问题探究】2已知一离散序列为 x k=Acos ok+Bcos0+）k）。用长度N=64的哈明窗对信号截短后近似计算其频谱。试用不同的A和B的值（如A和B近似相等，A和B近差距较大），确定用哈明2n窗能分辩的最小的谱峰间隔A w C今中c的值。N【题目分析】本题讨论的是 Hamming窗截断信号后信号频谱的频率分辨率。首先，加窗处理对频谱分析主要有两个方面的影响：O1频谱中出现多余的高频分量，这是窗函数突然截断而引起的，称为频率泄漏现象；O2谱线变成了有一定宽度的谱峰，谱峰的宽度与信号的长度成反比，所取信号的长度越长，谱峰的宽度就越窄，当信

6、号中两个不同频率分量的频率差？?小于谱峰宽度时，计算出的频谱可能显示不出两个明显的峰值。为使计算出的频谱能显示出相邻的谱 峰，则要求相邻频率分量的频率差?大于谱峰的有效宽度，即：?= ?勿2 ?所以，本题采用 Hamming窗截断信号在进行频谱分析。Hamming窗是以降低频率分辨率为代价来减少频率泄漏的，使用Hamming窗时要求能分辨的谱峰的间隔为：?=? /=?/?按照题目要求，选取 ？?=0.5,再分别令 A=1、B=1和A=1、B=0.1，不断改变c值，观察频谱的 变化，即可确定刚好可以分辨谱峰的c值。【仿真结果】当a=1、B=1时，c取1、2、3、4时的频谱如下图所示。其中w=?。

7、20w=0.29452 ; c=3w=0.19635 ; c=200.10.20.30.40.50.60.70.80.91归一化频率w=0.3927 ; c=42015151510值 幅2015140.40.5归一化频率可以看出，c=2时两谱峰不能分辨,w=0.29452 ; c=315c=3、c=4时可以分辨。w=0.3927 ; c=420w=?。20201510值 幅值 幅0.5归一化频率值10值1幅150.40.5归一化频率当A=1、B=0.1时，c取1、2、3、4时的频谱如下图所示。其中15elk11* 1r-1J-ii-1j w-111-0.10.20.30.40.50.60.70

8、.80.9归一化频率j 1 11 i 14 ”r1*.T.rw=0.098175 ; c=11210120.160400.20.30.40.50.60.70.8归一化频率00.9a11n110VP A-5-1 1 000.10.20.30.40.50.60.70.80.91归一化频率3J00.10.20.30.40.50.60.70.80.9w=0.19635 ; c=2归一化频率10w=0.3927 ; c=4可以看出，c=2时两谱峰不能分辨，c=3、c=4时可以分辨。【结果分析】将实验结果与教材中定义的哈明窗有效宽度相比较，发表你的看法。Hamming窗函数的幅值有中心向两端逐渐减弱，因而

9、其高频分量明显减小，频谱中旁瓣的幅度 较小，主瓣峰值与第一个旁瓣峰值相对衰减很大可达40dB, Hamming窗是以降低频率分辨率为代价来减少频率泄漏的。信号的时域加窗会导致在信号的频谱中产生频率泄漏。对同一个窗函数，增加长度N虽然可以减少主瓣宽度，但不减少旁瓣泄漏。为在主瓣宽度与旁瓣泄漏之间取得良好平衡，应根据实际信号 的特性采用合适的窗函数。35【自主学习内容】【阅读文献】【发现问题】（专题研讨或相关知识点学习中发现的问题):问题：如果使用矩形窗来对信号进行截短，信号的频谱将如何变化?解答：当A=1、B=1时，使用矩形窗截短信号的频谱如图。w=0.098175 ; c=1w=0.19635

10、 ; c=23530302525202015151050当A=1、B=0.1时，使用矩形窗截短信号的频谱如图。w=0.098175 ; c=13535302520151050w=0.29452 ; c=335302520151050值 幅值 幅值 幅值 幅0.5归一化频率0.40.5归一化频率0.5归一化频率0.5归一化频率0.5归一化频率0.40.5归一化频率阮 叩j 11 Iw=0.29452 ; c=335w=0.3927 ; c=4302520151050302520151050J1 HlL-11 1I-|11 j1 1AL-rIpJIIr 川1 T .|w=0.19635 ; c=2

11、E% 1：i30-i-25-1-20-1-15-111-10-1-5-JlI n,I. ILIl丿；.rr-01r11.0.10.20.30.40.50.60.70.80.9归一化频率值 幅w=0.3927 ; c=4可以看出，Haming窗可以有效减少旁瓣泄漏，但会造成分辨率下降。且当存在小信号也就是信 号中存在较弱的频率分量时，Haming窗比矩形窗更好。【问题探究】在离散序列频谱计算中为何要用窗函数？用不同的窗函数对计算结果有何影响？与矩形窗相比 哈明窗有何特点？如何选择窗函数？如果连续信号无限长，则离散后的序列也无限长，无法适用于DFT分析，需要对其进行加窗截短使之成为有限长序列，既可

12、以对信号进行DFT分析，也可以减少计算机的计算量，提取出我们想要的重点信息。对于不同的信号需要选择不同的窗函数，否则可能会因为窗函数选择的不合理造成较大的失真，从而无法得到符合真实情况的有用信息。与矩形窗相比，Haming窗主瓣宽度增加，频率分辨率下降，但是旁瓣的能量也随之下降，减少 了频率泄漏。但是 Hamming窗是以降低频率分辨率为代价来减少频率泄漏的。根据本题的研究结果，可知当存在小信号也就是信号中存在较弱的频率分量时，Haming窗比矩形窗更好。3 已知一离散序列为 xk=cos( ok)+0.75cos( ik), 0 k 63 其中 o=O.4 ,1= 0+ /64(1) 对xk

13、做64点FFT,画出此时信号的频谱。(2) 如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的64点信号补零而分辨出两个谱峰。通过编程进行证实，并解释其原因。(3) 给出一种能分辨出信号中两个谱峰的计算方案，并进行仿真实验。【题目分析】分析影响谱峰分辨率的主要因数，进一步认识补零在在频谱计算中的作用。本题目具体探究了关于信号分辨率的具体问题。进一步研究了关于信号处理中参数的选择问题。我们知道在进行信号的频域分析的时候需要一下4个参数1截取信号长度N2抽样频率fs3采样点数L4选择什么样的窗本题目中选择的矩形窗，在我们常用的三大类窗中(矩形窗，三角窗，指数窗)是主瓣宽度最 小的一种窗。而本

14、题目中给出的信号所需要的分辨频率/64此时c=0.51，而且本题目中两个信号的幅度差不多没有小信号所以如果修改窗函数类型没法达到预期目的。首先我们按照题目提示进行补零实验。改变采样点数L只能加大频域的取样点数补零起到对原X(K )做插值的作用一方面克服“栅栏“效应使谱的外观得到平滑另一方面由于对数据截断时所引起的频谱泄露有可能在频谱中出现一些难以确认的谱峰补零后有可能消除这种现象。但是这种方法只能改变频谱的可视分辨率。对于实际分辨率没有贡献。所以解决这个问题应该从截取信号的长度即窗的宽度下手，改变矩形窗主瓣宽度。【仿真结果】403530252015105-4011【结果分析】(1) 仿真结果如

15、图1所示，即信号频谱图。(2) 由(1)结果可知图中显示的谱不能分辨两个谱峰，进行补零后频谱图如图2所示，由图可知并不能使分辨率提高，依然分辨不出两个谱峰。(3) 解决方案：采用哈明窗，取N=64*2，仿真结果如图3所示，可见此时可清晰分辨出两个谱峰。【自主学习内容】【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：【问题探究】1、2、3题讨论的是离散信号频谱的计算问题。与连续信号频谱计算问题相比较，其计算误差 有何不同？4试用DFT近似计算高斯信号g(t) exp( dt2)的频谱抽样值。高斯信号频谱的理论值为rn2G(j) bexp(站)通过与理论值比较，讨论信号的时域截取

16、长度和抽样频率对计算误差的影响。【题目分析】连续非周期信号频谱计算的基本方法。计算中出现误差的主要原因及减小误差的方法。连续非周期信号频谱与Fourier的关系：112 n?e ?= - Xj (w nw sam) =X(j )T nT nT在计算中出现3个问题：1、混叠现象：当连续信号不是带限信号时，在连续信号离散化时，就会出现信号的频谱混叠。再用DFT分析连续信号频谱时，信号抽样频率f对分析时影响精度较大，因为其直接影响频谱的混叠程度。2、泄漏现象：由于矩形窗在其两个端点突然截断，使得频谱中存在许多高频分量，而增加信 号的长度不能减少频率泄漏，而且信号的谱峰宽度与信号的长度成正比。3、栅栏

17、现象：由于 DFT得出的是离散序列，无法反映连续信号频谱中抽样点的细节，为了改善这种现象，需要增加DET的点数。【仿真结果】1.5r10.50.50.5-0；02cuttime=1 fsam=4-20 -10 0 10 200-4-20.5-105cuttime=4 fsam=2-500-20 -10 0 10 200.5-0*-40-2002040cuttime=0.5 fsam=8cuttime=0.25 fsam=161.:0.5-40-20020400-100 -50050100cuttime=0.5 fsam=16cuttime=8 fsam=1cuttime=2 fsam=4cut

18、time=1 fsam=80.502410-20 -10-4-2-40-2002040-510 20L 1 - - 11-100 -50050100【结果分析】由于信号在时域和频谱都有理论表达式, 当截取时间增大时，谱线逐渐变为谱峰,在进行误差分析时希望给出一些定量的结果。而且谱峰宽度随截取时间在一定范围变化时内变窄。当抽样频率变化时，对DFT分析频谱的精度影响较大，以为其直接影响频谱混叠的精度。很容易能看出，增加截取时间和抽样频率都可以减小误差。【自主学习内容】连续非周期信号的频谱计算方法。【阅读文献】【发现问题】（专题研讨或相关知识点学习中发现的问题）：【问题探究】扩展题5本题研究连续周期

19、信号频谱的近似计算问题。周期为To的连续时间周期信号 x(t)可用Fourier级数表示为x(t) X(n )ejnotn其中1jn otX(n o)t0 x(t) e j odtToX(n o)称为连续时间周期信号x(t)的频谱函数。o 2 nTo 2仇称为信号的基频(基波)，n 称为信号的谐波。如果信号x(t)函数表达式已知，则可由积分得出信号的频谱。如果信号x(t)函数表达式未知，或者x(t)函数表达式非常复杂，则很难由积分得信号的频谱。本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。(1) 若在信号x(t)的一个周期To内抽样N个点，即ToNT， T为抽样周期 涧隔)，

20、可获得序列xkxk x(t)tgk 0,1, ,N 1试分析序列xk的DFT与连续时间周期信号 x(t)的频谱X(n o)的关系；(2) 由(1)的结论，给出由DFT近似计算周期信号频谱X(n o)的方案；(3) 周期信号x(t)的周期To=1，x(t)在区间0，1的表达式为x(t)=20t2(1 t)4cos(12 t)(a) 试画出信号x(t)在区间0，1的波形；(b) 若要用4次以内的谐波(包含 4次)近似表示x(t)，试给出计算方案，并计算出近似表示的误差。 讨论出现误差的原因及减小误差的方法。【题目分析】【理论推导】DFT计算所得结果Xm与连续周期信号频谱X(n o)的关系。【计算方

21、案】根据理论推导结果设计近似计算方案。分析产生误差的主要原因。【扩展分析】如果周期信号x(t)是带限信号，即信号的最高频率分量为M o(是正整数)，试确定在一个周期内的最少抽样点 N,使得在频谱的计算过程当中不存在混叠误差。与抽样定理给出的结论比较，发表 你的看法。【仿真结果】【结果分析】讨论DFT点数对近似计算的影响，讨论所取谐波项的多少对近似计算的影响。误差分析要给出定量的结果，如平均误差，最大误差等。 与连续非周期信号频谱计算过程中存在的误差相比较，连续周期信号频谱的计算计算误差有何 异同？【 自主学习内容 】阅读文献 】发现问题 】 （专题研讨或相关知识点学习中发现的问题）：问题探究

22、】【 仿真程序 】1、(1)k=0:31;x=sin(0.2*pi*k);x_32=fft(x,32);L=0:31;plot(L/32,abs(x_32);xlabel( Normalized frequency title( L=32 );(2)k=0:31;x=sin(0.2.*pi.*k);X_64=fft(x,64);L=0:63;subplot(2,2,1);plot(L/64,abs(X_64);xlabel( Normalized frequency title( L=64 ); X_128=fft(x,128);L=0:127;subplot(2,2,2);plot(L/12

23、8,abs(X_128);xlabel( Normalized frequency title( L=128 );X_256=fft(x,256);L=0:255;subplot(2,2,3);plot(L/256,abs(X_256););ylabel( Magnitude ););ylabel( Magnitude ););ylabel( Magnitude );xlabel( Normalized frequency title( L=256 );X_512=fft(x,512);L=0:511;subplot(2,2,4); plot(L/512,abs(X_512);xlabel(

24、Normalized frequency title( L=512 ););ylabel();ylabel(MagnitudeMagnitude););2、A=1;B=1;k=0:63;%e ?e ?D?0?3O ? e ?a64l=0:511;%c=input(c=); e ? e ?c?卩N=64;L=512;subplot(221);c=1;x=A*cos(0.5*pi*k)+B*cos(0.5*pi+c*2*pi/N)*k); h=(hamming(N);%33 ? e ?a64 卩?Hammingx=x.*h;%o ?Hamming ?D?o?X=fft(x,L); %?i ?u F

25、ourier 士 ? plot(l/L,abs(X);%?-3? i ?xlabel( 1 eo ? _ ?卩？ e ylabel( u ?卩)%e ?e ?D?o?w=2*c*pi/N;title(strcat(subplot(222);w= ,num2str(w),c= ,num2str(c);c=2; x=A*cos(0.5*pi*k)+B*cos(0.5*pi+c*2*pi/N)*k); h=(hamming(N);%33 ? e ?a64 卩?Hammingx=x.*h;%o ?Hamming ?D?o?X=fft(x,L); %?i ?u Fourier 士 ?%e ?e ?D?o

26、?plot(l/L,abs(X);%?-3? i ?xlabel( 1 eo ? _ ?卩? e； ylabel( u ?卩) w=2*c*pi/N;title(strcat(w= ,num2str(w),c= ,num2str(c);subplot(223);c=3;x=A*cos(0.5*pi*k)+B*cos(0.5*pi+c*2*pi/N)*k); h=(hamming(N);%33 ? e ?a64 卩?Hamming%e ?e ?D?o?x=x.*h;X=fft(x,L);%o ?Hamming ?D?o?%?i ?u Fourier 士 ?plot(l/L,abs(X);%?-3

27、? i ?xlabel(ylabel(1 eo ? _ ?卩?);U ? 1)w=2*c*pi/N;title(strcat(w= ,num2str(w), ; c=subplot(224);c=4; x=A*cos(0.5*pi*k)+B*cos(0.5*pi+c*2*pi/N)*k);h=(hamming(N);%33 ? e ?a64 卩?Hammingx=x.*h;%o ?Hamming ?D?o?X=fft(x,L);%?i ?u Fourier 士 ?plot(l/L,abs(X);%?-3? i ?xlabel( 1 eo ? _ ?卩?); ylabel( u ?卩) w=2*c*pi/N;,num2str(c);%) ?e ?D?o?title(strcat(w= ,num2str(w),c= ,num2str(c);3、(1)W0=0.4*pi;W1=W0+pi/64; k=0:63;x=cos(W0*k)+0.75*cos(W1*k);y=fft(x);yy=fftshift(abs(y); k=linspace(-pi,pi,length(yy); stem(k,yy,. );(2