c語言的FFT寫法

[給妳我所有的溫柔]

最近拿學長以前寫的2維FFT來修改為1維的轉換
可是發現一定要先做FFT才可以做IFFT
而不能先做IFFT後，再做FFT

我所測試的矩陣是一個大小為8內值為0，但picture[3]=1的矩陣
先做鏡射，成為result[16]，然後再做IFFT(1為FFT，-1為IFFT)
發現其結果不對

不知道有哪為高手可以幫我偵測一下，到底問題出在哪裡？
小弟真的是江郎才盡了
拜託拜託.....

底下為我修改後的1維轉換：
//---------------------------------------------------------------//

#include <iostream.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#define pi 3.141592654

const int N=8;
const int size=N*2;

class complex_var
{
public:
double real;
double imag;
};

complex_var f_uv[N];
double picture[N];
double result[size];

void fft(int forward_inverse,int size_x);

void main()
{
int freq_k=4,w=0;
for(int i=0;i<N;i++)
picture[i]=0;

picture[freq_k-1]=1;

for(i=0;i<N;i++)
{
w=(size-1)-i;
result[i]=picture[i];
result[w]=picture[i];
}

fft(-1,size);

fft(1,size);

}

void fft(int forward_inverse,int size_x)
{
int s_x,order_x;
int x,j,k;
int l,le,le1,ip;
double max_power,max_amp,temp;
complex_var u,w,t;

order_x=(int)(ceil(log(size_x)/log(2.0)));
s_x=(int)(pow(2,order_x));

if(forward_inverse==1)
{
for(x=0;x<s_x;x++)
{
f_uv[x].real=(double)result[x]*pow(-1,x);
f_uv[x].imag=0.0;
}
}
else
{
for(x=0;x<s_x;x++)
{
f_uv[x].imag=(-1.0)*f_uv[x].imag;
}
}

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

for(j=1,x=1;x<s_x;x++)
{
if(x<j)
{
t.real=f_uv[j-1].real;
t.imag=f_uv[j-1].imag;
f_uv[j-1].real=f_uv[x-1].real;
f_uv[j-1].imag=f_uv[x-1].imag;
f_uv[x-1].real=t.real;
f_uv[x-1].imag=t.imag;
}
k=s_x/2;
while(k<j)
{
j-=k;
k=k/2;
}
j+=k;
}

for(l=1;l<=order_x;l++)
{
le=pow(2,l);
le1=le/2;
u.real=1.0;
u.imag=0.0;
w.real=cos(pi/le1);
w.imag=-sin(pi/le1);
for(j=1;j<=le1;j++)
{
for(x=j;x<=s_x;x+=le)
{
ip=x+le1;
t.real=f_uv[ip-1].real*u.real-f_uv[ip-1].imag*u.imag;
t.imag=f_uv[ip-1].real*u.imag+f_uv[ip-1].imag*u.real;
f_uv[ip-1].real=f_uv[x-1].real-t.real;
f_uv[ip-1].imag=f_uv[x-1].imag-t.imag;
f_uv[x-1].real=f_uv[x-1].real+t.real;
f_uv[x-1].imag=f_uv[x-1].imag+t.imag;
}
t.real=w.real*u.real-w.imag*u.imag;
t.imag=w.real*u.imag+w.imag*u.real;
u.real=t.real;
u.imag=t.imag;
}
}

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

temp=s_x;

for(x=0;x<s_x;x++)
{
f_uv[x].real=f_uv[x].real/temp;
f_uv[x].imag=f_uv[x].imag/temp;
}

if(forward_inverse==1)
{
max_power=0;
for(x=0;x<s_x;x++)
{
temp=pow(f_uv[x].real,2)+pow(f_uv[x].imag,2);
if(temp>max_power)
max_power=temp;
}

max_amp=log10(1+sqrt(max_power));
for(x=0;x<s_x;x++)
{
temp=255*log10(1+sqrt(pow(f_uv[x].real,2)+pow(f_uv[x].imag,2)))/max_amp;
if(temp>=255)
temp=255;
result[x]=(double)temp;
}
}
else
{
for(x=0;x<size_x;x++)
{
f_uv[x].real=abs(f_uv[x].real);
if(f_uv[x].real>255)
f_uv[x].real=255;
result[x]=(double)(f_uv[x].real);
}
}
}
--
* Origin: ★ 交通大學資訊科學系 BBS ★ <bbs.cis.nctu.edu.tw: 140.113.23.3>