DSP编程

开发平台：
C/C++

ICETEK-DM642-C.c：源码内容
							#include "math.h"
#include <csl_cache.h>
#include "ICETEK-DM642-C.h"
//工作变量
unsigned char m_dbFrameBuffer[SIMGWIDTH*SIMGHEIGHT];
unsigned char m_dbTargetImage[TIMGWIDTH*TIMGHEIGHT];
unsigned char m_dbFFTImage[FFTIMGWIDTH*FFTIMGHEIGHT];
unsigned int m_uVideoStatus,m_bFreeze;
complex TD[TIMGWIDTH*TIMGHEIGHT];
complex FD[TIMGWIDTH*TIMGHEIGHT];
complex w[TIMGWIDTH];
complex x1[TIMGWIDTH];
complex x2[TIMGWIDTH];
void ICETEKDM642CInit()
{
	m_uVideoStatus=m_bFreeze=0;
}
void ICETEKDM642CFourier()
{
	// 中间变量
	float fTemp,fWork1,fWork2;
	unsigned char *pWork,*pWork3;
	int nWork,ii,jj;
	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<=TIMGWIDTH )
	{
		w*=2; wp++;
	}
	while ( h*2<=TIMGHEIGHT )
	{
		h*=2; hp++;
	}
	pWork1=TD;
	for ( i=0;i<h;i++ )	// 行
	{
		pWork=m_dbFrameBuffer+i*SIMGWIDTH;
		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=m_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+TIMGWIDTH*(TIMGHEIGHT-1-(i<h/2?i+h/2:i-h/2))+(j<w/2?j+w/2:j-w/2);
			(*pWork3)=(unsigned char)fTemp;
		}
	}
	pWork3=m_dbFFTImage;
	for ( i=0;i<FFTIMGHEIGHT;i++ )
	{
		for ( j=0;j<FFTIMGWIDTH;j++,pWork3++ )
		{
			pWork=m_dbTargetImage+TIMGHEIGHT/FFTIMGHEIGHT*i*TIMGHEIGHT+TIMGWIDTH/FFTIMGWIDTH*j;
			nWork=0;
			for ( ii=0;ii<TIMGHEIGHT/FFTIMGHEIGHT;ii++,pWork+=TIMGWIDTH )
				for ( jj=0;jj<TIMGWIDTH/FFTIMGWIDTH;jj++ )
					nWork+=pWork[jj];
			nWork/=((TIMGHEIGHT/FFTIMGHEIGHT)*(TIMGWIDTH/FFTIMGWIDTH));
			(*pWork3)=nWork;
		}
	}
	CACHE_clean(CACHE_L2ALL, 0, 0);
}
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;
	}
	CACHE_clean(CACHE_L2ALL, 0, 0);
}

						