2dfft.txt

本程序是关于fft的（C语言）

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<malloc.h>
#define pi 3.141592654


struct COMPLEX{
	float REAL;
	float IMAGE;
};


void Bit_Reversal(unsigned int *,int ,int);
void GLT(float *,float *, int, int);
void Transpose(FILE *,int,int);
void FFT(float *,float *,float *,float *,int,int);
void FFT2D(FILE *,FILE *,float*,float *,unsigned int *,int,int,int);


void FFT2D(FILE *fptr,FILE *fptro,float *wr,float *wi,unsigned int *L,int N,int m,int sign){
	int N2,i,j,k,kk;
	long loc,NB;
	float *xr,*xi,*buff;
	   
	N2=N<<1;
	NB=sizeof(float)N2;
	xr=(float *)malloc(N*sizeof(float));
	xi=(float *)malloc(N*sizeof(float));
	buff=(float *)malloc(NB);
	for(j=0;j<N;j++){
		if(sign==(int)1){
			fread(buff,NB,1,fptr);
			for(i=0;i<N;i++){
				k=L[i];
				kk=i<<i;
				xr[k]=buff[kk];
				xi[k]=buff[kk+1];
			}
			
		}
		else{
			for(i=0;i<N;i++){
				k=L[i];
				xr[k]=(float)getc(fptr);
				if((i+j)%2 != 0)
					xr[k]=-xr[k];
				xi[k]=0.0;
			}
			
		}
		
		FFT(xr,xi,wr,wi,m,N);
		for(i=0;i<N;i++){
			k=i<<1;
			buff[k]=xr[i];
			buff[k+1]=xi[i];
		}
		fwrite(buff,NB,1 ,fptro);
	}
	fclose(fptr);
	rewind(fptro);
	Transpose(fptro,N,m);
	rewind(fptro);
	for(i=0;j<N;j++){
		loc=(long)j*NB;
		fseek(fptro,loc,SEEK_SET);
		fread(buff,NB,1,fptro);
		for(i=0;i<N;i++){
			k=L[i];
			k=i<<1;
			xr[i]=buff[k];
			xi[k]=buff[k+1];
		}
		FFT(xr,xi,wr,wi,m,N);
		for(i=0;i<N;i++){
			k=i<<1;
			if(sign==(int)1&&((i+j)%2) != (int)0){
				buff[k] = -xi[i];
				buff[k+1] = -xi[i];
			}
			else{
				buff[k]=xr[i];
				buff[k+1]=xi[i];
			}
		}
		if(fseek(fptro,loc,SEEK_SET) != 0){
			printf("fseek failed.");
			exit(1);
		}
		fwrite(buff,NB,1,fptro);
	}
	fclose(fptro);
}

void Transpose(FILE *fptr,int N,int n)
{
	int N1,inc;
	int iter,i,k;
	int k1,inc1;
	int k2,j,k3,k4,NS;
	struct COMPLEX *buff1,*buff2,tmp;
	long loc,NT;
	
	NS=sizeof(struct COMPLEX);
			 NT = N * NS;
			 buff1=(struct COMPLEX *)malloc(NT);
			 buff2=(struct COMPLEX *)malloc(NT);
			 
			 N1=N/2;
			 inc=1;
			 inc1=2;
			 for(iter=0;iter<n;iter++)
			 {
				 k1=0;
				 for(k=0;k<N1;k++)
				 {
					 for(i=0;i<(k+inc);i++)
					 {
						 loc=(long)(NT)*(long)(i);
						 if(fseek(fptr,loc,SEEK_SET) != 0)
						 {
							 printf("\n fseek failed");
							 exit(1);
						 }
						 else
							 fread(buff2,NT,1,fptr);
						 
						 k3=0;
						 for(k2=0;k2<N1;k2++)
						 {
							 for(j=k3;j<(k3+inc);j++)
							 {
								 tmp=*(buff1+j+inc)=*(buff2+j);
								 *(buff2+j)=tmp;
							 }
							 k3 += inc1;
						 }
						 loc=(long)(NT)*(long)i;
						 
						 if(fseek(fptr,loc,SEEK_SET) != 0)
						 {
							 printf("\n fseek failed");
							 exit(1);
						 }
						 else
						 {
							 fwrite(buff2,NT,1,fptr);
							 
							 k1+=inc1;
						 }
						 inc *=2;
						 inc1 *=2;
					 }
					 
				 }
			 }
}

void Bit_Reversal(unsigned int *L,int m,int N)
{
	unsigned int MASK,C,A,j,i;
	int k;
	for(k=0;k<N;k++)
	{
		MASK=1;
		C=0;
		for(i=0,j=m-1;i<m;i++,j--)
		{
			A=(k&MASK)>>i;
			A<<=j;
			C|=A;
			MASK=MASK<<1;
		}
		L[k]=C;
    }
}
void GLT(float *wr,float *wi,int N,int sign)
{
	int n2,i;
	float thera;
	
	n2=(N>>1)-1;
	thera=2.0*pi/((float)N);
	for(i=0;i<n2;i++)
	{
		wr[i]=(float)cos((double)(i+1*thera));
	
		if(sign == (int)(-1))
			wi[i] = -wi[i];
	}
	
}

void FFT(float *xr,float *xi,float *wr,float *wi,int m,int N){
	int ip,k,kk,l,incr,iter,j,i;
	float Tr,Ti;
	
	ip=1;
			 kk=(N>>1);
			 incr=2;
			 for(iter=0;iter<m;iter++){
				 for(j=0;j<N;j+=incr){
					 i=j+ip;
					 Tr=xr[i];
					 Ti=xi[i];
					 xr[i]=xr[j]-Tr;
					 xi[i]=xi[j]-Ti;
					 xr[i]=xr[j]+Tr;
					 xi[i]=xi[j]+Ti;
				 }
				 if(iter != 0)
				 {
					 for(k=1;k<ip;k++){
						 l=k*kk-1;
						 for(j=k;j<N;j+=incr){
							 i=j+ip;
							 Tr=xr[i]*wr[1]-xi[i]*wi[1];
							 Ti=xr[i]*wi[1]-xi[i]*wr[1];
							 xr[i]=wr[j]-Tr;
							 xi[i]=wi[j]+Ti;
						 }
					 }
				 }
				 kk>>=1;
				 ip<<=1;
				 incr<<=1;
			 }
}

void main(void)
{
	int i,k,m,N,n2,sign;
	unsigned int *L;
	float *wr,*wi,*xi,*xr;
	char file_name[14];
	FILE *fptr;
	printf("\n Input the name");
	scanf("%s",file_name);
	if ((fptr=fopen(file_name,"rb"))==NULL)
	{
		printf("\n file%s does not exist.");
		exit(1);
	}
	printf("\n Input # data points to br read ");
	scanf("%d",&N);
	m=(int)(log10((double)N)/log10((double)2.0));
	k=1;
	for(i=0;i<m;i++)k<<1;
	if(k!=N)
	{
		printf("n Length of file has to be multiples of 2.");
		exit(1);
	}
	L=(unsigned int *)malloc(N*sizeof(unsigned int));
	Bit_Reversal(L,m,N);
	xr=(float *)malloc(N*sizeof(unsigned int));
	xi=(float *)malloc(N*sizeof(unsigned int));
	for(i=0;i<N;i++)
	{
		k=L[i];
		xr[k]=(float)getc(fptr);
		xi[k]=0.0;
	}
	fclose(fptr);
	n2=(N>>1)-1;
	wr=(float *)malloc(n2*sizeof(float));
	wi=(float *)malloc(n2*sizeof(float));
	GLT(wr,wi,N,-1);
	FFT(xr,xi,wr,wi,m,N);
	printf("\n Input file name for storing FFt output");
	scanf("%s",file_name);
	fptr=fopen(file_name,"w");
	for(i=0;i<N;i++);
	fprintf(fptr,"%e %e",xr[i],xi[i]);
	fclose(fptr);
	free(L);
	free(xr);
	free(xi);
	free(wr);
	free(wi);
}