1、FFT在功率谱密度计算中地应用FFT在功率谱密度计算中的应用一、FFT算法理论依据和编程思想FFT算法的基本思想： 考察DFT与IDFT的运算发现，利用以下两个特性可减少运算量： )系数 是一个周期函数，它的周期性和对称性可利用来改进运算，提高计算效率。 如: 因此 利用这些周期性和对称性，DFT运算中有些项可合并； ）利用WNnk的周期性和对称性，把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。因为DFT的计算量正比于N2，N小计算量也就小。 FFT算法正是基于这样的基本思想发展起来的。它有多种形式，下面是按时间抽取的FFT（N点DFT运算的分解） 先从一个特殊情况开始，假定N是
2、2的整数次方，N=2M，M：正整数 1.将N点的DFT分解为两个N/2点的DFT： 首先将序列x（n）分解为两组，一组为偶数项，一组为奇数项 r=0,1，N/2-1 将DFT运算也相应分为两组： 其中X1（k）和X2（k）分别是x1（r）和x2（r）的N/2点DFT。 可见，一个N点的DFT可以分解为两个N/2点的DFT，这两个N/2点的DFT再按照上面（1）式合成为一个N点DFT，注意到，X1（k），X2（k）有N/2个点，即k=0，1，N/2-1，由（1）式得到X（k）只有N/2点，而实际上X（k）有N个点，即k=0，1，N-1，要用X1（k），X2（k）表示全部X（K）值，还必须应用系数
3、w的周期性和对称性。2.X(k)的(N/2)N-1点表示： 由X(k)= X1（k）+wkN X2（k）, k=0,1,2,N/2-1 得： ，（2a）因为 , 且 同样 。 考虑到WNk对称性： 。 故 （2b） （2a）式表示了X（k）前半部分k=0N/2-1时的组成方式，（2b）式则表示了后半部分k=N/2N-1时的组成方式。这两式所表示的运算过程可用一个称作蝶形的信号流图来表示。 3.蝶形信号流图： 如图1（a）所示，图中左面两支为输入，中间以一个小圆圈表示加、减运算，右上支为相加输出，右下支为相减输出，如果在某一支路上信号需要进行乘法运算，则在该支路上标以箭头，并将相乘的系数标在简头
4、边，这样（2a），（2b）所表示的运算，可用图1（b）所表示的“蝶形结”来表示。采用这种表示法，可将以上以讨论的分解过程用计算流图来表示。 图2.6所示为N=23=8的例子。通过这样分解以后，每一个N/2点DFT只需要 图2.6 N点DFT分解为2个N/2点DFT（N=8）(N/2)2= N2/4次复数乘法运算，两个N/2点的DFT需要2(N/2)2= N2/2 次复乘，再加上将两个N/2点DFT合成为N点DFT时，在蝶形结前的N/2次复乘，共需要(N/2)2+N/2 N2/2次复乘，由此可见，经过这样的分解处理，运算量差不多节省了一倍。 4.将N/2点的DFT分解为两个N/4点的DFT： 既
5、然这样分解是有效的，由于N=2M，N/2仍然是偶数，因此可对两个N/2点的DFT再分别作进一步分解，例如对x1（r）和x2（r）可以再按其偶数部分及奇数部分分解为两个N/4点的DFT，既然这样分解是有效的，由于N=2M，N/2仍然是偶数，因此可对两个N/2点的DFT再分别作进一步分解，例如对x1（r）和x2（r）可以再按其偶数部分及奇数部分分解为两个N/4点的DFT， l=0,1,，N/4-1 而 同样X2(k)也可这样分解，并且将系数统一为 ，这样一个8点DFT就可分解为四个2点的DFT，如图2.7所示。 图2.7 N点DFT分解为4个N/4点的DFT(N=8)5.2个点DFT的表示： 最后
6、剩下的是2点DFT，它可以用一个蝶形结表示，例如，x（0），x（4）所组成的2点DFT就可表示式： 这样，一个8点的完整的按时间抽取运算的流图如图2.8所示。 由于这样的方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的，所以称为“按时间抽取法”或“时间抽取法”。 6.时间抽取法FFT运算特点： （1）蝶形运算 对任何一2的整数幂N=2M，总是可以通过M次分解最后完全成为2点的DFT运算。这样的M次分解，就构成从x(n)到X(k)的M级运算过程。从上面的流图可看到，每一级运算都由N/2个蝶形运算构成。因此每一级运算都需要 (N/2次复乘和N次复加(每个结作加、减各一次)，这样，经过时间
7、抽取后M级运算总共需要的运算： 复乘 复加 NM=Nlog2N 实际运算量与这个数字稍有出入，因为W 这几个系数实际上都不用乘法运算，因此在上面N=8的例子中，实际上只有两个系数WN1及WN3是需要乘法运算的。 用时间抽取法所需的计算量，不论是复乘还是复加都与Nlog2N成正比，而直接运算时则与N2成正比。例N=2048，N2=4194304，(N/2)log2N=11264，N2/(N/2)log2N=392.4倍。FFT显然要比直接法快得多。(2）原位计算 当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位计算。 例如，N=8的F
8、FT运算，输入x(0)，x(4)，x(2)，x(6)，x(7)可分别存入A(1)，A(2)，A(8)这9个存储单元中，在第一级运算中，首先是存储单元A(1)，A(2)中x(0)，x(4)进入蝶形运算，x(0)，x(4)输入运算器后，其数值不再需要保存，因此蝶形运算的结果可仍然送回存储单元A(1)，A(2)中保存，然后A(3)，A(4)中x(2)，x(6)再进入蝶形运算，其结果再送回A(3)，A(4)，一直到算完A(7)，A(8)，则完成了第一级运算过程。第二级运算仍可采用这种原位的方式，但是进入蝶形结的组合关系不同，首先进入蝶形结的是A(1)、A(3)存储单元中的数据，运算结果仍可送回A(1)
9、、A(3)保存，然后进入蝶形结的是A(2)、A(4)，依此类推，每一级运算均可在原位进行，这种原位运算结构可节省存储单元，降低设备成本，还可节省找地址的时间。 （3）序数重排 对按时间抽取FFT的原位运算结构，当运算完毕时，这种结构存储单元A(1)、A(2)，A(8)中正好顺序存放着X(0)，X(1)，X(2)，X(7)，因此可直接按顺序输出，但这种原位运算的输入x(n)却不能按这种自然顺序存入存储单元中，而是按X(0)，X(4)，X(2)，X(6)，X(7)的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。例如，原来的自然顺
10、序应是x(1)的地方，现在放着x(4)，用二进制码表示这一规律时，则是在x(0 0 1)处放着x(1 0 0)，x(0 1 1)处放着x(1 1 0)。即将自然顺序的二进制码位倒置过来，第一位码变成最末位码，这样倒置以后的顺序正是输入所需要的顺序。下表列出N=8时按码位倒置规律所得的顺序，其结果与按时间抽取FFT流图中的输入顺序是一致的。 表一 码位倒置顺序 自然顺序二进码表示码位倒置码位倒置顺序0000000010011004201001023011110641000011510110156110011371111117 在实际运算中，一般直接将输入数据x（n）按码位倒置的顺序排好输入很不方
11、便，总是先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储，然后进行FFT的原位计算。目前有许多通用DSP芯片支持这种码位倒置的寻址功能。（4）蝶形类型随迭代次数成倍增加 观察8点FFT的三次迭代运算 第一级迭代，只有一种类型的蝶形运算系数W08 第二级迭代，有二种类型的蝶形运算系数W08、W28，参加运算的两个数据点间隔为2。 第三级迭代，有四类蝶形运算系数W08、W18、W28、W38，参加运算的两个数据点间隔为4。 所以，每次迭代的蝶形类型比上一次蝶代增加一倍，数据点间隔也增大一倍。7.功率谱密度的计算根据相关定理与维纳-辛钦关系式可得随机信号序列x(n
12、)的功率谱密度 (1)其估计值 (2) 如果观察到序列x(n)的N个值，即x(0), x(1),x(N-1)，就可以通过FFT直接求得X(k),再按式(2)求得Sx(k),其计算过程如图2.9所示。XI(k)二、程序设计本设计在主程序中分别调用了输入get_in()、倒序re_order()和蝶式运算butterfly(),功率谱密度计算power(),绘图on_draw五个子函数。1.程序框图（1） 运算主程序框图(2)、整序程序流程图 图7 整序程序流程图（3）蝶形运算框图 图8 蝶形运算框图三、程序源代码主程序/ FFT.cpp : Defines the class behaviors
13、 for the application.#include stdafx.h#include FFT.h#include FFTDlg.h/自己添加的.#includemath.h#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE = _FILE_;#endif/ CFftDrawingAppBEGIN_MESSAGE_MAP(CFftDrawingApp, CWinApp) /AFX_MSG_MAP(CFftDrawingApp) / NOTE - the ClassWizard will add a
14、nd remove mapping macros here. / DO NOT EDIT what you see in these blocks of generated code! /AFX_MSG ON_COMMAND(ID_HELP, CWinApp:OnHelp)END_MESSAGE_MAP()/ CFftDrawingApp constructionCFftDrawingApp:CFftDrawingApp() / TODO: add construction code here, / Place all significant initialization in InitIns
15、tance/ The one and only CFftDrawingApp objectCFftDrawingApp theApp;/ CFftDrawingApp initializationBOOL CFftDrawingApp:InitInstance() AfxEnableControlContainer(); / Standard initialization / If you are not using these features and wish to reduce the size / of your final executable, you should remove
16、from the following / the specific initialization routines you do not need.#ifdef _AFXDLL Enable3dControls(); / Call this when using MFC in a shared DLL#else Enable3dControlsStatic(); / Call this when linking to MFC statically#endif CFftDrawingDlg dlg; m_pMainWnd = &dlg; int nResponse = dlg.DoModal()
17、; if (nResponse = IDOK) / TODO: Place code here to handle when the dialog is / dismissed with OK else if (nResponse = IDCANCEL) / TODO: Place code here to handle when the dialog is / dismissed with Cancel / Since the dialog has been closed, return FALSE so that we exit the / application, rather than
18、 start the applications message pump. return FALSE;/*/自己添加的.void CFftDrawingApp:Fft(double *shibu, double *xubu) double PI=3.1415926535; int M,N; double x512,xr512,xi512,t512,jiaodu; M=9; N=1M; jiaodu=(-2*PI)/N; /数据点输入: /caidian(); for (int i=0;iN;i+) if (i9) xi=1; else xi=0; /数据点输入结束. /将数据点重新排序,以实现
19、fft计算. /paixu(); int kuaishu=1,geshu; for( i=0;iN;i+) ti=xi; for(int pi=1;piM;pi+) kuaishu=1(pi); geshu=N/kuaishu; for(int fi=0;fikuaishu;fi+) if(fi%2=0) for(int fk=0;fkgeshu;fk+) xfi*geshu+fk=tfi*geshu+fk*2; if(fi%2=1) for(int fk=0;fkgeshu;fk+) xfi*geshu+fk=t(fi-1)*geshu+1+fk*2; for(int i=0;i512;i+
20、) ti=xi; /排序结束. /进行fft计算: /fft(); int zengzhi; double t1r,t2r,t1i,t2i; double jd; /计算的角度 int wz=1M; for( i=0;i=0;i1-) wz=wz/2; kuaishu=1(i1); geshu=N/kuaishu; zengzhi=kuaishu; for(int i2=0;i2kuaishu;i2+) for(int i3=0;i3(geshu/2);i3+) jd=jiaodu*zengzhi*i3; t1r=xrgeshu*i2+i3; t1i=xigeshu*i2+i3; t2r=xr
21、geshu*i2+i3+geshu/2; t2i=xigeshu*i2+i3+geshu/2; xrgeshu*i2+i3=t1r+t2r*cos(jd)-t2i*sin(jd); xigeshu*i2+i3=t1i+t2i*cos(jd)+t2r*sin(jd); xrgeshu*i2+i3+geshu/2=t1r-t2r*cos(jd)+t2i*sin(jd); xigeshu*i2+i3+geshu/2=t1i-t2i*cos(jd)-t2r*sin(jd); for(i=0;i512;i+) shibui=xri; xubui=xii; 四、计算实例采样函数：x(t)=A*cos(2*
22、PI*f1*t)+B*sin(2*PI*f2*t)采样点数：N=512计算结果：(打开程序运行时保存的文件fftdata.txt)操作过程：1. 打开生成的exe文件（FFT.exe），则会出现以下窗口：2. 点击软件界面的“参数导入”按钮，就可以计算完成了，生成FFTdata.Txt文件保存在所在文件夹。计算结果为：X0=17.4235+0 j X1=17.4211+0.606382 j X2=17.4135+1.21888 j X3=17.4004+1.84378 j X4=17.381+2.48781 j X5=17.3543+3.15828 j X6=17.3188+3.86342 j
23、 X7=17.2725+4.61272 j X8=17.2128+5.41733 j X9=17.1364+6.2907 j X10=17.0389+7.2493 j X11=16.9144+8.31383 j X12=16.7553+9.5107 j X13=16.5514+10.8745 j X14=16.2885+12.4513 j X15=15.9471+14.3045 j X16=15.4992+16.5237 j X17=14.9032+19.2398 j X18=14.0953+22.6524 j X19=12.9727+27.0821 j X20=11.3582+33.0782
24、 j X21=8.91951+41.6675 j X22=4.95384+55.0175 j X23=-2.32497+78.6307 j X24=-19.1725+131.765 j X25=-93.5606+362.07 j X26=206.277+-559.291 jX27=78.496+-164.482 j X28=57.6634+-98.7266 j X29=49.4477+-71.6923 jX30=45.2766+-56.9855 j X31=42.9308+-47.7638 j X32=41.5812+-41.462 j X33=40.8499+-36.9021 j X34=4
25、0.5424+-33.4674j X35=40.5503+-30.8044 jX36=40.8107+-28.6962 j X37=41.2871+-27.0029j X38=41.9592+-25.6305 jX39=42.8182+-24.514 j X40=43.8635+-23.6076 j X41=45.1019+-22.8788 jX42=46.5468+-22.3045 j X43=48.2183+-21.8689j X44=50.1444+-21.5618 jX45=52.362+-21.3781 j X46=54.9198+-21.3171 j X47=57.8815+-21
26、.383 j X48=61.3308+-21.5851 j X49=65.3799+-21.9394j X50=70.1808+-22.4703 jX51=75.9449+-23.2143 j X52=82.9737+-24.2255j X53=91.7124+-25.5861 j X54=102.846+-27.4244 j X55=117.487+-29.9503j X56=137.564+-33.5302 j X57=166.744+-38.8621 j X58=212.964+-47.4581j X59=297.164+-63.3121 j X60=498.461+-101.515 j
27、 X61=1615.21+-314.361j X62=-1256.78+233.677 j X63=-446.848+79.3409 j X64=-269.743+45.715 j X65=-192.181+31.0702 j X66=-148.672+22.9142 j X67=-120.84+17.7421 j X68=-101.515+14.1866 j X69=-87.3207+11.604 j X70=-76.4576+9.65156 j X71=-67.8801+8.13005 j X72=-60.9384+6.91586 j X73=-55.2075+5.92821j X74=-50.3979+5.11215 j X75=-46.3055+4.42897 j X76=-42.7823+3.85067j X77=-39.7182+3.35648 j X78=-37.0299+2.93069 j X79=-34.653+2.5612 j X80=-32.537+2.23853 j X81=-30.6418+1.95519 j X82=-28.9351+1.70515j X83=-27.3904+1.48351 j X84=-25
