实现一下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;
}
