fourier.cpp

一般的论文 毕业用的 有很大的用处

//1 首先定义一个复数的结构

struct CplexNum
{
	double re;
	double im;
};

//2 申请两块内存，pTd存放输入数据，pFd存放Fourier变换的数据
CplexNum *pTd = new CplexNum[sizeof(CplexNum)*width * height];// 分配内存
CplexNum *pFd = new CplexNum[sizeof(CplexNum)*width * height];

//3 将pTd数组的实部赋以图像数据的值，调用下面的函数进行处理后得到pFd的值，
// 不一定大家一定要遵照函数的参数要求，可以根据需要修改

//4 得到pFd的数据后，根据求每个复数的模，并将该模赋给相应的象素，不过要小心数据超界，所以
//最好是归一后在乘以相应的图像最大灰度级。

//用于复数运算
CplexNum Add(CplexNum c1,CplexNum c2)
{
	CplexNum c;
	c.re=c1.re+c2.re;
	c.im=c1.im+c2.im;
	return c;
}
CplexNum Sub(CplexNum c1,CplexNum c2)
{
	CplexNum c;
	c.re=c1.re-c2.re;
	c.im=c1.im-c2.im;
	return c;
}
CplexNum Mul(CplexNum c1,CplexNum c2)
{
	CplexNum c;
	c.re=c1.re*c2.re-c1.im*c2.im;
	c.im=c1.re*c2.im+c2.re*c1.im;
	return c;
}

/*************************************************************************
 * 函数名称：Fourier(CplexNum * pTd, int lWidth, int lHeight, CplexNum * pFd)
 * 函数参数:
 *   CplexNum * pTd，指向时域值的指针
 *   int    lWidth，图象宽度
 *   int    lHeight，图象高度
 *   CplexNum * pFd	，指向频域值的指针
 *
 * 函数类型:void
 * 函数功能:二维快速傅立叶变换
 *************************************************************************/
void  Fourier(CplexNum * pTd, int lWidth, int lHeight, CplexNum * pFd)
{
	
	// 循环控制变量
	int	j;
	int	i;
	// 进行傅立叶变换的宽度和高度，（2的整数次幂）
	// 图象的宽度和高度不一定为2的整数次幂
	int		wid=1;
	int 	hei=1;
	int		widpor=0,heiPor=0;//2的幂数

	while(wid * 2 <= lWidth)// 计算进行付立叶变换的宽度和高度（2的整数次方）
	{
		wid *= 2;
		widpor++;
	}	
	while(hei * 2 <= lHeight)
	{
		hei *= 2;
		heiPor++;
	}	
	
	for(i = 0; i < hei; i++)
	{
		// x方向进行快速傅立叶变换
		FastFourierTran(&pTd[wid * i], &pFd[wid * i], widpor);
	}
	
	// pFd中目前存储了pTd经过行变换的结果
	// 为了直接利用FastFourierTran，需要把pFd的二维数据转置，再一次利用FastFourierTran进行
	// 傅立叶行变换（实际上相当于对列进行傅立叶变换）
	for(i = 0; i < hei; i++)
	{
		for(j = 0; j < wid; j++)
		{
			pTd[hei * j + i] = pFd[wid * i + j];
		}
	}
	
	for(j = 0; j < wid; j++)
	{
		// 对x方向进行快速傅立叶变换，实际上相当于对原来的图象数据进行列方向的
		// 傅立叶变换
		FastFourierTran(&pTd[j * hei], &pFd[j * hei], heiPor);
	}

	// pFd中目前存储了pTd经过二维傅立叶变换的结果，但是为了方便列方向
	// 的傅立叶变换，对其进行了转置，现在把结果转置回来
	for(i = 0; i < hei; i++)
	{
		for(j = 0; j < wid; j++)
		{
			pTd[wid * i + j] = pFd[hei * j + i];
		}
	}

	memcpy(pTd, pFd, sizeof(CplexNum) * hei * wid );
}





/*************************************************************************
 * 函数名称：FastFourierTran(CplexNum * pTd, CplexNum* pFd, int power)
 * 函数参数:
 *   CplexNum * pTd,指向时域数组的指针
 *   CplexNum * pFd,指向频域数组的指针
 *   int             power,2的幂数，即迭代次数
 * 函数类型:void 
 函数功能:用来实现快速付立叶变换
************************************************************************/
void  FastFourierTran(CplexNum * pTd, CplexNum * pFd, int power)
{	
	long i;                 //行循环变量
	long j;                 //列循环变量
 			
	long	dotCount;// 付立叶变换点数		
	int		k;// 循环变量		
	int		bfsize,p;// 中间变量		
	double	angle;// 角度	
	CplexNum *pWn,*temReg1,*temReg2,*temReg;	
	
	dotCount= 1 <<power;// 计算付立叶变换点数		
	pWn= new CplexNum[sizeof(CplexNum)*dotCount/ 2];// 分配运算所需存储器
	temReg1 = new CplexNum[sizeof(CplexNum)*dotCount];
	temReg2 = new CplexNum[sizeof(CplexNum)*dotCount];		
	for(i = 0; i < dotCount/ 2; i++)// 计算加权系数
	{
		angle = -i * pi* 2 / dotCount;
		pWn[i].re = cos(angle);
		pWn[i].im=sin(angle);
	}	
	memcpy(temReg1, pTd, sizeof(CplexNum)*dotCount);// 将时域点写入temReg1		
	for(k = 0; k < power; k++)// 采用蝶形算法进行快速付立叶变换
	{
		for(j = 0; j < 1 << k; j++)
		{
			bfsize = 1 << (power-k);
			for(i = 0; i < bfsize / 2; i++)
			{
				p = j * bfsize;
				temReg2[i+p]=Add(temReg1[i+p],temReg1[i+p+bfsize/2]);
				temReg2[i+p+bfsize/2]=Mul(Sub(temReg1[i+p],temReg1[i+p+bfsize/2]),
				pWn[i*(1<<k)]);
			}
		}
		temReg  = temReg1;
		temReg1 = temReg2;
		temReg2 = temReg;
	}		
	for(j = 0; j <dotCount; j++)// 重新排序
	{
		p = 0;
		for(i = 0; i <power; i++)
		{
			if (j&(1<<i))
			{
				p+=1<<(power-i-1);
			}
		}
		pFd[j]=temReg1[p];
	}		
	delete pWn;// 释放内存
	delete temReg1;
	delete temReg2;
}