main.cc

计算f t)=exp -|t|)的FFT

#include <iostream>
#include "matrix.h"
#include "main.h"
#include <math.h>
#include <time.h>
using namespace std;
char *out;
Complex *a=new Complex[N];

int parseCommandLine(int argc, char *argv[]){
    if(argc<2){      //没有参数，给出程序运行方法及相关信息 
         cout << "程序作者：王攀文" << endl;
         cout << "学号：S200715025" << endl; 
         cout << "程序信息：求f(t)=e^(-|t|)的FFT" << endl;
         cout << "程序运行方式：" << endl;
         cout << "程序文件名 <输出文件名>" << endl; 
         return 0;
    } 
	out=argv[1];
}

double fx(double x){                    	//原函数 
	return exp(-fabs(x));
}

void getA(FILE *fp){
	fprintf(fp,"原始数据:\n");
	for(int m=0;m<N;m++){
		a[m]=(fx(m*DELTA_T)+fx((m-N)*DELTA_T))*DELTA_T;
		fprintf(fp,"%4d: ",m+1);
		a[m].fprint(fp);
		fprintf(fp,"\t");
		if((m-NUM_PER_LINE+1)%NUM_PER_LINE==0)fprintf(fp,"\n");
	}
}

int converse(int src,int size){       	//进行二进制换位     
	int tmp=src; 
	int des=0; 
	for (int i=size-1; i>=0; i--){ 
	   des=((tmp & 0x1) << i) | des; 
	   tmp=tmp >> 1; 
	} 
	return des; 
}

Complex *convArr(Complex *a,int n){		// 得到二进制换位后的数据数组 
	int i,b;
	float t=log(n)/log(2);
	b=int(t);
	Complex *temp=new Complex[n];
	for(i=0;i<n;i++){
		temp[converse(i,b)]=a[i];
	}
	return temp;
}

double dft(){
	double start=clock();
	Complex *C=new Complex[N];
	Complex *W=new Complex[N*N];
	Matrix *T=new Matrix(N);
	int m,n;
    
    W[0]=1;     //计算w^0～w^N*N-1 
    W[1]=Complex(cos(_2PI_/(-1*N)),sin(_2PI_/(-1*N)));
    for(m=1;m<=(N-1)*(N-1);m++){
		W[m]=Complex(W[m-1].getReal()*W[1].getReal()-W[m-1].getImage()*W[1].getImage(),
					W[1].getImage()*W[m-1].getReal()+W[1].getReal()*W[m-1].getImage());
	}  

	for(m=0;m<N;m++){//得到　T 
		for(n=0;n<N;n++){
			T->setElement(W[m*n],m,n);
		}
	}

	C=(*T)*a;
	
	return double(clock() - start)/CLOCKS_PER_SEC;
}

double fft(FILE *fp){
	double start=clock();
	//变量声明 
	int m,n,nn,mn,k;
	
	Complex *C=new Complex[N];
	Complex *A=new Complex[N];	 
	Complex *aa=new Complex[N];
	aa=convArr(a,N);//换位,得到A' (aa) 
	
    //开始变换 
	for(m=1;m<=M;m++){
		nn=int(pow(2,m));
		Complex *W=new Complex[nn];
		W[0]=1;
		W[1]=Complex(cos(_2PI_/(-1*nn)),sin(_2PI_/(-1*nn)));
		for(n=2;n<nn/2;n++){
			W[n]=Complex(W[n-1].getReal()*W[1].getReal()-W[n-1].getImage()*W[1].getImage(),
						W[1].getImage()*W[n-1].getReal()+W[1].getReal()*W[n-1].getImage());
		}
		
        for(n=0;n<N/nn;n++){
      		for(mn=n*nn;mn<n*nn+nn;mn++){
				if(mn<n*nn+nn/2)A[mn]=aa[mn]+aa[mn+nn/2]*W[mn-n*nn];	
				else A[mn]=aa[mn-nn/2]-aa[mn]*W[mn-n*nn-nn/2];
			}
        }
        for(n=0;n<N;n++)aa[n]=A[n];
        delete [] W;
	}
	for(m=0;m<N;m++)C[m]=aa[m];
		
	//输出到文件
	fprintf(fp,"\n\nFFT 变换后的数据:\n");
	for(m=0;m<N;m++){
		fprintf(fp,"%4d: ",m+1);
		C[m].fprint(fp);
		fprintf(fp,"\t");
		if((m-NUM_PER_LINE+1)%NUM_PER_LINE==0)fprintf(fp,"\n");
	} 	
	return double(clock() - start)/CLOCKS_PER_SEC;
}

int main(int argc, char *argv[]){
	//从程序参数里得到输出文件地址 
	FILE *fp_dft,*fp_fft;
	if(!parseCommandLine(argc,argv)) return 0; 
	if(!(fp_fft=fopen(out,"w"))){
		cerr << "Cannot open output file: " << out << endl;
		return 0;
	}
	getA(fp_fft);
	
	//dft
	cout << "DFT 程序运行时间: " << dft()  << " 秒" << endl;
	
	//fft	
	cout << "FFT 程序运行时间: " << fft(fp_fft)  << " 秒" << endl;
	
	return 0;
}