FFT快速入门

只看该作者 · 2011-12-25 12:51

FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。
虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。

现在圈圈就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就
不在此罗嗦了。

采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。

假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。
假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。
好了，说了半天，看着公式也晕，下面圈圈以一个实际的信号来做说明。

假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)
式中cos参数为弧度，所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？我们来看看FFT的结果的模值如图所示。

图1FFT结果
从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看：
1点：512+0i
2点：-2.6195E-14-1.4162E-13i
3点：-2.8586E-14-1.1898E-13i

50点：-6.2076E-13-2.1713E-12i
51点：332.55-192i
52点：-1.6707E-12-1.5241E-12i

75点：-2.2199E-13-1.0076E-12i
76点：3.4315E-12+192i
77点：-3.0263E-14+7.5609E-13i

很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0，即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。分别计算这三个点的模值，结果如下：
1点：512
51点：384
76点：192
按照公式，可以计算出直流分量为：512/N=512/256=2；50Hz信号的幅度为：384/(N/2)=384/(256/2)=3；75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。可见，从频谱分析出来的幅度是正确的。
然后再来计算相位信息。直流信号没有相位可言，不用管它。先计算50Hz信号的相位，atan2(-192,332.55)=-0.5236,结果是弧度，换算为角度就是180*(-0.5236)/pi=-30.0001。再计算75Hz信号的相位，atan2(192,3.4315E-12)=1.5708弧度，换算成角度就是180*1.5708/pi=90.0002。可见，相位也是对的。根据FFT结果以及上面的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。

总结：假设采样频率为Fs，采样点数为N，做FFT之后，某一点n（n从1开始）表示的频率为：Fn=(n-1)*Fs/N；该点的模值除以N/2就是对应该频率下的信号的幅度（对于直流信号是除以N）；该点的相位即是对应该频率下的信号的相位。相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi到pi。要精确到xHz，则需要采样长度为1/x秒的信号，并做FFT。要提高频率分辨率，就需要增加采样点数，这在一些实际的应用中是不现实的，需要在较短的时间内完成分析。解决这个问题的方法有频率细分法，比较简单的方法是采样比较短时间的信号，然后在后面补充一定数量的0，使其长度达到需要的点数，再做FFT，这在一定程度上能够提高频率分辨力。
具体的频率细分法可参考相关文献。

只看该作者 · 2011-12-25 13:02

哈!经过连续几个晚上的奋战,终于弄懂了FFT推导过程及实现! HappyJ
基2 FFT总的思想是将输入信号对半分割,再对半分割,再再对半分割(以下省略10000个再再...J)直至分割到2点.

两点DFT简化
假设输入为x[0],x[1];输出为X[0],X[1]. 伪代码如下
:

// ------------------------------------------------------------------

#define
N 2

#define
PI 3.1415926

// ------------------------------------------------------------------

int
i, j

for(i=0, X=0.0; i<N; i++)

for(j=0; j<N; j++)

   X += x[j] * ( cos(2*PI*i*j/N) - sin(2*PI*i*j/N) );

注意到(我想Audio编解码很多时候都是对cos,sin进行优化!)

[size=10.5000pt]

j=0[size=10.5000pt]

j=1[size=10.5000pt]

[size=10.5000pt]

i=0[size=10.5000pt]

cos 1[size=10.5000pt]
sin 0[size=10.5000pt]

[size=10.5000pt]
tw  1[size=10.5000pt]
cos 1[size=10.5000pt]
sin 0[size=10.5000pt]

[size=10.5000pt]
tw  1[size=10.5000pt]

[size=10.5000pt]

i=1[size=10.5000pt]

cos 1[size=10.5000pt]
Sin 0[size=10.5000pt]

[size=10.5000pt]
tw  1[size=10.5000pt]
cos -1[size=10.5000pt]
sin 0[size=10.5000pt]

[size=10.5000pt]
tw  -1[size=10.5000pt]

X[0]
=  x[0]*(1-0) + x[1]*(1-0) = x[0] + 1*x[1];
X[1]
=  x[0]*(1-0) + x[1]*(-1-0) = x[0] - 1*x[1];

这就是单个2点蝶形算法.

FFT实现流程图分析(N=8,以8点信号为例)

FFT implementation of an 8-point DFT as two 4-point DFTs and four 2-point DFTs

只看该作者 · 2011-12-25 13:03

file:///C:\DOCUME~1\WHOAMI~1\LOCALS~1\Temp\ksohtml\wps_clip_image-9595.pngfile:///C:\DOCUME~1\WHOAMI~1\LOCALS~1\Temp\ksohtml\wps_clip_image-9039.png

8点FFT流程图(Layer表示层, gr表示当前层的颗粒)

只看该作者 · 2011-12-25 13:03

下面以LayerI为例.

file:///C:\DOCUME~1\WHOAMI~1\LOCALS~1\Temp\ksohtml\wps_clip_image-9689.png

LayerI部分,具有4个颗粒,每个颗粒2个输入

(注意2个输入的来源,由时域信号友情提供,感谢感谢J)

我们将输入x[k]分为两部分x_r[k], x_i[k].具有实部和虚部,时域信号本没有虚部的,因此可以让x_i[k]为0.那么为什么还要画蛇添足分为实部和虚部呢?这是因为LayerII, LayerIII的输入是复数,为了编码统一而强行分的.当然你编码时可以判断当前层是否为1来决定是否分.但是我想每个人最后都会倾向分的.
旋转因子
tw = cos(2*PI*k/N)-j*sin(2*PI*k/N);也可以分为实部和虚部,令其为tw_r, tw_i;
则tw = tw_r - j*tw_i;

X[k] = (x_r[k] + j*x_i[k]) + (tw_r–j*tw_i) * (x_r[k+N/2]+j*x_i[k+N/2])
则

X_R[k] = x_r[k] + tw_r*x_r[k+N/2] + tw_i*x_i[k+N/2];

X_I[k] = x_i[k] - tw_i*x_r[k+N/2] + tw_r*x_i[k+N/2];

file:///C:\DOCUME~1\WHOAMI~1\LOCALS~1\Temp\ksohtml\wps_clip_image-25831.png

LayerII部分,具有2个颗粒,每个颗粒4个输入

(注意4个输入的来源,由LayerI友情提供,感谢感谢J)

只看该作者 · 2011-12-25 13:04

file:///C:\DOCUME~1\WHOAMI~1\LOCALS~1\Temp\ksohtml\wps_clip_image-9752.png

LayerIII部分,具有1个颗粒,每个颗粒8个输入

(注意8个输入的来源,由LayerII友情提供,感谢感谢J)

LayerI, LayerII, LayerIII从左往右,蝶形信号运算流非常明显!

假令输入为x[k], x[k+N/2],输出为X[k], X[k+N/2]. x[k]分解为x_r[k], x_i[k]部分
则该蝶形运算为
X[k]
= (x_r[k]-j*x_i[k]) + (x_r[k+N/2]-j*x_i[k+N/2])*(cos(2*PI*k/N)-j*sin(2*PI*k/N));
再令cos(2*PI*k/N)为tw1, sin(2*PI*k/N)为tw2则
X[k] = (x_r[k]-j*x_i[k]) + (x_r[k+N/2]-j*x_i[k+N/2])*(tw1-j*tw2);

X_R[k] = x_r[k] + x_r[k+N/2]*tw1 - x_i[k+N/2]*tw2;
X_I[K] = x_i[k]

x_r[k] = x_r[k] + x_r[k+b]*tw1 + x_i[k+b]*tw2;

x_i[k] = x_i[k] - x_r[k+b]*tw2 + x_i[k+b]*tw1;

譬如8点输入x[8]
1. 先分割成2部分: x[0], x[2], x[4], x[6]和 x[1], x[3], x[5], x[7]
2. 信号x[0], x[2], x[4], x[6]再分割成x[0], x[4]和 x[2], x[6]

信号x[1], x[3], x[5], x[7]再分割成x[1], x[5]和 x[3], x[7]
3.无法分割了,已经分割成2点了J.

如上图:
在LayerI的时候,我们是对2点进行DFT.(一共4次DFT )
输入为 x[0]&x[4]; x[2]&x[6]; x[1]&x[5]; x[3]&x[7]
输出为 y[0],y[1]; Y[2],y[3]; Y[4],y[5]; Y[6],y[7];

流程:
I.希望将输入直接转换为x[0], x[4], x[2], x[6], x[1], x[5], x[3], x[7]的顺序
II.对转换顺序后的信号进行4次DFT

步骤I代码实现

/**

*反转算法.
这个算法效率比较低!先用起来在说,之后需要进行优化.

*/

static
void
bitrev(void
)

{

int p=1, q, i;

int bit_rev[ N ];

float xx_r[ N ];

只看该作者 · 2011-12-25 13:04

bit_rev[ 0 ] = 0;

while( p < N )

{

for(q=0; q<p; q++)

{

bit_rev[ q ] = bit_rev[ q ] * 2;

bit_rev[ q + p ] = bit_rev[ q ] + 1;

}

p *= 2;

}

for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];

for(i=0; i<N; i++) x_r = xx_r[ bit_rev ];

}
// ------------------------此刻序列x重排完毕------------------------

步骤II代码实现
int j;
float TR; //临时变量
float tw1; //旋转因子

/*两点DFT */

for(k=0; k<N; k+=2)

{

//两点DFT简化告诉我们tw1=1

TR = x_r[k]; // TR就是A, x_r[k+b]就是B.

x_r[k]  = TR + tw1*x_r[k+b];

x_r[k+b] = TR - tw1*x_r[k+b];

}

在LayerII的时候,我们希望得到z,就需要对y进行DFT.
y[0],y[2]; y[1],y[3]; y[4],y[6]; y[5],y[7];
z[0],       z[1]; z[2],z[3]; z[4],z[5]; z[6],z[7];

在LayerIII的时候,我们希望得到v,就需要对z进行DFT.
z[0],z[4]; z[1],z[5]; z[2],z[6]; z[3],z[7];
v[0],v[1]; v[2],v[3]; v[4],v[5]; v[6],v[7];

准备

令输入为x~~, x[s+N/2],输出为y~~, y[s+N/2]~~~~

这个N绝对不是上面的8,这个N是当前颗粒的输入样本总量

对于LayerI而言N是2;对于LayerII而言N是4;对于LayerIII而言N是8

复数乘法:(a+j*b) * (c+j*d)

实部= a*c – bd;

虚部= ad + bc;

旋转因子:

实现(C描述)

#include
<stdio.h>

#include
<math.h>

#include
<stdlib.h>

//#include "complex.h"

// --------------------------------------------------------------------------

#define  N 8 //64

#define  M 3 //6 //2^m=N

#define  PI 3.1415926

// --------------------------------------------------------------------------

float  twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};

float  x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};

float  x_i[N];                      //N=8

/*

float  twiddle[N/2] = {1,    0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,

                     0.7075, 0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,

                     0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,

                  -0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951};    //N=64

float  x_r[N]={1,1,1,1,1,1,1,1,

            1,1,1,1,1,1,1,1,

            1,1,1,1,1,1,1,1,

                           1,1,1,1,1,1,1,1,

                           0,0,0,0,0,0,0,0,

                           0,0,0,0,0,0,0,0,

                           0,0,0,0,0,0,0,0,

                           0,0,0,0,0,0,0,0,};

float  x_i[N];

*/

FILE *fp;

// ----------------------------------- func -----------------------------------

/**

*初始化输出虚部

*/

static
void
fft_init(void
)

{

int
i;

for(i=0; i<N; i++)  x_i = 0.0;

}

/**

*反转算法.将时域信号重新排序.

*这个算法有改进的空间

*/

static
void
bitrev(void
)

{

int p=1, q, i;

int bit_rev[ N ]; //

float  xx_r[ N ]; //

bit_rev[ 0 ] = 0;

while( p < N )

{

for(q=0; q<p; q++)

   {

      bit_rev[ q ] = bit_rev[ q ] * 2;

      bit_rev[ q + p ] = bit_rev[ q ] + 1;

   }

   p *= 2;

}

for(i=0; i<N; i++)  xx_r[ i ] = x_r[ i ];

for(i=0; i<N; i++)  x_r = xx_r[ bit_rev ];

}

/* ------------ add by sshc625 ------------ */

static
void
bitrev2(void
)

{

return
;

}

/* */

void
display(void
)

{

printf("\n\n");

int  i;

for(i=0; i<N; i++)

   printf("%f\t%f\n", x_r, x_i);

}

/**

*

*/

void
fft1(void
)

{ fp = fopen("log1.txt","a+");

int L, i, b, j, p, k, tx1, tx2;

float  TR, TI, temp; //临时变量

float  tw1, tw2;

/*深M.对层进行循环. L为当前层,总层数为M. */

for(L=1; L<=M; L++)

{

   fprintf(fp,"----------Layer=%d----------\n", L);

/* b的意义非常重大,b表示当前层的颗粒具有的输入样本点数 */

   b = 1;

   i = L - 1;

while(i > 0)

   {

      b *= 2;

      i--;

   }

// --------------是否外层对颗粒循环,内层对样本点循环逻辑性更强一些呢! --------------

/*

   * outter对参与DFT的样本点进行循环

   * L=1,循环了1次(4个颗粒,每个颗粒2个样本点)

   * L=2,循环了2次(2个颗粒,每个颗粒4个样本点)

   * L=3,循环了4次(1个颗粒,每个颗粒8个样本点)

   */

for(j=0; j<b; j++)

   {

/*求旋转因子tw1 */

      p = 1;

      i = M - L; // M是为总层数, L为当前层.

while(i > 0)

      {

         p = p*2;

         i--;

      }

      p  = p * j;

      tx1 = p % N;

      tx2 = tx1 + 3*N/4;

      tx2 = tx2 % N;

// tw1是cos部分,实部; tw2是sin部分,虚数部分.

      tw1 = ( tx1>=N/2)? -twiddle[tx1-N/2]  : twiddle[ tx1 ];

      tw2 = ( tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];

/*

         * inner对颗粒进行循环

         * L=1,循环了4次(4个颗粒,每个颗粒2个输入)

         * L=2,循环了2次(2个颗粒,每个颗粒4个输入)

         * L=3,循环了1次(1个颗粒,每个颗粒8个输入)

         */

for(k=j; k<N; k=k+2*b)

      {

         TR = x_r[k]; // TR就是A, x_r[k+b]就是B.

         TI = x_i[k];

         temp = x_r[k+b];

/*

            *如果复习一下
(a+j*b)(c+j*d)两个复数相乘后的实部虚部分别是什么

            *就能理解为什么会如下运算了,只有在L=1时候输入才是实数,之后层的

            *输入都是复数,为了让所有的层的输入都是复数,我们只好让L=1时候的

            *输入虚部为0

            * x_i[k+b]*tw2是两个虚数相乘

            */

         fprintf(fp,"tw1=%f, tw2=%f\n", tw1, tw2);

         x_r[k]  = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;

         x_i[k]  = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;

         x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;

         x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;

只看该作者 · 2011-12-25 13:05

fprintf(fp,"k=%d, x_r[k]=%f, x_i[k]=%f\n", k, x_r[k], x_i[k]);

fprintf(fp,"k=%d, x_r[k]=%f, x_i[k]=%f\n", k+b, x_r[k+b], x_i[k+b]);

}//

}

/**

* ------------ add by sshc625 ------------

*该实现的流程为

* for( Layer )

* for( Granule )

* for( Sample )

*

*/

void
fft2(void
)

{ fp = fopen("log2.txt","a+");

int cur_layer, gr_num, i, k, p;

float tmp_real, tmp_imag, temp; //临时变量,记录实部

float tw1, tw2;//旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.

int
step; //步进

int
sample_num; //颗粒的样本总数(各层不同,因为各层颗粒的输入不同)

/*对层循环
*/

for(cur_layer=1; cur_layer<=M; cur_layer++)

{

/*求当前层拥有多少个颗粒(gr_num) */

gr_num = 1;

i = M - cur_layer;

while(i > 0)

{

i--;

gr_num *= 2;

}

/*每个颗粒的输入样本数N' */

sample_num = (int)pow(2, cur_layer);

/*步进.步进是N'/2 */

step = sample_num/2;

/* */

k = 0;

/*对颗粒进行循环
*/

for(i=0; i<gr_num; i++)

{

/*

*对样本点进行循环,注意上限和步进

*/

for(p=0; p<sample_num/2; p++)

{

//旋转因子,需要优化...

tw1 = cos(2*PI*p/pow(2, cur_layer));

tw2 = -sin(2*PI*p/pow(2, cur_layer));

tmp_real = x_r[k+p];

tmp_imag = x_i[k+p];

temp = x_r[k+p+step];

/*(tw1+jtw2)(x_r[k]+jx_i[k])

*

* real : tw1*x_r[k] - tw2*x_i[k]

* imag : tw1*x_i[k] + tw2*x_r[k]

*我想不抽象出一个

* typedef struct {

* double real; //实部

* double imag; //虚部

* } complex;以及针对complex的操作

*来简化复数运算是否是因为效率上的考虑!

*/

/*蝶形算法
*/

x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );

x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );

/* X[k] = A(k)+WB(k)

* X[k+N/2] = A(k)-WB(k)的性质可以优化这里*/

//旋转因子,需要优化...

tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));

tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));

x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );

x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );

printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);

printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);

}

/*开跳!:) */

k += 2*step;

}

/*

*后记:

*究竟是颗粒在外层循环还是样本输入在外层,好象也差不多,复杂度完全一样.

*但以我资质愚钝花费了不少时间才弄明白这数十行代码.

*从中我发现一个于我非常有帮助的教训,很久以前我写过一部分算法,其中绝大多数都是递归.

*将数据量减少,减少再减少,用归纳的方式来找出数据量加大代码的规律

*比如FFT

* 1.先写死LayerI的代码;然后再把LayerI的输出作为LayerII的输入,又写死代码; ......

* 大约3层就可以统计出规律来.这和递归也是一样,先写死一两层,自然就出来了!

* 2.有的功能可以写伪代码,不急于求出结果,降低复杂性,把逻辑结果定出来后再添加.

* 比如旋转因子就可以写死,就写1.0.流程出来后再写旋转因子.

*寥寥数语,我可真是流了不少汗! Happy!

*/

void
dft(void
)

{

int i, n, k, tx1, tx2;

float tw1,tw2;

float xx_r[N],xx_i[N];

/*

* clear any data in Real and Imaginary result arrays prior to DFT

*/

for(k=0; k<=N-1; k++)

xx_r[k] = xx_i[k] = x_i[k] = 0.0;

// caculate the DFT

for(k=0; k<=(N-1); k++)

{

for(n=0; n<=(N-1); n++)

{

tx1 = (n*k);

tx2 = tx1+(3*N)/4;

tx1 = tx1%(N);

tx2 = tx2%(N);

if(tx1 >= (N/2))

tw1 = -twiddle[tx1-(N/2)];

else

tw1 = twiddle[tx1];

if(tx2 >= (N/2))

tw2 = -twiddle[tx2-(N/2)];

else

tw2 = twiddle[tx2];

xx_r[k] = xx_r[k]+x_r[n]*tw1;

xx_i[k] = xx_i[k]+x_r[n]*tw2;

}

xx_i[k] = -xx_i[k];

}

// display

for(i=0; i<N; i++)

printf("%f\t%f\n", xx_r, xx_i);

}

// ---------------------------------------------------------------------------

int
main(void
)

{

fft_init( );

bitrev( );

// bitrev2( );

//fft1( );

fft2( );

display( );

system("pause"
);

// dft();

return
1;

}