2dfft.txt

本程序是触摸屏相关的驱动

#include <stdio.h>
#include<math.h>
#define M 64
#define N 256
#define dx 10     //单位是米//
#define dt 0.002  //单位是秒//
void fft(float *xr ,float *xi, int nr, float T);
int IBIT(int j, int m);

void main()
{
	int i,j;
	//float dk,df;
	float temp;
	float f[M][N],k[M][N];
	float x[M][N],xr[M][N],xi[M][N];
	float yr[N],yi[N],zr[M],zi[M]; 
	float A[M][N],B[M][N];
	FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6;
	fp1=fopen("data","rb");
    fp2=fopen("da.txt","w");  
	fp3=fopen("asp.txt","w"); //ap=amplitude spectrum//
	fp4=fopen("a8888","wb");	
	fp5=fopen("amp.txt","w");
	fp6=fopen("99","wb");


	//dk=1/(dx*M);
	//df=1/(dt*N);
//	printf("%f\t%f\n",dk,df);
	
	for(i=0;i<M;i++)
		for(j=0;j<N;j++)
		{
			fread(&x[i][j],4,1,fp1);
			fprintf(fp2,"%d\t%d\t%f\n",i,j,x[i][j]);
			//printf("%d\t%d\t%f\n",i,j,x[i][j]);
		}
		fclose(fp1);	
		fclose(fp2);
		
	for(i=0;i<M;i++)
		for(j=0;j<N;j++)
		{
			xr[i][j]=0;
			xi[i][j]=0;
		}
	
	for(i=0;i<M;i++)
		for(j=0;j<N;j++)
		{
			xr[i][j]=x[i][j];
			//printf("%d\t%d\t%f\n",i,j,xr[i][j]);
		}



//***********************进行二维傅立叶变换**************************//	
			for(i=0;i<M;i++)
			{
				for(j=0;j<N;j++)
				{
					yr[j]=xr[i][j];
					yi[j]=xi[i][j];
				}	
				
				fft(yr,yi,8,1.0);
				
				for(j=0;j<N;j++)
				{
					xr[i][j]=yr[j];
					xi[i][j]=yi[j];
					//fprintf(fp2,"%d\t%d\t%f %f\n",i,j,xr[i][j],xi[i][j]);
				}
			}	
			
			//fclose(fp2);
			
			for(j=0;j<N;j++)
			{
				for(i=0;i<M;i++)
				{
					zr[i]=xr[i][j];
					zi[i]=xi[i][j];
				}
				fft(zr,zi,6,1.0);
				for(i=0;i<M;i++)
				{
					xr[i][j]=zr[i];
					xi[i][j]=zi[i];
					//fprintf(fp2,"%d\t%d\t%f %f\n",i,j,xr[i][j],xi[i][j]);
				}
			}	
			
			
			//********************计算二维傅立叶变换后的振幅谱********************//
			for(i=0;i<M;i++)
				for(j=0;j<N;j++)
				{
					A[i][j]=sqrt(xr[i][j]*xr[i][j]+xi[i][j]*xi[i][j]);
					fprintf(fp3,"%d\t%d\t%f\n",i,j,A[i][j]);
					//fwrite(&i,sizeof(int),1,fp4);
                    //fwrite(&j,sizeof(int),1,fp4);
					fwrite(&A[i][j],sizeof(float),1,fp4);
				}
				fclose(fp3);
				fclose(fp4);
		
				
            for(i=0;i<M;i++)
					for(j=0;j<N;j++) 
					{
						k[i][j]=1.0/(dx*M)*i*10000;
					    f[i][j]=1.0/(dt*N)*j;
						A[i][j]=sqrt(xr[i][j]*xr[i][j]+xi[i][j]*xi[i][j]);	
						fprintf(fp5,"%f\t%f\t%f\n",k[i][j],f[i][j],A[i][j]);
					    //fwrite(&k[i][j],sizeof(float),1,fp4);
                        //fwrite(&f[i][j],sizeof(float),1,fp4);
					    fwrite(&A[i][j],sizeof(float),1,fp4);
					}
                    fclose(fp5);

//******************************象限颠倒************************************//
					for(i=0;i<M;i++)
						for(j=0;j<N;j++)
						{
							B[i][N-1-j]=A[i][j]; 
//					fwrite(&A[i][j],sizeof(float),1,fp6);
						}for(i=0;i<M;i++)
						for(j=0;j<N;j++)
						{fwrite(&B[i][j],sizeof(float),1,fp6);
						}
						fclose(fp6);
}



//------- 一维快速傅立叶(反)变换 -----------//
int IBIT(int j, int m)
{
	int i, it, j1, j2;
    it=0;
    j1=j;
    for(i=1; i<=m; i++)
	{
		j2=j1/2;
		it=it*2+(j1-2*j2);
		j1=j2;
	}
    return it;
}
void fft(float *xr ,float *xi, int nr, float T)
{
	int   i, j, k, l, n, n2, nr1, i1, j1, k2, k1n2;
	float c, s, p, tr, ti, trc, tic, ars;

    n=(int)(pow(2,nr));
    n2=n/2;
    nr1=nr-1;

    k=0;
    for(l=1; l<=nr; l++)
	{
loop1:  for(j=1; j<=n2; j++)
		{
	    k2=k/(int)(pow(2,nr1));
	    p=(double)(IBIT(k2,nr));
	    ars=(double)(6.2831852*p/(double)(n));
	    c=(double)(cos(ars));
	    s=(double)(sin(ars));
            k1n2=k+n2;
            tr=xr[k1n2]*c+xi[k1n2]*s*T;
            ti=xi[k1n2]*c-xr[k1n2]*s*T;
            xr[k1n2]=xr[k]-tr;
            xi[k1n2]=xi[k]-ti;
            xr[k]=xr[k]+tr;
            xi[k]=xi[k]+ti;
            k++;		
		}
		k+=n2;
        if(k<n)
		{  goto loop1;  }
		else
		{
			k=0;
            nr1=nr1-1;
            n2 /=2;
		}
	}
    for(j=1; j<=n; j++)
	{
		i=IBIT(j-1,nr)+1;
		if(i>j)
		{
			j1=j-1;
			i1=i-1;
			trc   =xr[j1];
            tic   =xi[j1];
            xr[j1]=xr[i1];
            xi[j1]=xi[i1];
            xr[i1]=trc;
            xi[i1]=tic;
		}
	}
    if(T<0.0)
    {
		for(j=0; j<n; j++)
		{
			xr[j]/=n;
		    xi[j]/=n;
		}
    }
}