typedef strUCt tagComplex{
float Re; //復數的實部
float Im; //復數的虛部
}Complex;int NV2=N/2;
int NM1=N-1;
int I,J,K=0;
Complex T;//用于中介的復數變量T
I=J=1;
while(I<=NM1)
{
if(I<J)
{
T=A[J-1];//將A[J-1]的內容和A[I-1]的內容互換，借助于中間變量T
A[J-1]=A[I-1];
A[I-1]=T;
}
K=NV2;
while(K<J)
{
J-=K;
K/=2;
}
J+=K;
I++;
}……
Complex U,W,T;
int LE,LE1,I,J,ip;
int N=(int)pow(2,M);
//在此采用的是時間抽選奇偶分解方式，所以在參加運算前首先要對時間序列進行倒序
ReverSEOrder(A,N);
int L=1;
while(L<=M)
{
LE=(int)pow(2,L);
LE1=LE/2;
U.Re=1.0f;
U.Im=0.0f;
W.Re=(float)cos(PI/(1.0*LE1));//計算W算子的值
W.Im=(float)-1.0*sin(PI/(1.0*LE1));
if(abs(W.Re)<1.0e-12)
W.Re=0.0f;
if(abs(W.Im)<1.0e-12)
W.Im=0.0f;
J=1;
while(J<=LE1)
{
I=J;
while(I<=N)
{
IP=I+LE1;
T.Re=(float)A[IP-1].Re*U.Re-A[IP-1].Im*U.Im;//計算復數運算A*U
T.Im=(float)A[IP-1].Re*U.Im+A[IP-1].Im*U.Re;
A[IP-1].Re=(float)A[I-1].Re-T.Re;//計算復數運算A-T
A[IP-1].Im=(float)A[I-1].Im-T.Im;
A[I-1].Re+=T.Re;//計算復數運算A+T
A[I-1].Im+=T.Im;
I+=LE;
}
float temp=U.Re;
U.Re=(float)U.Re*W.Re-U.Im*W.Im;//計算復數運算U*W
U.Im=(float)temp*W.Im+U.Im*W.Re;
J++;
}
L++;
}
……新聞熱點
疑難解答
