實現一下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;
}