fft.cpp

DSP图像处理程序源码

#include "stdafx.h"
#include "fft.h"
int CFFT::IS2N(int N){									//judge if N is 2^n
	if(   (N&(N-1))==0   ) 
        return 1;
else 
        return 0;
}

void CFFT::changeOrder(complex *a,complex *b,int N){
		int r,i,j,p;
		r=(int)(log(N+1.0)/log(2));             //power
		for (j=0;j<N;j++){
			p=0;
			for (i=0;i<r;i++){
				if (j&(1<<i))
					p+=1<<(r-i-1);
			}
			*(a+j)=*(b+p);
		}
		
}

void CFFT::FFT(double *xn,complex *Xk,int N){
	if (!IS2N(N))
		CDFT::DFT(xn,Xk,N);
	else FFT_RUN(xn,Xk,N);
}
void CFFT::IFFT(complex *Xk,double *xn,int N){
	if (!IS2N(N))
		CDFT::IDFT(Xk,xn,N);
	else IFFT_RUN(Xk,xn,N);
}

void CFFT::FFT_RUN(double *xn,complex *Xk,int N){
	int thisTime,totalTime,pointDistance,rowNo,lineNo,p,i;
	complex wn,xtemp;
	complex *Xtemp=(complex *)malloc(sizeof(complex)*N);
	for (i=0;i<N;i++)
		Xtemp[i]=complex(xn[i]);
	complex *Xn=(complex *)malloc(sizeof(complex)*N);
	changeOrder(Xn,Xtemp,N);
	totalTime=(int)(log(N+1)/log(2));//The total time of butterfly
	wn=exp(-complex(0.0,1.0)*2.0*M_PI/(double)N);
	for (thisTime=1;thisTime<=totalTime;thisTime++){		//This time is now butterfly time	
			pointDistance=(1<<(thisTime-1));				//the distance between butterfly
			for (rowNo=0;rowNo<pointDistance;rowNo++){
				p=(1<<(totalTime-thisTime))*rowNo;
				for (lineNo=rowNo;lineNo<N;lineNo+=(1<<thisTime)){
					xtemp=Xn[lineNo];
					Xn[lineNo]=Xn[lineNo]+Xn[lineNo+pointDistance]*(wn^((double)p));
					Xn[lineNo+pointDistance]=xtemp-Xn[lineNo+pointDistance]*(wn^((double)p));
				}
			}
	}
	for (i=0;i<N;i++)
		Xk[i]=Xn[i];

}			
	
void CFFT::IFFT_RUN(complex *Xk,double *xn,int N){
	int thisTime,totalTime,pointDistance,rowNo,lineNo,p,i;
	complex wn,xtemp;
	complex *Xn=(complex *)malloc(sizeof(complex)*N);
	complex *Xtemp=(complex *)malloc(sizeof(complex)*N);
	for (i=0;i<N;i++)
		Xtemp[i]=Xk[i];
	changeOrder(Xn,Xtemp,N);
	totalTime=(int)(log(N+1)/log(2));//The total time of butterfly
	wn=exp(-complex(0.0,1.0)*2.0*M_PI/(double)N);
	for (thisTime=1;thisTime<=totalTime;thisTime++){		//This time is now butterfly time	
			pointDistance=(1<<(thisTime-1));				//the distance between butterfly
			for (rowNo=0;rowNo<pointDistance;rowNo++){
				p=(1<<(totalTime-thisTime))*rowNo;
				for (lineNo=rowNo;lineNo<N;lineNo+=(1<<thisTime)){
					xtemp=Xn[lineNo];
					Xn[lineNo]=(Xn[lineNo]+Xn[lineNo+pointDistance]*(wn^((double)(-p))))*0.5;
					Xn[lineNo+pointDistance]=(xtemp-Xn[lineNo+pointDistance]*(wn^((double)(-p))))*0.5;
				}
			}
	}
	for (i=0;i<N;i++)
		xn[i]=Xn[i].real;

}


/*debug 
#include "elapsetime.h"
#include <conio.h>
void main(){
	double xn1[128],xn[128]={1.0,2.0,3.0,0.0};
	complex Xk[128];
	CycleCountStart();
	CFFT::FFT(xn,Xk,128);
	CycleCountElapse();
	int i;
	for (i=0;i<128;i++){
		Xk[i].print();
	}
	CFFT::IFFT(Xk,xn1,128);
	for (i=0;i<128;i++){
		cout<<xn1[i]<<endl;
	}
	getch();
	CycleCountStart();
	CDFT::DFT(xn,Xk,128);
	CycleCountElapse();
	for (i=0;i<128;i++){
		Xk[i].print();
	}
	CDFT::IDFT(Xk,xn1,128);
	for (i=0;i<128;i++){
		cout<<xn1[i]<<endl;
	}

}/*
*From the debug we can see when N=128 FFT used 5ms 
*while DFT used 80ms
*
*It tolds us FFT is more efficient than DFT
*Powered By Luda ^_^
//*/