dif-fft.cpp

计算序列的DFT的快速算法

# include <vector>
# include <iostream>
# include <complex>
# include <math.h>
# include <iomanip>
using namespace std;
const double Pi=acos(-1.0);
int rev(int index,int NumberOfBit)
{
        int ans=0;
        for(int i=0;i<NumberOfBit;i++,index>>=1) ans=(ans<<1)|(index&1);
        return ans;
}
int NBits(int TOT)
{
        int count=0;
        while(TOT){
                TOT>>=1;
                count++;
        }
        return count-1;
}
vector<complex<double> > FFT(vector<complex<double> > a)
{
        int n=a.size();
        int k,m,j;
        if(n==1) return a;
        complex<double> w_m,w,t,u;
        vector<complex<double> >A(n);
        for(k=0;k<n;k++) A[rev(k,NBits(n))]=a[k];
        for(m=2;m<=n;m*=2){
                w_m=complex<double>(cos(2*Pi/m),-sin(2*Pi/m));
                for(k=0;k<n;k+=m){
                        w=1;
                        for(j=0;j<m/2;j++) {
                                t=w*A[k+j+m/2];
                                u=A[k+j];
                                A[k+j]=u+t;
                                A[k+j+m/2]=u-t;
                                w*=w_m;
                        }
                }
        }
        return A;
}
void print(vector<complex<double> >a)
{
	for(int i=0;i<a.size();i++) cout<<a[i]<<endl;
}
int main()
{
	vector<complex<double> > a;
	double tmp;
	cout<<"Note: the number of input real numbers must be power of 2!\n      a number >=1000 is the end."<<endl;
	//freopen("input_f.txt","r",stdin);
	while(cin>>tmp,fabs(tmp)<1000)
		a.push_back(complex<double>(tmp,0));
	vector<complex<double> > ans=FFT(a);
	print(ans);
	return 0;
}