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

《数字信号处理》课程研究性学习报告姓名学号同组成员指导教师时间DFT近似计算信号频谱专题研讨【目的】(1)掌握利用DFT近似计算不同类型信号频谱的原理和方法。(2)理解误差产生的原因及减小误差的方法。(3)培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】基本题1．利用DFT分析x(t)=Acos(2f1t)+Bcos(2f2t)的频谱，其中f1=100Hz，f2=120Hz。(1)A=B=1;(2)A=1,B=0.2。要求选择不同的窗函数。【题目分析】由DFT结果可得：通过对不同抽样频率，不同的窗函数对信号的DFT结果可以看出：在对信号做DFT时，由于对信号进行截短，因此会产生频谱泄漏，要想从频谱中很好的分辨出个频率分量，需要考虑时域抽样频率，所加的窗函数，窗函数的长度，以及DFT的点数等参数对结果的影响：因为f1f2,所以为满足抽样定理，应使fsam=2f2即fsam=240Hz△f=f2-f1=20N=cfsam/△f（1）A=B=1x(t)=cos(2f1t)+cos(2f2t)矩形窗1：条件：fsam=240HzN=20L=512矩形窗2：条件：fsam=600HzN=40L=512矩形窗3：fsam=1200HzN=80L=512Hamming窗1：N=40;L=512;fs=600;Hamming窗2：N=120;L=512;fs=600;(2)A=1,B=0.2x(t)=cos(2f1t)+0.2cos(2f2t)矩形窗：N=100;L=512;fs=600;hamming窗：N=100;L=512;fs=600;【仿真结果】（1）A=B=1x(t)=cos(2f1t)+cos(2f2t)矩形窗：Hamming窗：（2）A=1,B=0.2x(t)=cos(2f1t)+0.2cos(2f2t)矩形窗：【结果分析】在（1）中进行矩形窗仿真时，我们选择了不同的fsam，分别为240,600,1200它们均满足抽样定理，但是我们仍旧发现，在240hz时出现了混叠现象。所以，在实际应用中抽样频率应大于最低抽样频率3-5倍才能有更好的结果。进行hamming窗仿真时，在保证抽样频率相同的条件下，取不同的长度也40,120。其中，N=40不满足N=60的要求，我们可以看到出现了混叠，而N=120时，仿真效果良好。在（2）的条件下进行仿真时，我们选取了相同的N、L、fsam值，但是分别使用了矩形窗和hamming窗。使用矩形窗时，幅度较小的峰值与旁瓣的幅度接近，甚至难以区分，效果不理想。使用hamming窗后，泄露现象被有效遏制，所以可以清楚区分主瓣、旁瓣。综上所述，在选择参数进行DFT变换时，应该保证抽样频率满足抽样定理，并且能大于最小抽样值3-5倍。长度选择保证N=cfsam/△f。为防止泄露现象，特别是峰值之间差异较大时，应该选择加特殊的窗，如hamming窗。【自主学习内容】【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：按照理论分析最小抽样频率只需要满足2fmax就可以满足抽样定理，但在仿真中发现该频率无法满足要求，频谱发生严重的混叠。所以抽样频率应为最小抽样频率3-5倍。另外，在使用哈明窗作为窗函数时，按照理论分析，哈明窗长度满足N2fs/f，结合题中条件，当fs=600时N为60时恰好可以分出频谱，而实际中N=60时无法分出两个频率分量，当N=90时则可以分出。因此在做DFT是窗函数长度应大于最小长度。【问题探究】【仿真程序】（1）矩形窗1程序：N=20;L=512;f1=100;f2=120;fs=240;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);X=fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X));ylabel('矩形窗1')矩形窗2程序：N=40;L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);X=fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X));ylabel('矩形窗2')hamming1:N=40;L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);wh=(hamming(N))';x=x.*wh;X=fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X));ylabel('hamming1')hamming2N=120;L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);wh=(hamming(N))';x=x.*wh;X=fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X));ylabel('hamming1')（2）矩形窗：N=100;L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);X=fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X));ylabel('(2)矩形窗')hamming窗：N=100;L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);wh=(hamming(N))';x=x.*wh;X=fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(X));ylabel('(2)hamming´°')【研讨题目】基本题2．已知一离散序列为31,,1,0),π2.0sin(][kkkx(1)用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；(2)对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱峰的频率；(3)讨论所获得的结果，给出你的结论。该结论对序列的频谱计算有何指导意义？【题目分析】对离散非周期序列做DFT时，所得结果相当于对序列做傅里叶变换，再在一个周期内进行频域抽样，而DFT的点数为频域抽样的点数【温磬提示】在计算离散非周期序列频谱时常用/作为横坐标，称/为归一化频率normalizedfrequency)。在画频谱时需给出横坐标。每幅图下都需给出简要的文字说明。由于离散非周期序列频谱是周期的，所以在计算时不需要用fftshift函数对fft计算的结果进行重新排列。【序列频谱计算的基本方法】在仿真软件中fft函数可以进行离散序列的DFT其调用形式为：fft(x,N)计算序列的N点DFT，若序列长度为M，当MN时则将原序列截短为N点序列再进行计算，若NM则将序列补零后在做N点的DFT.【仿真结果】【结果分析】对序列进行32点的DFT得到的频谱如图所示，谱峰的频率分别为0.09pi和0.91pi（-0.09pi）对序列补零后再做DFT相当于增加了DFT的点数，即频域抽样的点数，而原离散序列没有改变，其傅里叶变换结果也没有改变，同时若DFT点数太少则获得的频谱信息过少，有可能会使得重要的频率信息丢失当1，仿真结果为：L=32时谱峰的频率为0.2pi和1.8pi由结果可知，DFT点数越多，产生的离散谱中含有的信息也就越多，得到的频谱能更好的反应原连续谱中的信息在对离散序列用DFT做谱分析时，应当适当增加DFT的点数，以减小栅栏效应【自主学习内容】【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：【问题探究】【仿真程序】k=0:31x=sin(0.2*pi*k)k1=0:31x_32=fft(x,32)subplot(5,1,1)plot(2*k1/32,abs(x_32),'g')title('L=32')k2=0:63x_64=fft(x,64)subplot(5,1,2)plot(2*k2/64,abs(x_64),'r')title('L=64')k3=0:127x_128=fft(x,128)subplot(5,1,3)plot(2*k3/128,abs(x_128),'m')title('L=128')k4=0:255x_256=fft(x,256)subplot(5,1,4)plot(2*k4/256,abs(x_256),'y')title('L=256')k5=0:511x_512=fft(x,512)subplot(5,1,5)plot(2*k5/512,abs(x_512),'k')title('L=512')figurek=0:31x=sin(0.2*pi*k)L=0:1023X=fft(x,1024)subplot(2,1,1)plot(L/1024,abs(X),'r')holdon;k1=0:31x_32=fft(x,32)stem(2*k1/32,abs(x_32),'x','b')title('L=32')L=512;f1=0.1;fs=0.4;T=1/fs;ws=2*pi*fs;k=(0:1:L-1)*T;x=sin(0.2*pi*k);X=fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w/pi,abs(X));figurek=0:31x=sin(0.2*pi*k)L=0:1023X=fft(x,1024)subplot(2,1,1)plot(2*L/1024,abs(X),'r')holdon;k1=0:31x_32=fft(x,32)stem(2*k1/32,abs(x_32),'x','b')title('L=32')k=0:31x=sin(0.2*pi*k)L=0:1023X=fft(x,1024)subplot(2,1,2)plot(2*L/1024,abs(X),'r')holdonk3=0:127x_128=fft(x,128)stem(2*k3/128,abs(x_128),'*','b')title('L=1024')【研讨题目】基本题3．已知一离散序列为x[k]=Acos0k+Bcos0+)k)。用长度N=64的哈明窗对信号截短后近似计算其频谱。试用不同的A和B的取值，确定用哈明窗能分辩的最小的谱峰间隔Ncπ2Δw中c的值。（M2-3）【题目分析】Hamming窗函数的幅值有中心向两端逐渐减弱，因而其高频分量明显减小，频谱中旁瓣的幅度较小，主瓣峰值与第一个旁瓣峰值相对衰减很大，hamming窗以增加主瓣宽度来降低旁瓣能量，用hamming窗极端频谱时要求能分辨的谱峰的间隔Δƒ≥c/Tp=c*fs/N。为使两个峰值间有明显差别，取A=5,B=1。【仿真结果】【结果分析】由于信号中存在一个较弱的频率分量，可以用hamming窗来减小旁瓣引起的频率泄露，但当c=2时并不能分辨出两个谱峰，由图形可以看出能分辨出图形的c的最小值为c=2.5，随着A与B比例的增加，即两个幅值分量的差距增大，频谱更加难以分辨，频率泄露越大，失真越大。【自主学习内容】【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：【问题探究】在离散序列频谱计算中为何要用窗函数？用不同的窗函数对计算结果有何影响？与矩形窗相比哈明窗有何特点?如何选择窗函数？【仿真程序】cleark=0:63;A