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++;
}
……
新闻热点
疑难解答