fft_ifft.h

数值计算利用fft分解法来求解线性方程组,简明方便

//说明：MARK的取值1为FFT，-1为IFFT
void changeorder(complex *p,int M) 
{ int LH,N,N1,i,j,k;float T;
N=pow(2,M);
LH=N/2;j=LH;N1=N-2;
for(i=1;i<=N1;i++)
{if(i<j)
{
T=(p+i)->real;(p+i)->real=(p+j)->real;(p+j)->real =T;
T=(p+i)->imag;(p+i)->imag=(p+j)->imag;(p+j)->imag =T;
}
k=LH;
while(j>=k)
{
	j=j-k;
	k=(int)(k/2+0.5);
}
j=j+k;
}
}
void FFT_IFFT(complex *p,int M,int mark)
{
    int N,l,b,h,k,i,j;
    N=pow(2,M);
	complex *w=new complex[N];	
    float pI=8.0*atan(1.0);complex m;
    changeorder(p,M);
    for(i=0;i<N;i++)                              
	{
	(w+i)->real=cos(i*pI/N);
	(w+i)->imag=sin(i*pI/N)*mark;
	}
    for(l=1;l<=M;l++)                          
	{
	b=(int)(pow(2,(l-1))+0.5);
	for(j=0;j<=b-1;j++)
	{
		h=j*((int)(pow(2,(M-l))+0.5));
		for(k=j;k<=N-1;k+=(int)(pow(2,l)+0.5))
		{m.real=((p+k+b)->real)*((w+h)->real)-((p+k+b)->imag)*((w+h)->imag);
         m.imag=((p+k+b)->imag)*((w+h)->real)+((p+k+b)->real)*((w+h)->imag);
		 (p+k+b)->real=(p+k)->real-m.real;
		(p+k+b)->imag=(p+k)->imag-m.imag;
		 (p+k)->real=(p+k)->real+m.real;
		(p+k)->imag=(p+k)->imag+m.imag;
		}
	}
}
delete w;
}