电子信息工程数字谱分析实践及误差讨论

1、汕头大学工学院三级项目报告课程名称：数 字 信 号 处 理课程设计题目：数字谱分析实践及误差讨论指导教师：系别：电子工程系专业：电子信息工程学号：姓名：合作者完成时间：2010年10月6日至10月26日成绩：评阅人：一 项目意义与目标数字谱分析是DSP应用系统最基础算法。如语音压缩（MP3、无线通信）用到FFT，频谱分析仪用到FFT，4G的关键技术OFDM用到FFT在实际DSP系统中能够正确运用数字谱分析算法。通过亲自编程、绘频谱图，掌握数字谱分析有关算法DTFT、DFT和FFT，并能够恰当举例说明数字谱分析算法产生误差原因及程度，找到减少误差的方法。理解数字谱分析中物理分辨力和计算分辨力的概

2、念，讨论两者间的关系。二 项目内容1. 举例说明时窗宽度、类型与物理分辨率的关系。得出相应结论，给出提高物理分辨率的方法。2. 举例说明时窗类型与频谱泄露量的关系，说明频谱泄露的危害。给出降低频谱泄露量的方法。3. 举例说明物理分辨率与计算分辨率的关系，如何避免频率分量的漏判与误判。4. 进行课堂演示及讨论(PPT)。5. 分组讨论，独立撰写项目报告，提交Word电子文档。三 项目报告正文1. 数字谱分析的意义信号特征提取，语音、图象处理，通信信号处理等。应用范围：家用电器，仪器仪表，医疗仪器，有线通信，无线通信，生物生理信息处理，地球物理信息探测，雷达，电子对抗，航空航天，宇宙探索，动力控制

3、系统几乎所有电子设备和信息语音压缩（MP3、无线通信）用到FFT，频谱分析仪用到FFT，4G的关键技术OFDM用到FFT2. DTFT、DFT和FFT算法概述用Equation编辑器编辑公式，下同。3. 时窗效应的仿真分析3.1 矩形窗和海明窗各自与物理分辨率的关系3.1.1Rectangular Window当时窗宽度L=64时，设定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代码为：clear;close all;L=64;n=0:L-1;f1=1000;f2=1050;f3=1100;fs=5000;w1=2*pi*f1/fs;w2=2*pi*f2/fs;

4、w3=2*pi*f3/fs;xL=sin(w1*n)+sin(w2*n)+sin(w3*n);w=0:2*pi/300:2*pi;yL=dtft(xL,w);plot(w/2/pi*fs,abs(yL),'k')axis(500,2500,0,L/1.5)xlabel('f/Hz')gtext('L=64');仿真结果如下：当时窗宽度L=200时，设定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代码：（略）仿真结果如下：仿真结果分析：由物理分辨率计算公式知：（c=1），则当L=64，fs=5kHz时=78.125

5、Hz， 又因为3个信号的谱峰间隔为50Hz，则>.物理分辨率大于信号频率间隔，故不能分辨出3个信号的频谱；当L=200，fs=5kHz时，=25Hz，则<物理分辨率小于信号频率间隔，故能分辨出3个信号的频谱3.1.2 Hamming window当时窗宽度L=100时，设定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代码为：clear;close all;L=100;n=0:L-1;fs=5000;f1=1000;f2=1050;f3=1100;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*f3/fs;xL=sin(w1*n

6、)+sin(w2*n)+sin(w3*n);w=0:2*pi/500:2*pi;wh=0.54-0.46*cos(2*pi*n/(L-1);xL_H=xL.*wh;yL=dtft(xL_H,w);plot(w/2/pi*fs,abs(yL),'k')axis(500,1500,0,L/3)xlabel('f/Hz')gtext('L=100');仿真结果如下：当时窗宽度L=300时，设定f1=1kHz，f2=1.05kHz，f3=1.10kHz；fs=5kHz。代码：clear;close all;L=300;n=0:L-1;fs=5000;f1

7、=1000;f2=1050;f3=1100;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*f3/fs;xL=sin(w1*n)+sin(w2*n)+sin(w3*n);w=0:2*pi/1000:2*pi;wh=0.54-0.46*cos(2*pi*n/(L-1);xL_H=xL.*wh;yL=dtft(xL_H,w);plot(w/2/pi*fs,abs(yL),'k')axis(500,1500,0,L/3)xlabel('f/Hz')gtext('L=300');仿真结果如下：仿真结果分析：由物理分辨率计算公式知：

8、 ，（c=2）则当L=100，fs=5kHz时=100 Hz， 又因为3个信号的谱峰间隔为50Hz，则>.物理分辨率大于信号频率间隔，故不能分辨出3个信号的频谱；当L=300，fs=5kHz时，=33.3Hz，则<物理分辨率小于信号频率间隔，故能分辨出3个信号的频谱结论：取相同的时窗宽度L，则矩形窗的物理分辨率比海明窗高。由物理分辨率计算公式可知，在采样频率fs一定时，增加时窗宽度L，两种类型的时窗的物理分辨率均提高。提高物理分辨率的方法是：在采样率不变的前提下增加时窗宽度L或者使用矩形窗。3.2矩形窗和海明窗各自与频谱泄露量的关系 Rectangular leakage设时窗宽度

9、L=100时，设定f1=0.5kHz，f2=1kHz，f3=1.5kHz；fs=5kHz。代码如下：clear;close all;L=100;n=0:L-1;f1=500;f2=1000;f3=1500;fs=5000;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*f3/fs;xL=sin(w1*n)+0.1*sin(w2*n)+sin(w3*n);w=0:2*pi/500:2*pi;yL=dtft(xL,w);plot(w/2/pi*fs,abs(yL),'k')axis(0,2000,0,L/1.5)xlabel('f/Hz')g

10、text('L=100');仿真结果如下：仿真结果分析：由图可知，信号以矩形窗截取，有较大的频谱泄露。频率为0.5kHz及1.5kHz的信号频谱，最靠近谱峰的副瓣幅值比频率为1kHz的信号的谱峰值还要高，导致难以根据谱峰值区分信号，从而引致信号的缺失。Hamming leakage设时窗宽度L=200时，设定f1=0.5kHz，f2=1kHz，f3=1.5kHz；fs=5kHz。代码为：clear;close all;L=200;n=0:L-1;fs=5000;f1=500;f2=1000;f3=1500;w1=2*pi*f1/fs;w2=2*pi*f2/fs;w3=2*pi*

11、f3/fs;xL=sin(w1*n)+0.1*sin(w2*n)+sin(w3*n);w=0:2*pi/500:2*pi;wh=0.54-0.46*cos(2*pi*n/(L-1);xL_H=xL.*wh;yL=dtft(xL_H,w);plot(w/2/pi*fs,abs(yL),'k')axis(0,2000,0,L/3)xlabel('f/Hz')gtext('L=200');仿真结果如下：仿真结果分析：由图可知，信号以海明窗截取，频谱泄露大为减小。结论：在满足物理分辨率的前提下，用海明窗截取的信号的频谱泄露远比用矩形窗截取的少。频谱泄露不

12、仅会降低信号功率而且还会掩盖一些微弱信号，导致信号缺失。我们可以通过增加采样的长度或者加海明窗来降低频谱泄露量。4. DFT/DTFT计算分辨率与物理分辨率关系的仿真分析4.1 物理分辨率仿真 代码见附录。仿真结果如下： 图14.2 计算分辨率仿真：（代码见附录）仿真结果如下： 图2仿真结果分析：由上图1可知，当物理分辨率低于信号频率间隔时，无论计算分辨率再大，也无法从频谱图中把信号的频率分量分辨出来；相反，由图2知，当物理分辨率高于信号频率间隔时，随着计算分辨率的加大频谱图像取点更丰富，N=16时，3个频率分量中只有频率f=2.5KHz有对应的频谱点，很可能造成频率分量的漏判；N=64时，基

13、本可以判断有3个频率分量，除了频率f=2.5KHz能正确对应频谱，其余两个分量的频谱点与其对应频率都有偏差，这时候就很可能造成误判；N=256时，清楚显示3个频率分量及其对应频率。综上所述，在满足物理分辨率的前提下，增加计算分辨率，可以有效避免频率分量的漏判与误判。5. 数字谱分析中存在的误差及减小误差的方法Aliasing, Frequency resolution, Frequency leakage四 总结 通过撰写这次三级项目的报告，我深深体会到我们CDIO理念中的Conceive和Implement的意义。随着这次项目的开展，我不断的自我增值，独立撰写报告使我的自学研究能力有所提高，

14、和同学们的讨论在增强我的沟通及语言表达能力的同时也加强了我的团队合作精神。开始时，为了更深入了解数字谱分析，我翻阅英文教材以及研读老师的课件。为了实现数字谱分析的MATLAB仿真，我从图书馆借来MATLAB教程自学MATLAB软件，并不断向其他精通MATLAB的同学请教。从开始构思如何分析数字谱到最终通过MATLAB软件仿真来实现构想，期间我不但加深了对数字谱分析的理解，还学会了运用MATLAB来对某些理论进行仿真验证，真是可谓获益良多。五 附录部分代码1.DTFT函数代码：function X=dtft(x,w)%x为输入离散序列X = x*exp(-j*1:length(x)'*w

15、); %w为频率数组2.物理分辨率DFT/DTFT仿真clear allclfj=sqrt(-1);f1=2; f2=2.5; f3=3;fs=10;w1=2*pi*f1/fs; w2=2*pi*f2/fs; w3=2*pi*f3/fs;L1=10; L=800;N1=16; N2=64; N3=256; N=800;n=0:L1-1;x=sin(w1*n)+sin(w2*n)+sin(w3*n);vct_w1=0: 2*pi/N1: 2*pi-2*pi/N1;vct_w2=0: 2*pi/N2: 2*pi-2*pi/N2;vct_w3=0: 2*pi/N3: 2*pi-2*pi/N3;vct

16、_w =0: 2*pi/N: 2*pi-2*pi/N;X1=dtft(x,vct_w1); X2=dtft(x,vct_w2);X3=dtft(x,vct_w3);X =dtft(x,vct_w);subplot(3,1,1)axis(0 fs/2 0 12)hold onbar (0:N1-1)/N1*fs,abs(X1),0) plot(0:N1-1)/N1*fs,abs(X1),'b*') plot(fs*vct_w/pi/2,abs(X),'b:')legend('DFT, L=10 N=16','DTFT',0)subp

17、lot(3,1,2)axis(0 fs/2 0 12)hold onbar (0:N2-1)/N2*fs,abs(X2),0) plot(0:N2-1)/N2*fs,abs(X2),'b*') plot(fs*vct_w/pi/2,abs(X),'b:')legend('DFT, L=10 N=64','DTFT',0)subplot(3,1,3)axis(0 fs/2 0 12)hold onbar (0:N3-1)/N3*fs,abs(X3),0) plot(0:N3-1)/N3*fs,abs(X3),'b*'

18、) plot(fs*vct_w/pi/2,abs(X),'b:')legend('DFT, L=10 N=256','DTFT',0)3.计算分辨率DFT/DTFT仿真clfL2=100; L=800;N1=16;N2=64;N3=256;n=0:L2-1;x=cos(w1*n)+cos(w2*n)+cos(w3*n);vct_w1=0: 2*pi/N1: pi;vct_w2=0: 2*pi/N2: pi;vct_w3=0: 2*pi/N3: pi;vct_w =0: 2*pi/N: pi;X1=dtft(x,vct_w1);X2=dtft(x,vct_w2);X3=dtft(x,vct_w3);X =dtft(x,vct_w);subplot(3,1,1)axis(0 fs/2 0 60)hold onNN HN1=size(abs(X1);bar (0:HN1-1)/N1*fs,abs(X1),