实现一下FFT 与 IFFT 反变换

 //发现网上有好多FFT变换，感觉有一半都是互抄的，

  //那我也抄一下，加入了自己的理解！应该有部分与网上实现有点点不同！

  // 实际上FFT 与IFFT 可以在一个函数中实现就可以了！不过不想改了！

#define PI (3.1415926535)

typedef struct _plural
{
double real;
double img;
} plural;

plural EE( plural b1, plural b2)
{
plural b3;
b3.real=b1.real*b2.real- b1.img*b2.img;
b3.img=b1.real*b2.img+ b1.img*b2.real;
return(b3);
}

void iOrder(plural *pData,int levelN)
{
plural tmp;
int k,j=1;

//变变址运算，实现数据运算前的位倒序
for(int i=1;i< levelN;i++)
{

if(i>j)
{
tmp=pData[j-1];
pData[j-1]=pData[i-1];
pData[i-1]=tmp;
}

k=levelN/2;
while(k<j)
{
j=j-k;
k=k/2;
}
j=j+k;
}

}

plural *pFftW = NULL;

void FFT_Init(int N)
{
int i;
pFftW=new plural[N];

for(i=0;i<N;i++)
{
pFftW[i].real=cos(2*i*PI/N);
pFftW[i].img=-sin(2*i*PI/N);
printf("%4d %11.7f %11.7f \n", i, pFftW[i].real, pFftW[i].img);

}

}

void FFT_DFT(plural *pData, int N)
{
int m = 0;

int f = N;
for(m=1;(f=f/2)!=1;m++) ;

printf("\n====================\n");
printf("=== 倒序 \n");

// 倒序
iOrder(pData,N);

plural v,w;
int len = 1;
for(int k = 0; k<m; k++) // 一共得计算多少次
{
len = len*2;

printf("\n====================\n");
printf("=== 【 %d 】 %d %d \n", k, len);

v.real=1.0; v.img=0.0;

int tmp = len/2;

//此处是把360度 划分成 2^m 份 -->此处只取一份
w.real=cos(PI/tmp);
w.img=-sin(PI/tmp);

for (int j = 0; j<len/2; j++)
{
for (int i=0; i< N; i+=len)
{
//取出二对复数
plural tmp1 = pData[i+j];
plural tmp2 = pData[i+j + len/2];

// *W(0,N)
// 复数相乘 (a+bj)(c+dj) = (ac-bd) + (bc+ad)j
plural t=EE(pData[i+j + len/2],v);

pData[i+j + len/2].real = pData[i+j].real - t.real;
pData[i+j + len/2].img = pData[i+j].img - t.img;
pData[i+j].real = pData[i+j].real + t.real;
pData[i+j].img = pData[i+j].img + t.img;

}
v=EE(v,w); //此处实际上是把Z坐标角度相加
}

}
}

void IFFT_DFT(plural *pData, int N)
{
int m = 0;

int f = N;
for(m=1;(f=f/2)!=1;m++) ;

printf("\n====================\n");
printf("=== 倒序 \n");

// 倒序
iOrder(pData,N);

plural v,w;
int len = 1;
for(int k = 0; k<m; k++) // 一共得计算多少次
{
len = len*2;

printf("\n====================\n");
printf("=== 【 %d 】 %d %d \n", k, len);

v.real=1.0; v.img=0.0;

int tmp = len/2;

//此处是把360度 划分成 2^m 份 -->此处只取一份
w.real=cos(-PI/tmp);
w.img=-sin(-PI/tmp);

for (int j = 0; j<len/2; j++)
{
for (int i=0; i< N; i+=len)
{
//取出二对复数
plural tmp1 = pData[i+j];
plural tmp2 = pData[i+j + len/2];

// *W(0,N)
// 复数相乘 (a+bj)(c+dj) = (ac-bd) + (bc+ad)j
plural t=EE(pData[i+j + len/2],v);

pData[i+j + len/2].real = (pData[i+j].real - t.real);
pData[i+j + len/2].img = (pData[i+j].img - t.img);
pData[i+j].real = (pData[i+j].real + t.real);
pData[i+j].img = (pData[i+j].img + t.img);

}
v=EE(v,w); //此处实际上是把Z坐标角度相加
}

for (int i=0; i< N; i++)
{
pData[i].real /= 2;
pData[i].img /= 2;
}

}

}

#define FFT_N (8)
plural pData[FFT_N];

int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])
{
int nRetCode = 0;

// 初始化 MFC 并在失败时显示错误
if (!AfxWinInit(::GetModuleHandle(NULL), NULL, ::GetCommandLine(), 0))
{
// TODO: 更改错误代码以符合您的需要
_tprintf(_T("错误: MFC 初始化失败\n"));
return 1;
}

//pData = new plural[FFT_N];
for (int i = 0; i<FFT_N; i++)
{
pData[i].real = i*1.0;
pData[i].img = 0;
}

for (int i = 0; i<FFT_N; i++)
{
printf("%4d ---> %7.3f %7.3f \n", i, pData[i].real, pData[i].img);
}
printf("\n");

FFT_Init(FFT_N);
FFT_DFT(pData,FFT_N);

printf("\nFFT 输出\n");

for (int i = 0; i<FFT_N; i++)
{
printf("%4d ---> %7.3f %7.3f \n", i, pData[i].real, pData[i].img);
}

IFFT_DFT(pData,FFT_N);

printf("\nIFFT 输出\n");

for (int i = 0; i<FFT_N; i++)
{
printf("%4d ---> %7.3f %7.3f \n", i, pData[i].real, pData[i].img);
}

// delete pData;

system("pause");

return nRetCode;
}