利用DFT进行频谱分析

利用DFT进行频谱分析内容与要求利用DFT对多种信号（例如由多个正弦信号组成的信号）进行频谱分析，并研究不同数据长度、补零、加窗等对频率分辨率的影响。方法原理1、引入当数字计算机对信号进行频谱分析时，要求信号必须以离散值作为输入，而计算机输出所得的频谱值自然也是离散的。因此，要使信号是时间的连续函数、频谱是频率的连续函数或者信号及频谱二者都是变量的连续函数这三种形式的信号能用数字计算机进行计算，必须针对每一种形式的具体情况，或者在时域与频域上取样，或者在时域上取样，或者在频域上取样。信号在时域上取样导致频率的周期函数，在频域上取样导致时域的周期函数，最后都将使原时间函数和频率函数二者都成为周期离散的函数。我们采用DFT（离散傅里叶变换）来对连续时间信号的傅里叶变换进行逼近，进而分析连续时间信号的频谱。离散傅里叶变换是有限长序列的傅里叶变换，它相当于把信号的傅里叶变换进行等频率间隔采样，并且有限长序列的离散傅里叶变换和周期序列的离散傅里叶级数本质是一样的。2、推导离散傅里叶级数定义为nkjNkppekxNnxN210)(1)(π将上式两端乘以nmjNeπ2并对n在0~N-1求和可得10)(1101010)(10N2N2N2)()(1)(NnmknjNNkpNnNkmknjpNnnmjpekXekXNenxπππ因为mk1mk0)(N)(10)(N2N2N2-1-1N11mkjmkjNnmknjeeeNπππ所以1010)()()(N2NkpNnnmjpmkkXenxπ这样10N2)()(NnnmjppenxmXπ用k代替m得10N2)()(NnnkjpPenxkXπ令N2πjNeW，则DFS10)()()(NnnkNpppWnxkXnxIDFS10)(1)()(NnnkNpppWkXNnxkX其中)()(kXnxpp、都是周期为N的周期序列，DFS[·]表示离散傅里叶级数正变换，IDFS[·]表示离散傅里叶级数反变换。习惯上，对于长为N的周期序列，把0nN-1区间称为主值区，把)1(~)0(Nxxpp称为)(nxp的主值序列，同样也称)1(~)0(NXXpp为)(kXp的主值序列。由于)()()(nRnxnxNp，对于周期序列)(nxp仅有N个独立样值，对于任何一个周期进行研究就可以得到它的全部信息。在主值区研究)(nxp与)(nx是等价的，因此在主值区计算DFS和DFT是相等的，所以DFT计算公式形式与DFS基本相同。其关系为)()()(nRnxnxNp，)()()(kRkXkXNp所以离散傅里叶正变换WnkNNnnxnxDFTkX100kN-13、定义DFT：设有限长序列x(n)长为N（0nN-1），其离散傅里叶变换是一个长为N的频率有限长序列（0kN-1），其正变换为WnkNNnnxnxDFTkX100kN-1DFT的分辨率：指其能够分辨的最小频率间隔。频率分辨率主要由数据截断的长度决定，即时间长度的倒数。也可以说由时间窗函数的傅里叶变换，即谱窗的主瓣宽度决定。不同的谱窗的主瓣宽度不同。矩形窗的主瓣宽度最窄，但其副瓣最高（不利于对频率相邻弱信号的分辨），其它常用的窗函数的主瓣宽度与其副瓣高度近似存在反比关系。主瓣窄，副瓣高，有利于相邻强信号的分辨，但不利于相邻弱信号的分辨。主瓣宽，副瓣低不利于相邻强信号的分辨，但可能有利于相邻弱信号的分辨。4、实质把有限长序列当做周期序列的主值序列进行DFS变换，x(n)、X(k)的长度均为N，都是N个独立值，因此二者具有的信息量是相等的。已知x(n)可以唯一确定X(k),已知X(k)可以唯一确定x(n)。虽然离散傅里叶变换是两个有限长序列之间的变化，但它们是利用DFS关系推导出来的，因而隐含着周期性作业内容1、几种信号的频谱分析（1）自定义DFT函数functionxk=dft_1(xn)N=length(xn);WN=exp(-1i*2*pi/N);n=0:1:N-1;k=0:1:N-1;nk=k'*n;WNnk=WN.^(nk);xk=xn*WNnk;end（2）对信号进行频谱分析N=input('N=');n=0:1:N-1;xn=input('xn=');Xk=dft_1(xn1,N);subplot(3,1,1)stem(n,xn,'.k');xlabel('n');axis([0,N,-2.5,2.5]);w=2*pi*(0:1:2047)/2048;Xw=xn*exp(-1i*n'*w);subplot(3,1,2);plot(w/pi,abs(Xw));xlabel('w');axis([0,1,0,N]);subplot(3,1,3)k1=0:1:N-1;w1=2*pi/N*k1;stem(w1/pi,abs(Xk),'.k');xlabel('w');axis([0,1,0,N]);（3）几种不同信号的频谱分析图（N=100,wn=boxcar(N)矩形窗）①xn=cos(0.4*pi*n)+sin(0.6*pi*n)②xn=0.02*n③xn=heaviside(n)2、几种因素对频率分辨率的影响（1）数据长度（取xn=cos(0.4*pi*n)+sin(0.6*pi*n)，wn=boxcar(N)矩形窗）①N=10②N=20③N=100结论：由图可见，数据长度的增长改变了频谱混叠作用，提高了物理分辨率。（2）补零（取xn=cos(0.4*pi*n)+sin(0.6*pi*n),N=100,wn=boxcar(N)矩形窗）①不补零②补零至N=300结论：由图可见，补零只改变了Xk的密度，截断函数的频谱混叠作用没有改变。这说明，补零仅仅是提高了计算分辨率，得到的是高密度频谱，而得不到高分辨率谱。（3）加窗（取xn=cos(0.*4*pi*n)+sin(0.6*pi*n),N=100）①矩形窗wn=boxcar(N)i）窗长度40ii）窗长度100②三角形窗wn=triang(N)i）窗长度40ii）窗长度100③汉宁窗wn=hanning(N)i）窗长度40ii）窗长度100④海明窗wn=hamming(N)i）窗长度40ii）窗长度100⑤布拉克曼窗wn=blackman(N)i）窗长度40ii）窗长度100结论：不同的谱窗的主瓣宽度不同。矩形窗的主瓣宽度最窄。窗长度为数据长度时分辨率最高。3.与MATLAB自带函数对比（1）幅频特性①代码N=100;n=0:1:N-1;xn=cos(0.4*pi*n)+sin(0.6*pi*n);subplot(3,1,1)stem(n,xn,'.k');title('时域序列图xn');xlabel('n');axis([0,N,-2.5,2.5]);Xw=fft(xn,N);subplot(3,1,2)k1=0:1:N-1;w1=2*pi/N*k1;stem(w1/pi,abs(Xw),'.k');title('fft幅频特性');xlabel('频率');axis([0,1,0,N]);Xk=dft_1(xn,N);subplot(3,1,3)k2=0:1:N-1;w2=2*pi/N*k2;stem(w2/pi,abs(Xk),'.k');title('dft幅频特性');xlabel('频率');axis([0,1,0,N]);②结果（2）相频特性①代码N=100;n=0:1:N-1;xn=cos(0.4*pi*n)+sin(0.6*pi*n);subplot(3,1,1)stem(n,xn,'.k');title('时域序列图xn');xlabel('n');axis([0,N,-2.5,2.5]);Xw=fft(xn,N);subplot(3,1,2)k2=0:1:N-1;w2=2*pi/N*k2;stem(w2/pi,angle(Xw),'.k');title('fft相频特性');xlabel('频率');axis([0,1,0,N]);Xk=dft_1(xn,N);subplot(3,1,3)k1=0:1:N-1;w1=2*pi/N*k1;stem(w1/pi,angle(Xk),'.k');title('dft相频特性');xlabel('频率');axis([0,1,0,N]);②结果结论：dft与MATLAB自带函数相比，相频特性不同，幅频特性无显著区别。