复件 (2) fft.c

适用于TI公司DM642的FFT例程 很好的资料

#include <stdlib.h>
#include <stdio.h>
#include <math.h>


// 常数π
#define	PI 			3.1415926535
#define DLP_Q0064 	0x80200004
#define PCR1 		*(volatile unsigned int  *)0x01900024

short 	fftdata[2048];
short 	fftdatan[2048];
float 	col[2048];
float 	sine[512] ;
float 	cosine[512];
#pragma 		DATA_SECTION(newbuf,"H263VCBUF");
unsigned char 	newbuf[512*512];
#pragma 		DATA_SECTION(tempbuf,"H263VCBUF");
unsigned char 	tempbuf[512*512];
#pragma 		DATA_SECTION(fftbuf,"H263VCBUF");
float 			fftbuf[512*512*2];


volatile unsigned char *ptr = (unsigned char *)0x80200004;
void initEMIF(void)
{
	#define EMIF_GCTL       0x01800000
	#define EMIF_CE1        0x01800004
	#define EMIF_CE0        0x01800008
	#define EMIF_CE2        0x01800010
	#define EMIF_CE3        0x01800014
	#define EMIF_SDRAMCTL   0x01800018
	#define EMIF_SDRAMTIMING  0x0180001C
  	#define EMIF_SDRAMEXT     0x01800020


  	*(int *)EMIF_GCTL = 0x00003040;/* EMIF global control register         */
  	*(int *)EMIF_CE1 = 0xffffff23; /* CE1 - 32-bit asynch access after boot*/
  	*(int *)EMIF_CE0 = 0xFFFFFF30; /* CE0 - SDRAM                          */
  	*(int *)EMIF_CE2 = 0x21200420; /* CE2 - 32-bit asynch on daughterboard */
  	*(int *)EMIF_CE3 = 0xffffff23; /* CE3 - 32-bit asynch on daughterboard */
  	*(int *)EMIF_SDRAMCTL = 0x57116000; /* SDRAM control register (100 MHz)*/
  	*(int *)EMIF_SDRAMTIMING = 0x00000410; /* SDRAM Timing register        */
}           
///////////////////////////////////////////
// Function name : BitReverse
// Description : 二进制倒序操作
// Return type : int 
// Argument : int src 待倒读的数
// Argument : int size 二进制位数
int BitReverse(int src, int size)
{
	int i;
	int tmp = src;
	int des = 0;
	for (i=size-1; i>=0; i--)
	{
		des = ((tmp & 0x1) << i) | des;
		tmp = tmp >> 1;
	}
	return des;
}

//////////////////////////////////////////////////
// Function name : Reorder
// Description : 数据二进制整序
// Return type : void 
// Argument :  x[MAX_N] Comp待整序数组
// Argument : int N FFT点数
// Argument : int M 点数的2的幂次
void Reorder(float x[], int N, int M)
{
	int i,j;
	for (i=0; i<N*2; i+=2)
	{
		j = BitReverse(i>>1, M)<<1;
		fftdatan[i] = x[j];
		fftdatan[i+1] = x[j+1];
	}	
	// 重新存入原数据中(已经是二进制整序过了的数据)
	for (i=0; i<N*2; i++)
	x[i] = fftdatan[i];
}

void FFT_1D(short m, float xy[])                            
{                                                                      
	
        float xt,yt,c,s;
        double angle;
   		short n1,n2,ie,ia,i,j,k,l;                                                
        short n =1<<m;                                                                      
        n2 = n;                                                              
        ie = 1;                    
        for(i = 0; i < n / 2; i++)
		{
			angle = i * PI * 2 / n;
			sine[i] =(float)sin(angle);
			cosine[i] =(float)cos(angle);
		}                                          
        for (k=n; k > 1; k = (k >> 1) ) {                                    
        	n1 = n2;                                                           
           	n2 = n2>>1;                                                        
           	ia = 0;                                                            
           	for (j=0; j < n2; j++) {                                           
             		c = -cosine[ia];                                                     
             		s = -sine[ia];                                                   
             		ia = ia + ie;                                                    
             		for (i=j; i < n; i += n1) {                                      
               			l = i + n2;                                                    
               			xt      = xy[2*l] - xy[2*i];                                   
               			xy[2*i] = xy[2*i] + xy[2*l];                                   
               			yt      = xy[2*l+1] - xy[2*i+1];                               
               			xy[2*i+1] = xy[2*i+1] + xy[2*l+1];                             
               			xy[2*l]   = (float)((c*xt + s*yt));                                 
               			xy[2*l+1] = (float)((c*yt - s*xt));                                 
             		}                                                                
           	}                                                                  
           	ie = ie<<1;                                                        
         }         
         Reorder(xy, n, m);                                                           
	return;
}                                                                      
void FFT_2D(float *x,int M)
{
	int i,j;
	int N = (1 << M);
	// 再逐列进行 1D-FFT
	for (j=0; j<N*2; j+=2)
	{
		// 取得第j列的数据
		for (i=0; i<N*2; i+=2)
		{
			col[i] = x[j + i*N];
			col[i+1] = x[j+1 + i*N];
		}	
		
		// 对第j列数据进行 1D-FFT
		FFT_1D(M,col); // <--- 计算结果在数组col中
		
		// 将结果放回矩阵第j列中
		for(i=0; i<N*2; i+=2)
		{
			x[j + i*N] = col[i];
			x[j + 1 + i*N] = col[i+1];
		}	
	}

	// 先逐行进行 1D-FFT
	for (i=0; i<N; i++)
	{
		for(j=0;j<N*2;j++)
		{
			col[j] = x[i*N*2 + j];
		}
		FFT_1D(M,col); // <--- 计算结果再存入矩阵x中
		for(j=0;j<N*2;j++)
		{
			x[i*N*2 + j] = col[j];
		}
	}	



}
void main()
{
  	int   	i;
  	int  row,col;
  	unsigned char *ptr1;
  	unsigned int j ;
  	
//  	short   w[10]={-32767,0,-23170,-23170,0,-32767,23170,-23170,32767,0};
	initEMIF();
#if 0	
  	for(i=0;i<1024;i++)
    		fftbuf[i] = 0.0;
  	for(i=0;i<4;i++)
    		fftbuf[i*2] = 255;
    	FFT_1D(9,fftbuf);
	
    asm("\tnop");
    asm("\tnop");
    asm("\tnop");
    asm("\tnop");
    asm("\tnop");
    

	for(ptr1 =(unsigned char *)ptr,i=0;i<512*512;i++)   
	{
		ptr1[i] = 0;
	}
	for(ptr1=(unsigned char *)ptr,j=128;j<384;j++)
	{
		for(i=128;i<384;i++)
		{
			ptr1[j*512+i] = 0xFF;
		}
	}	
	asm("\tnop");
    asm("\tnop");
    asm("\tnop");
    asm("\tnop");
    asm("\tnop");
    for(i=0;i<512*512;i++)
	{
		((unsigned char *)tempbuf)[i] = 0x80;
		((unsigned char *)newbuf)[i] = 0x00;
	}
	for(i=0;i<512*512;i++)
	{
		fftbuf[i*2] = (float)(*ptr++);
		fftbuf[i*2+1] = 0;
	}	
	FFT_2D(fftbuf,9);
#endif	
	for(i = 0,row = 0;row<512;row++)
	{
		for(col = 0; col<512;col++,i++)
		{
			j =(unsigned int)sqrt(fftbuf[i*2]*fftbuf[i*2]+fftbuf[i*2+1]*fftbuf[i*2+1]);
			if(j>255) j=255;
			newbuf[((256+col)&0x1ff) +((256+row)&0x1ff)*512] = (unsigned char)j;		
		}
	}
	
	asm("\tnop");
	asm("\tnop");
	asm("\tnop");

	asm("\tnop");
	asm("\tnop");
	return;
}