数字信号处理实验报告3 DSP信号与系统实验报告 信号加窗及谱分析 电子科技大学 2018版.docx

《数字信号处理实验报告3 DSP信号与系统实验报告 信号加窗及谱分析 电子科技大学 2018版.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验报告3 DSP信号与系统实验报告 信号加窗及谱分析 电子科技大学 2018版.docx（18页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、优质文本电 子 科 技 大 学实 验 报 告学生姓名：Nickel 学号：20XXXXXXXXX 指导教师：杨錬一、实验室名称：数字信号处理实验室二、实验工程名称：信号加窗及谱分析三、实验原理：1、信号的时域加窗自然界的信号大多是无限长的随时间无限延伸，而实际的数字信号处理系统只能处理有限长的信号，所以在对它们进行处理之前，必须对输入信号进行分段，一段段放入系统中进行处理。具体做法是从信号中截取一个时间片段，然后用截取的信号时间片段进行周期延拓处理，得到虚拟的无限长的信号，然后就可以对信号进行傅里叶变换、相关分析等数学处理。其中对信号分段的过程称为“时域加窗。时域加窗的实质为 3.1其中，为分

2、段后的有限长信号，为原始的无限长或很长的信号，为窗函数。1.1 时域加窗对信号频域的改变时域加窗后，根据DTFT的时域相乘频域相卷积性质，变换信号的频域上表现为周期卷积，即 3.2这种卷积在一定程度上，会改变信号原频谱的特性，图3.1给出了理想低通滤波器在时域发生截断，频谱的卷积过程。图3.1 理想低通滤波的频域卷积过程1.2 窗的类型通常，我们用得最多的是矩形窗如上面例如中采用的窗，矩形窗就好似我们屋子里的窗口一样，直接对你想观察的数据进行截取。实际的信号处理过程中，矩形窗会在其边缘处突然将信号截断，窗外时域信息全部消失，导致在频域增加了频率分量，即频谱泄漏如图3.1所示，理想的低通滤波器频

3、谱中通带内和阻带内由于周期卷积产生了其他频率成分。防止泄漏的最正确方法是满足整周期采样条件，但实际中是不可能做到的。对于非整周期采样的情况，必须考虑如何减少加窗时造成的泄漏误差，主要的措施是使用合理的加窗函数，使信号截断时的锐角钝化，从而使频谱的扩散减到最少。频谱泄漏与窗函数频谱的两侧旁瓣有关，如果两侧瓣的高度趋于零，而使能量相对集中在主瓣，就可以较为接近于真实的频谱，为此，在时域中可采用不同的窗函数来截断信号。因此在矩形窗的根底上，派生出了一系列窗，最常见的就是滚降余弦窗，包括Hann窗，Hamming窗，Blackman窗，此外还有Bartlett窗和可调窗如Dolph-Chebyshev

4、窗和Kaiser窗等。常见的几种窗频谱特性如图3.2所示。 a矩形窗 bBartlett窗 cHann窗 dHamming窗 dBlackman窗图3.2 几种常见的定窗的频谱特性几种窗的时域函数分别为（1） 矩形窗 3.3（2） Bartlett窗 3.4（3） Hann窗 3.5（4） Hamming窗 3.6（5） Blackman窗 3.72、 信号的谱分析基于计算机对信号的谱分析常常采用的工具是FFT。在FFT中，要求时域序列是有限长的，所以会进行信号截断，而信号的截断产生了能量泄漏，而用FFT算法计算频谱又会产生栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同

5、的窗函数对它们的影响进行抑制。2.1 FFT算法思想：1DFT的定义：对于有限长离散数字信号xn，0 n N-1，其离散谱xk可以由离散付氏变换DFT求得。DFT和IDFT的定义为： 3.8 3.9通常令，称为旋转因子。2直接计算DFT的问题及FFT的根本思想：由DFT的定义可以看出，在xn为复数序列的情况下，完全直接运算N点DFT需要N2次复数乘法和N(N-1)次复数加法。因此，对于一些相当大的N值如1024来说，直接计算它的DFT所作的计算量是很大的。FFT的根本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以简单地组合起来得到原序列的DFT。例如，假设N为偶数，将原有的

6、N点序列分成两个N/2点序列，那么计算N点DFT将只需要约(N/2)2 2=N2/2次复数乘法。即比直接计算少作一半乘法。因子N/22表示直接计算N/2点DFT所需要的乘法次数，而乘数2代表必须完成两个N/2DFT。上述处理方法可以反复使用，即N/2点的DFT计算也可以化成两个N/4点的DFT假定N/2为偶数，从而又少作一半的乘法。这样一级一级地划分下去一直到最后就划分成两点DFT运算的情况。a基2按时间抽取DIT的FFT算法思想：设序列长度为，L为整数如果序列长度不满足此条件，可通过在后面补零让其满足。将长度为的序列，先按n的奇偶分成两组： DFT化为： 3.10上式中利用了旋转因子的可约性

7、，即：。又令 ， 那么式3.10可以写成： 3.11可以看出，分别为从中取出的N/2点偶数点和N/2点奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从式3.11可以看出，这样的组合仅表示出了前N/2点的DFT值，还需要继续利用表示的后半段本算法推导才完整。利用旋转因子的周期性，有：，那么后半段的DFT值表达式： 同样， k=0,1,N/2-1，所以后半段k=N/2,N-1的DFT值可以用前半段k值表达式获得，中间还利用到，得到后半段的值表达式为： 3.12 这样，通过计算两个N/2点序列的N/2点DFT，可以组合得到N点序列的DFT值，其

8、组合过程如图3.3所示：图3.3 两个N/2点DFT组合成一个N点DFT比方，一个N = 8点的FFT运算按照这种方法来计算FFT可以用图3.4来表示：图3.4 8点基2按时间抽取DIT算法流程图b基2按频率抽取DIF的FFT算法思想：设序列长度为，L为整数如果序列长度不满足此条件，通过在后面补零让其满足。在把按k的奇偶分组之前，把输入按n的顺序分成前后两半：因为，那么有，所以：按k的奇偶来讨论，k为偶数时：k为奇数时：前面已经推导过，所以上面的两个等式可以写为：通过上面的推导，的偶数点值和奇数点值分别可以由组合而成的N/2点的序列来求得，其中偶数点值为输入xn的前半段和后半段之和序列的N/2

9、点DFT值，奇数点值为输入xn的前半段和后半段之差再与相乘序列的N/2点DFT值。 令，那么有：这样，也可以用两个N/2点DFT来组合成一个N点DFT。3在FFT计算中使用到的MATLAB命令：函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)那么计算R点序列的N点DFT，假设RN，那么直接截取R点DFT的前N点，假设RN，那么x先进行补零扩展为N点序列再求N点DFT。函数ifft(X)可以计算R点的谱序列的R点IDFT值；而ifft(X,N)同fft(x,N)的情况。四、实验目的：1、理解信号时域加窗后频谱的所发生的变化。2、观察时域选择加不同的窗频域的变化情况，从而分析各种

10、窗的特性。3、深刻理解FFT根本思想，学习基2 DIT和基2 DIF算法。重点理解FFT算法在计算DFT方面的时间优化。4、观察时域选择加不同的窗后利用FFT绘出的该信号频域的变化情况，再次确认各种窗的特性。五、实验内容：1、计算两个给定的加窗的实数序列的DTFT，其中窗函数分别取为矩形窗，Bartlett窗，Hann窗，Hamming窗和Blackman窗，比照各自的频谱特性。2、对以上两个序列，利用DFT定义式，编程直接计算2个要求序列的DFT值直接加矩形窗对其截断。3、利用MATLAB中提供的FFT函数，计算以上两个序列的DFT值，比照与直接利用DFT定义式计算的结果与所花费的计算时间。

11、4、对以上两个序列，分别采用不同的窗Bartlett窗，Hann窗，Hamming窗和Blackman窗，利用FFT函数计算DFT值，比照其结果。六、实验器材设备、元器件：安装MATLAB软件的PC机一台，DSP实验演示系统一套。七、实验步骤：1、计算实数序列的DTFT，其中分别取为矩形窗，Bartlett窗，Hann窗，Hamming窗和Blackman窗，比照各自的频谱特性。2、计算周期为1kHz的方波序列占空比为50，幅度取为/-512，采样频率为25kHz，取256点长度的DTFT，对其进行截断分别取矩形窗，Bartlett窗，Hann窗，Hamming窗和Blackman窗，比照各自

12、的频谱特性。3、利用DFT定义式，编程直接计算上面两个序列的DFT值。4、利用MATLAB中提供的FFT函数，计算2个要求序列的DFT值。调用MATLAB中提供的计算程序运行时间函数，比照与直接利用DFT定义式计算所花费的时间。5、对以上两个序列，分别采用不同的窗Bartlett窗，Hann窗，Hamming窗和Blackman窗，利用FFT函数计算DFT值，比照其结果。6、(拓展要求)不改变序列的点数，仅改变DFT计算点数如变为计算1024点DFT值，观察画出来的频谱与前面频谱的差异，并解释这种差异。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。八、实验数据及结果分析：需

13、要的程序代码：（1） 对要求的2个序列加不同窗后计算DTFT，比照其各自的频谱。hold ongrid onn = 0:255;Fs = 1;L = length(n);x = cos(5*pi*n/16);w1 = bartlett(L);w2 = hann(L);w3 = hamming(L);w4 = blackman(L);LINES = 5; ROWS = 2; ID = -1; STEP = 2;ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w1)ID =

14、 ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w2)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w3)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w4)ID = 0;ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x.*w1,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x.*w2,L,Fs)ID =

15、 ID+STEP; subplot(LINES,ROWS,ID)dtft(x.*w3,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x.*w4,L,Fs)hold ongrid on n = 0:255;Fs = 25000;L = length(n);x = square(2*pi*1000*n/Fs,50); w1 = bartlett(L);w2 = hann(L);w3 = hamming(L);w4 = blackman(L); LINES = 5; ROWS = 2; ID = -1; STEP = 2;ID = ID+STEP; su

16、bplot(LINES,ROWS,ID)stem(x)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w1)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w2)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w3)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w4) ID = 0;ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x,L,Fs)ID = ID+STEP; subplot(LINES,R

17、OWS,ID)dtft(x.*w1,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x.*w2,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x.*w3,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)dtft(x.*w4,L,Fs)（2） 对要求的2个序列直接用定义式进行DFT计算的程序。以下是自己写的dft函数，不调用MATLAB的fft函数：% Author = Nie H.T.% Stu ID = 2016010903009% Date = Mar 27, 201

18、8% Calculate DFT without internal function fft.%function X = dft(x)L = length(x);X = zeros(L);for k = 1:L for n = 1:L X(k) = X(k)+x(n)*exp(-1i*2*pi*k*n/L); endendn = 0:255;L = length(n); Fs1 = 1;x1 = cos(5*pi*n/16); Fs2 = 25000;x2 = square(2*pi*1000*n/Fs2,50); ticX1 = dft(x1);tocticX2 = dft(x2);toc

19、LINES = 2; ROWS = 1; ID = 0; STEP = 1;ID = ID+STEP; subplot(LINES,ROWS,ID)stem(X1)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(X2)（3） 利用MATLAB中提供的FFT函数，计算2个要求序列的DFT值，并调用MATLAB中提供的计算程序运行时间函数，比照与直接利用DFT定义式计算所花费的时间的程序。直接调用理论上和上面的效果相同，只要把dft改成fft即可。但是时间上差距较大：代码如下：n = 0:255;L = length(n); Fs1 = 1;x1 = cos(5*

20、pi*n/16); Fs2 = 25000;x2 = square(2*pi*1000*n/Fs2,50); ticX1 = fft(x1);tocticX2 = fft(x2);toc LINES = 2; ROWS = 1; ID = 0; STEP = 1;ID = ID+STEP; subplot(LINES,ROWS,ID)stem(X1)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(X2)（4） 对要求的2个序列，分别采用不同的窗Bartlett窗，Hann窗，Hamming窗和Blackman窗，利用FFT函数计算DFT值，比照其结果的程序。直

21、接调用理论上和上面的效果相同，只要把dft改成fft即可。hold ongrid on n = 0:255;Fs = 1;L = length(n);x = cos(5*pi*n/16); w1 = bartlett(L);w2 = hann(L);w3 = hamming(L);w4 = blackman(L); LINES = 5; ROWS = 2; ID = -1; STEP = 2;ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w1)ID = ID+STEP

22、; subplot(LINES,ROWS,ID)stem(x.*w2)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w3)ID = ID+STEP; subplot(LINES,ROWS,ID)stem(x.*w4) ID = 0;ID = ID+STEP; subplot(LINES,ROWS,ID)fdtft(x,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)fdtft(x.*w1,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)fdtft(x.*w2,L,Fs)ID = ID+

23、STEP; subplot(LINES,ROWS,ID)fdtft(x.*w3,L,Fs)ID = ID+STEP; subplot(LINES,ROWS,ID)fdtft(x.*w4,L,Fs)5拓展要求分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。直接调用理论上和上面的效果相同，只要把dft改成fft即可。但是时间上差距较大。另外，fft函数会对序列自动补零，会使DFT结果更加细致平滑。6拓展要求保持序列点数不变，改变DFT计算点数变为1024点，观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。左边是原始序列长度256，右边是变为1024

24、点序列之后的DTFT效果，可见时域补零相当于频域插值，从而对频谱分辨率有着提高的效果。 九、 实验结论：直接调用fft函数结果和自己编写的dft函数效果相同，但是时间上差距较大。另外，fft函数会对序列自动补零，会使DFT结果更加细致平滑。时域补零相当于频域插值，从而对频谱分辨率有着提高的效果。十、 总结及心得体会：MATLAB为我们提供了高效的函数用以进行数字信号处理的根本操作，如fft等。这些函数都可以在Mathworks官网找到。十一、对本实验过程及方法、手段的改良建议：可以多介绍一些Simulink以及dsp.类函数的使用。使用过程中发现此类函数可以大大减少编程时间。 报告评分： 指导教师签字：