求fft代码；

[lanren_me]

手头有部分代码：但是不知是干吗的。所以可以的话来点注释;

[dcyu]

Fast Fourier Transform 是信号处理常用的简化计算的算法

[dcyu]

72. 函数名：_FFT

功能：FFT变换

用法：void _FFT(double *fr, double *fi, int n, int flag)

参数：fr：采样点的实数表，fi：采样点的虚数表，n：采样点个数，flag：flag=0表示求Fourier变换，flag=1表示求逆Fourier变换

返回：无

void _FFT(double *fr, double *fi, int n, int flag)
{
int mp,arg,q,cntr,p1,p2;
int i,j,a,b,k,m;
double sign,pr,pi,harm,x,y,t;
double *ca,*sa;

ca=(double *)calloc(n,sizeof(double));
sa=(double *)calloc(n,sizeof(double));
if(ca==NULL||sa==NULL) Error("calloc error in _FFT!");
j=0;
if(flag!=0)
{
sign=1.0;
for(i=0;i<=n-1;i++)
{
fr[i]=fr[i]/n;
fi[i]=fi[i]/n;
}
}
else
sign=-1.0;
for(i=0;i<=n-2;i++)
{
if(i<j)
{
t=fr[i];
fr[i]=fr[j];
fr[j]=t;
t=fi[i];
fi[i]=fi[j];
fi[j]=t;
}
k=n/2;
while(k<=j)
{
j-=k;
k/=2;
}
j+=k;
}
mp=0;
i=n;
while(i!=1)
{
mp+=1;
i/=2;
}
harm=2*M_PI/n;
for(i=0;i<=n-1;i++)
{
sa[i]=sign*sin(harm*i);
ca[i]=cos(harm*i);
}
a=2;
b=1;
for(cntr=1;cntr<=mp;cntr++)
{
p1=n/a;
p2=0;
for(k=0;k<=b-1;k++)
{
i=k;
while(i<n)
{
arg=i+b;
if(k==0)
{
pr=fr[arg];
pi=fi[arg];
}
else
{
pr=fr[arg]*ca[p2]-fi[arg]*sa[p2];
pi=fr[arg]*sa[p2]+fi[arg]*ca[p2];
}
fr[arg]=fr[i]-pr;
fi[arg]=fi[i]-pi;
fr[i]+=pr;
fi[i]+=pi;
i+=a;
}
p2+=p1;
}
a*=2;
b*=2;
}
free(ca);
free(sa);
}

[dcyu]

在频域方面会有复频率出现，因此就需要建立n个复数的抽样点，
实数表，虚数表是这n个点各自的实部和虚部。
比如四个采样点为：1+j2 , 2+j3 , 4+j5 , 67+j12
它们的实数表，虚数表是两个一维数组：
{1,2,4,67},{2,3,5,12}

[lanren_me]

能否解释一下：实数表，虚数表

[eastsun]

void HPFAST_OneDimension_FFT_Trans(
_complex_digital_array the_src_sign_array,
_complex_digital_array& result_sign_array )
{
_complex_digital_array src_sign_array;
src_sign_array = the_src_sign_array;

int nInputSignAmount = src_sign_array.GetSize( );
result_sign_array.RemoveAll( );

if( ! Is2MiNumber( nInputSignAmount ) )
return;
int nInputSignLog2Mi = GetLog2( nInputSignAmount );
if( nInputSignLog2Mi <= 0 )
{
result_sign_array = src_sign_array;
return;
}

result_sign_array.InitBySize( nInputSignAmount );
////////////////////////////////////////////////////////////

int nMiIndex = 0;
////////////////////////////////////////////////////////////
int nDotIndex = 0;
for( nDotIndex = 0; nDotIndex < nInputSignAmount; nDotIndex ++ )
{
int m = nDotIndex;
int is = 0;
for( nMiIndex = 0; nMiIndex < nInputSignLog2Mi; nMiIndex ++ )
{
int j = m / 2;
is = 2 * is + ( m - 2 * j );
m = j;
}
_complex_digital temp_des_digital;
_complex_digital temp_src_digital;
temp_src_digital = src_sign_array.GetAt( is );
temp_des_digital = temp_src_digital;
result_sign_array.SetAt( nDotIndex, temp_des_digital );
}//这个循环完成置换。

double dGuiyiCoefficient = 2.0 * PI / nInputSignAmount;
_complex_digital temp_digital0;
temp_digital0.HP_SetRealAndVirtualValue( 1.0, 0.0 );

_complex_digital temp_digital1;
temp_digital1.HP_SetRealAndVirtualValue( cos( dGuiyiCoefficient ), - sin( dGuiyiCoefficient ) );
{
src_sign_array.SetAt( 0, temp_digital0 );
src_sign_array.SetAt( 1, temp_digital1 );
}

_complex_digital base_trans_digital = temp_digital1;
for( nDotIndex = 2; nDotIndex < nInputSignAmount; nDotIndex ++ )
{
_complex_digital the_pre_digital = src_sign_array.GetAt( nDotIndex - 1 );

_complex_digital current_digital;
current_digital = the_pre_digital * base_trans_digital;

src_sign_array.SetAt( nDotIndex, current_digital );
}//这个循环完成基本的单点变换。

for( nDotIndex = 0; nDotIndex < nInputSignAmount - 1; nDotIndex += 2 )
{
_complex_digital current_digital = result_sign_array.GetAt( nDotIndex );
_complex_digital next_digital = result_sign_array.GetAt( nDotIndex + 1 );

_complex_digital origen_current_digital = current_digital;
current_digital += next_digital;
next_digital = origen_current_digital - next_digital;

result_sign_array.SetAt( nDotIndex, current_digital );
result_sign_array.SetAt( nDotIndex + 1, next_digital );
}

int m = nInputSignAmount / 2;
int nv = 2;
for( nMiIndex = nInputSignLog2Mi - 2; nMiIndex >= 0; nMiIndex -- )
{
m = m / 2;
nv = 2 * nv;
for( nDotIndex = 0; nDotIndex <= ( m - 1 ) * nv; nDotIndex += nv )
{
int jjs = 0;
for( jjs = 0; jjs < ( nv / 2 ); jjs ++ )
{
int nOffset = nDotIndex + jjs + nv / 2;
_complex_digital temp_src_mj_digital = src_sign_array.GetAt( m * jjs );
_complex_digital temp_result_offset_digital = result_sign_array.GetAt( nOffset );

_complex_digital the_product_digital = temp_src_mj_digital * temp_result_offset_digital;

_complex_digital the_result_index_add_jjs_digital = result_sign_array.GetAt( nDotIndex + jjs );

temp_result_offset_digital = the_result_index_add_jjs_digital - the_product_digital;
the_result_index_add_jjs_digital += the_product_digital;

result_sign_array.SetAt( nOffset, temp_result_offset_digital );
result_sign_array.SetAt( nDotIndex + jjs, the_result_index_add_jjs_digital );
}
}
}
result_sign_array /= ( 1.0 * nInputSignAmount );
}

[strafer]

wa, 楼上大哥的程序风格像是大公司的， 变量都选那么长的！！！

[lanren_me]

大大：

问题是：比如：采样N个点，每个点的值为float值;

那么这个float值是实数还是虚数？？？怎么转化为复数啊???