C语言实现傅里叶变换函数dft,idft,fft,ifft

自定义结构体complex（复数）

typedef struct{
	double re; //实部 
	double im; //虚部 
}complex;

用于为离散序列倒码排序的函数

int  reverse(int n,int M){
	int i,m,dev,temp=0;
	for(i=0,m=1,dev=1;i<M;i++,m*=2,dev*=2){
		temp=(temp<<1)+(n&m)/dev;
	}
	return temp;
}

void reSequence(complex*src, complex*dst,int N){
	int i,M=(int)log2(N);
	for(i=0;i<N;i++){
		dst[i].re=src[reverse(i,M)].re;
		dst[i].im=src[reverse(i,M)].im;
	}
}

dft、idft、fft、ifft实现代码

void dft(complex*src,complex*dst,int N){
	int i,j; 
	double temp;
	for(i=0;i<N;i++){
		dst[i].re=0;
		dst[i].im=0;
		for(j=0;j<N;j++){
			temp=-2*pi*j*i/N;
			dst[i].re+=src[j].re*cos(temp)-src[j].im*sin(temp);
			dst[i].im+=src[j].re*sin(temp)+src[j].im*cos(temp);
		}
	}
}

void idft(complex*src,complex*dst,int N){
	int i,j;
	double temp;
	for(i=0;i<N;i++){
		dst[i].re=0;
		dst[i].im=0;
		for(j=0;j<N;j++){
			temp=2*pi*i*j/N;
			dst[i].re+=(src[j].re*cos(temp)-src[j].im*sin(temp))/N;
			dst[i].im+=(src[j].re*sin(temp)+src[j].im*cos(temp))/N;
		}
	}
}

void fft(complex*src,complex*dst,int N){
	int i,j;
	double temp;
	int g,n,gap;
	complex x1,x2;
	reSequence(src,dst,N);
	
	for(g=N/2;g>0;g/=2){	//g为蝴蝶单元群的数量 
		n=N/g;                 //n为每一蝴蝶单元群中蝴蝶单元的数量 
		gap=n/2;			//gap为蝴蝶单元上两个节点的间距 
		for(i=0;i<N;i+=n){
			for(j=i;j<i+n/2;j++){
				temp=-2*pi*g*(j-i)/double(N);
				x1.re=dst[j].re+dst[j+gap].re*cos(temp)-dst[j+gap].im*sin(temp);
				x1.im=dst[j].im+dst[j+gap].re*sin(temp)+dst[j+gap].im*cos(temp);
				x2.re=dst[j].re-dst[j+gap].re*cos(temp)+dst[j+gap].im*sin(temp);
				x2.im=dst[j].im-dst[j+gap].re*sin(temp)-dst[j+gap].im*cos(temp);
				dst[j].re=x1.re;dst[j].im=x1.im;dst[j+gap].re=x2.re;dst[j+gap].im=x2.im;
			}
		}
	}
}

void ifft(complex*src,complex*dst,int N){
	int i,j;
	double temp;
	int g,n,gap;
	complex x1,x2;
	reSequence(src,dst,N);
	
	for(g=N/2;g>0;g/=2){	//g为蝴蝶单元群的数量 
		n=N/g;                 //n为每一蝴蝶单元群中蝴蝶单元的数量 
		gap=n/2;			//gap为蝴蝶单元上两个节点的间距 
		for(i=0;i<N;i+=n){
			for(j=i;j<i+n/2;j++){
				temp=2*pi*g*(j-i)/double(N);
				x1.re=(dst[j].re+dst[j+gap].re*cos(temp)-dst[j+gap].im*sin(temp));
				x1.im=(dst[j].im+dst[j+gap].re*sin(temp)+dst[j+gap].im*cos(temp));
				x2.re=(dst[j].re-dst[j+gap].re*cos(temp)+dst[j+gap].im*sin(temp));
				x2.im=(dst[j].im-dst[j+gap].re*sin(temp)-dst[j+gap].im*cos(temp));
				dst[j].re=x1.re;dst[j].im=x1.im;dst[j+gap].re=x2.re;dst[j+gap].im=x2.im;
			}
		}
	}
	
	for(i=0;i<N;i++){
		dst[i].re/=N;dst[i].im/=N;
	}
}