fourier.c

tms320_6713上实现傅立叶变换的程序

#include<c6x.h>
#include<math.h>

#define IMAGEWIDTH 128
#define IMAGEHEIGHT 128
#define MODEBLOCK 6
#define MODEINCLINE 7
#define MODEFPHOTO 8

#define PI 3.14159265

typedef struct complex_struct
{
	float real;
	float img;
} complex;

void InitImage(unsigned int nMode,unsigned char *pImage,int nWidth,int nHeight);
void FFT(complex * TD,complex * FD,int r);
void Fourier(unsigned char *pImage,int nWidth,int nHeight);

unsigned char dbImage[IMAGEWIDTH*IMAGEHEIGHT];
unsigned char dbTargetImage[IMAGEWIDTH*IMAGEHEIGHT];
complex TD[IMAGEWIDTH*IMAGEHEIGHT];
complex FD[IMAGEWIDTH*IMAGEHEIGHT];
complex w[IMAGEWIDTH];
complex x1[IMAGEWIDTH];
complex x2[IMAGEWIDTH];

/* 傅立叶变换实验程序 */
int main()
{
	InitImage(MODEBLOCK,dbImage,IMAGEWIDTH,IMAGEHEIGHT);
	Fourier(dbImage,IMAGEWIDTH,IMAGEHEIGHT);
	InitImage(MODEINCLINE,dbImage,IMAGEWIDTH,IMAGEHEIGHT);	//BreakPoint
	Fourier(dbImage,IMAGEWIDTH,IMAGEHEIGHT);
	InitImage(MODEFPHOTO,dbImage,IMAGEWIDTH,IMAGEHEIGHT);	//BreakPoint
	Fourier(dbImage,IMAGEWIDTH,IMAGEHEIGHT);
	while ( 1 );											//BreakPoint
}		

void Fourier(unsigned char *pImage,int lWidth, int lHeight)
{
	// 中间变量
	float fTemp,fWork1,fWork2;
	unsigned char *pWork,*pWork3;
	complex *pWork1,*pWork2;

	// 循环变量
	int	i;
	int	j;
	
	// 进行付立叶变换的宽度和高度（2的整数次方）
	int	w;
	int	h;
	
	int	wp;
	int	hp;
	
	// 赋初值
	w = 1;
	h = 1;
	wp = 0;
	hp = 0;
	
	// 计算进行付立叶变换的宽度和高度（2的整数次方）
	while ( w*2<=lWidth )
	{
		w*=2; wp++;
	}
	while ( h*2<=lHeight )
	{
		h*=2; hp++;
	}
	pWork=pImage; pWork1=TD;
	for ( i=0;i<h;i++ )	// 行
		for ( j=0;j<w;j++,pWork++,pWork1++ )	// 列
		{
			pWork1->real=(*pWork);
			pWork1->img=0;
		}
	pWork1=TD; pWork2=FD;
	// 对y方向进行快速付立叶变换
	for ( i=0;i<h;i++ )	FFT(pWork1+w*i,pWork2+w*i,wp);
	// 保存变换结果
	for ( i=0;i<h;i++ )
		for ( j=0;j<w;j++ )
			TD[j*h+i]=FD[i*w+j];

	// 对x方向进行快速付立叶变换
	for ( i=0;i<w;i++ )	FFT(pWork1+i*h,pWork2+i*h,hp);
	for ( i=0;i<w;i++ )
		for ( j=0;j<h;j++ )
			TD[j*w+i]=FD[i*h+j];
	pWork1=TD; pWork2=FD; pWork=dbTargetImage;
	for ( i=0;i<h;i++ )	// 行
	{
		for ( j=0;j<w;j++ )	// 列
		{
			// 计算频谱
			fWork1=(pWork2+j*h+i)->real;
			fWork2=(pWork2+j*h+i)->img;
			fTemp=sqrt(fWork1*fWork1+fWork2*fWork2)/100;	//将幅度适当缩小后转换成灰度显示
			if ( fTemp>255 )	// 判断是否超过255
				fTemp = 255;	// 对于超过的，直接设置为255
			pWork3=pWork+lWidth*(lHeight-1-(i<h/2?i+h/2:i-h/2))+(j<w/2?j+w/2:j-w/2);
			(*pWork3)=(unsigned char)fTemp;
		}
	}
}

void FFT(complex *TD,complex *FD,int r)
{
	// 付立叶变换点数
	int	count;
	// 循环变量
	int i,j,k;
	// 中间变量
	int bfsize,p;
	// 角度
	float angle;
	complex *W,*X1,*X2,*X;
	
	W=w; X1=x1; X2=x2;
	// 计算付立叶变换点数
	count=1<<r;
	
	// 计算加权系数
	for ( i=0;i<count/2;i++ )
	{
		angle=-i*PI*2/count;
		(W+i)->real=cos(angle);
		(W+i)->img=sin(angle);	
	}
	
	// 将时域点写入X1
	for ( i=0;i<count;i++ )
	{
		(X1+i)->real=(TD+i)->real;
		(X1+i)->img=(TD+i)->img;
	}
	
	// 采用蝶形算法进行快速付立叶变换
	for ( k=0;k<r;k++ )
	{
		for ( j=0;j<1<<k;j++ )
		{
			bfsize=1<<(r-k);
			for ( i=0;i<bfsize/2;i++ )
			{
				p=j*bfsize;
				(X2+i+p)->real=(X1+i+p)->real+(X1+i+p+bfsize/2)->real;
				(X2+i+p)->img=(X1+i+p)->img+(X1+i+p+bfsize/2)->img;
				(X2+i+p+bfsize/2)->real=((X1+i+p)->real-(X1+i+p+bfsize/2)->real)*(W+i*(1<<k))->real 
				                          -((X1+i+p)->img-(X1+i+p+bfsize/2)->img)*(W+i*(1<<k))->img;
				(X2+i+p+bfsize/2)->img=((X1+i+p)->real-(X1+i+p+bfsize/2)->real)*(W+i*(1<<k))->img 
				                          +((X1+i+p)->img-(X1+i+p+bfsize/2)->img)*(W+i*(1<<k))->real;             		
			}
		}
		X=X1; X1=X2; X2=X;
	}
	
	// 重新排序
	for ( j=0;j<count;j++ )
	{
		p = 0;
		for ( i=0;i<r;i++ )
			if ( j&(1<<i) )
				p+=1<<(r-i-1);
		(FD+j)->real=(X1+p)->real;
		(FD+j)->img=(X1+p)->img;
	}
}