用FFT做谱分析

实验二用FFT做谱分析实验项目名称：用FFT做谱分析实验项目性质：综合性实验所属课程名称：数字信号处理实验计划学时：4一.实验目的(1)进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。(2)熟悉FFT算法原理和FFT子程序的应用。(3)学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。二.实验内容和要求(1)对下列典型信号进行谱分析：ttttxnnxnnxnnnnnnxnnnnnnxnRnx20cos16cos8cos)(8sin)(4cos)(7430,0,3,4)(7430,0,8,1)()()(6543241其他其他，cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t)这里给出针对各信号的FFT变换区间N以及对连续信号tx6的采样频率sf，供实验时参考：643216,64:168654321，，，：，，，，NHzftxNnxnxnxnxnxs(2)令)()()(54nxnxnx，用FFT计算8点和16点离散傅立叶变换，)]([)(nxDFTkX并根据DFT的对称性，由)(kX求出)]([)(44nxDFTkX和)]([)(55nxDFTkX（1）中所得结果进行比较。[提示：取16N时，nNxnx44，nNxnx55。](3)令)()()(54njxnxnx，重复（2）。三.实验主要仪器设备和材料计算机，MATLAB6.5或以上版本四.实验方法、步骤及结果测试(1)复习DFT的定义、性质和用DFT做谱分析的有关内容。(2)复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。(3)编制信号产生子程序，产生下列典型信号供谱分析用。1423456745845()()1,03()8,470,4,03()3,470,()cos4()sin8()cos8cos16cos20()()()()()()xnRnnnxnnnnnnxnnnnxnnxnnxttttxnxnxnxnxnjxn其他其他应当注意，如果给出的是连续信号txa，则首先要根据其最高频率确定采样速率sf以及由频率分辨率选择采样点数N，然后对其进行软件采样（即计算nTxnxa，10Nn），产生对应序列nx。对信号tx6，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。请实验者根据DFT的隐含周期性思考这个问题。(4)编写主程序。图2-1给出了主程序框图，供参考。本实验的主程序比较简单，直接根据上图给出的框图编写主程序即可，编程中的难点是FFT子程序。不过，各种语言的FFT子程序都可以在有关的信号处理的程序库中找到，由于C语言当前最普及，所以为了为实验提供方便，下面给出C语言FFT函数，供参考/*DIT-FFT函数（C语言）*/fft—基2DIT—FFT函数要求：指向复数数组指针X，FFT长度为2m，m为正整数FFT输出结果放在输入复数数组中。开始读入长度N调用信号产生子程序产生实验信号调用绘图子程序（函数）绘制时间序列波形图调用FFT子程序（函数）计算信号的DFT调用绘图子程序（函数）绘制KX曲线结束图2-1主程序框图/*计算N点FFT子程序*//*xr：=信号序列实部，xi：=信号序列虚部，N：=FFT变换区间长度N=2^M*//*如果信号长度小于N，应该给xr，xi后面补0*//*计算如果X（K）的实部和虚部分别储存在数组xr和xi中*/VoidFft(doublexr[],doublexi[],intN,intM){intL,B,J,P,k,i;doublerPartKB,iPartKB;doublerCf[128],iCf[128]/*计算旋转因子*/doublePI2=8.0*atan(1.0);for(i=0;iN;i++){rCf[i]=cos(i*PI2/N);iCf[i]=sin(i*PI2/N);}ChangeOrder(xr,xi,N);/*计算各级蝶形*/for(L=1;L=M;L++){B=(int)(pow(2,(L-1))+0.5);for(J=0;J=B-1;J++){P=J*((int)(pow(2,(M-L))+0.5));for(k=J;k=N-1;k+=(int)(pow(2,L)+0.5)){rPartKB=xr[k+B]*rCf[P]-xi[k+B]*iCf[P];iPartKB=xi[k+B]*rCf[P]+xr[k+B]*iCf[P]xr[k+B]=xr[k]-rPartKB;xi[k+B]=xi[k]-iPartKB;xr[k]=xr[k]+rPartKB;xi[k]=xi[k]+iPartKB;}}}}/*倒序子程序*/voidChangeOrdor(doublexr[],doublexi[],intN){intLH,N1,I,J,K;doubleT;LH=N/2;J=LH;N1=N–2;for(I=1;I=N1;I++){if(IJ){T=xr[I];xr[I]=xr[J];xr[J]=T;T=xi[I];xi[I]=xi[J];xi[J]=T;}K=LH;while(J=K){J=J-K;K=(int)(K/2+0.5);}J=J+K;}}(5)按实验内容要求，上机实验并写出实验报告。本实验采用的是MATLAB语言，因此FFT子程序直接调用MATLAB语言中的FFT函数就可以实现。下面给出完整的MATLAB程序%实验二，用FFT做谱分析b=menu('请选择信号x1(n)--x8(n)','x1(n)','x2(n)','x3(n)','x4(n)','x5(n)','x6(n)','x7=x4+x5','x8=x4+jx5','Exit');ifb==9b=0;endi=0;closeall;while(b)ifb==6temp=menu('请选择FFT变换区间长度N','N=16','N=32','N=64');iftemp==1N=16;elseiftemp==2N=32;elseN=64;endfs=64;n=0:N-1;x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);elsetemp=menu('请选择FFT变换区间长度N','N=8','N=16','N=32');iftemp==1N=8;elseiftemp==2N=16;elseN=32;endifb==1x=[11110000];elseifb==2x=[12344321];elseifb==3x=[43211234];elseifb==4n=0:N-1;x=cos(0.25*pi*n);elseifb==5n=0:N-1;x=sin((pi*n)/8);elseifb==7n=0:N-1;x=cos(n*pi/4)+sin(n*pi/8);elseifb==8n=0:N-1;x=cos(n*pi/4)+j*sin(n*pi/8);endendendendendendendend%%TOCalculateFFTf=fft(x,N);i=i+1;figure(i);printf(x,abs(f),abs(N),abs(b));ifN==16ifb==7k=conj(f);x4=(f+k)/2;%Re[X7(k)=x4(k)figure(i+2);subplot(2,2,1);stem(abs(x4),'.');xlabel('k');ylabel('|X4(k)|');title('恢复后的X4(k)');x5=(f-k)/2;%jIm[X7(k)=X5(k)subplot(2,2,3);Stem(abs(x5),'.');xlabel('k');ylabel('|X5(k)|');title('恢复后的X5(k)');endifb==8k(1)=conj(f(1));form=2:Nk(m)=conj(f(N-m+2));endfe=(x+k)/2;%求X8(k)的共轭对称分量fo=(x-k)/2;%求X8(k)的共轭反对称分量xr=ifft(fe,N);%xr=x4(n)b=4;figure(i+1)printf(xr,abs(fe),abs(N),abs(b));xi=ifft(fo,N)/j;%xi=x5(n)b=5;figure(i+2)printf(xi,abs(f),abs(N),abs(b));endendb=menu('请选择信号x1(n)--x8(n)','x1(n)','x2(n)','x3(n)','x4(n)','x5(n)','x6(n)','x7=x4+x5','x8=x4+jx5','Exit');ifb==9b=0;endcloseall;end程序运行结果：按照实验内容选择~x6(n)以及x7(n)=x4(n)+x5(n)、x8(n)=x4(n)+jx5(n)进行谱分析，输出x1(n)~x5(n)的波形以及8点DFT和16点DFT如图2-2到图2-11所示，输出x6(n)的16点、32点和64点采样序列及其DFT如图2-12、13、14所示。x1(n)8点：0123456700.51时域信号0123456701234频域信号16点：05101500.51时域信号05101501234频域信号X2(n)8点0123456701234时域信号0123456705101520频域信号16点05101501234时域信号05101505101520频域信号X3(n)8点0123456701234时域信号0123456705101520频域信号16点05101501234时域信号05101505101520频域信号X4(n)8点01234567-1-0.500.51时域信号0123456701234频域信号16点051015-1-0.500.51时域信号05101502468频域信号X5(n)8点0123456700.51时域信号012345670246频域信号16点051015-1-0.500.51时域信号05101502468频域信号X6(n)8点16点选7时，首先计算并图示x7(n)=[x4(n)+x5(n)]*R8(n)和x7(n)=[x4(n)+x5(n)]*R16(n)及其DFT，如图2-15、2-16所示；然后程序自动计算并绘图验证DFT的共轭对称性。当N=16时，x4(n)=x4(N-n)，x5(n)=-x5(N-n)。即x4(n)为x7(n)的共轭对称分量，而x5(n)是x7(n)的共轭反对称分量。根据DFT的共轭对称性，应有如下结果：X7(k)=DFT[x7(n)]16点=Re[X7(k)]+jIm[X7(k)]X4(k)=DFT[x4(n)]16点=Re[X7(k)]X5(K)=DFT[x5(n)]16点=jIm[X7(k)]图2-17绘出了Re[X7(k)和jIm[X7(k)]的模，而他们正是图2-9和图2-11中的16点|X4(k)|和|X5(k)|。图2-15x7(n)和8点DFT图2-16x7(n)和16点DFT选8时，计算并图示X8(n)=[x4(n)+jx5(n)]R8(n)和x8(n)=[x4(n)+jx5(n)]*R16(n)及其DFT，如图2-18、2-19所示。然后程序自动计算并绘图验证DFT的共轭对称性的第二种形式：如果x(n)=xr(n)+jxi(n)，X(k)=DFT[x(n)]=Xep(k)+Xop(k)，则Xep(k)=DFT[xr(n)]，Xop(k