一维快速傅里叶变换与反变换

/*************************************************************************  *  * \函数名称：  *   FFT_1D()  *  * \输入参数:  *   complex<double> * pCTData	- 指向时域数据的指针，输入的需要变换的数据  *   complex<double> * pCFData	- 指向频域数据的指针，输出的经过变换的数据  *   nLevel						－傅立叶变换蝶形算法的级数，2的幂数，  *  * \返回值:  *   无  *  * \说明:  *   一维快速傅立叶变换。  *  *************************************************************************  */  VOID WINAPI FFT_1D(complex<double> * pCTData, complex<double> * pCFData, int nLevel) {	 	// 循环控制变量 	int		i; 	int     j; 	int     k;  	double PI = 3.1415926;   	// 傅立叶变换点数 	int	nCount =0 ;  	// 计算傅立叶变换点数 	nCount =(int)pow(2,nLevel) ; 	 	// 某一级的长度 	int		nBtFlyLen; 	nBtFlyLen = 0 ; 	 	// 变换系数的角度 ＝2 * PI * i / nCount 	double	dAngle; 	 	complex<double> *pCW ; 	 	// 分配内存，存储傅立叶变化需要的系数表 	pCW  = new complex<double>[nCount / 2];      // 计算傅立叶变换的系数 	for(i = 0; i < nCount / 2; i++) 	{ 		dAngle = -2 * PI * i / nCount; 		pCW[i] = complex<double> ( cos(dAngle), sin(dAngle) ); 	}  	// 变换需要的工作空间 	complex<double> *pCWork1,*pCWork2;  	 	// 分配工作空间 	pCWork1 = new complex<double>[nCount];  	pCWork2 = new complex<double>[nCount];  	 	// 临时变量 	complex<double> *pCTmp; 	 	// 初始化，写入数据 	memcpy(pCWork1, pCTData, sizeof(complex<double>) * nCount);  	// 临时变量 	int nInter;  	nInter = 0;  	// 蝶形算法进行快速傅立叶变换 	for(k = 0; k < nLevel; k++) 	{ 		for(j = 0; j < (int)pow(2,k); j++) 		{ 			//计算长度 			nBtFlyLen = (int)pow( 2,(nLevel-k) ); 			 			//倒序重排，加权计算 			for(i = 0; i < nBtFlyLen/2; i++) 			{ 				nInter = j * nBtFlyLen; 				pCWork2[i + nInter] =  					pCWork1[i + nInter] + pCWork1[i + nInter + nBtFlyLen / 2]; 				pCWork2[i + nInter + nBtFlyLen / 2] = 					(pCWork1[i + nInter] - pCWork1[i + nInter + nBtFlyLen / 2])  					* pCW[(int)(i * pow(2,k))]; 			} 		}  		// 交换 pCWork1和pCWork2的数据 		pCTmp   = pCWork1	; 		pCWork1 = pCWork2	; 		pCWork2 = pCTmp		; 	} 	 	// 重新排序 	for(j = 0; j < nCount; j++) 	{ 		nInter = 0; 		for(i = 0; i < nLevel; i++) 		{ 			if ( j&(1<<i) ) 			{ 				nInter += 1<<(nLevel-i-1); 			} 		} 		pCFData[j]=pCWork1[nInter]; 	} 	 	// 释放内存空间 	delete pCW; 	delete pCWork1; 	delete pCWork2; 	pCW		=	NULL	; 	pCWork1 =	NULL	; 	pCWork2 =	NULL	;  }  /*************************************************************************  *  * \函数名称：  *    IFFT_1D()  *  * \输入参数:  *   complex<double> * pCTData	- 指向时域数据的指针，输入的需要反变换的数据  *   complex<double> * pCFData	- 指向频域数据的指针，输出的经过反变换的数据  *   nLevel						－傅立叶变换蝶形算法的级数，2的幂数，  *  * \返回值:  *   无  *  * \说明:  *   一维快速傅立叶反变换。  *  ************************************************************************  */ VOID WINAPI IFFT_1D(complex<double> * pCFData, complex<double> * pCTData, int nLevel) { 	 	// 循环控制变量 	int		i;  	// 傅立叶反变换点数 	int nCount;  	// 计算傅立叶变换点数 	nCount = (int)pow(2,nLevel) ; 	 	// 变换需要的工作空间 	complex<double> *pCWork;	 	 	// 分配工作空间 	pCWork = new complex<double>[nCount]; 	 	// 将需要反变换的数据写入工作空间pCWork 	memcpy(pCWork, pCFData, sizeof(complex<double>) * nCount); 	 	// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭 	// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭 	for(i = 0; i < nCount; i++) 	{ 		pCWork[i] = complex<double> (pCWork[i].real(), -pCWork[i].imag()); 	} 	 	// 调用快速傅立叶变换实现反变换，结果存储在pCTData中 	FFT_1D(pCWork, pCTData, nLevel); 	 	// 求时域点的共轭，求得最终结果 	// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据 	// 相差一个系数nCount 	for(i = 0; i < nCount; i++) 	{ 		pCTData[i] =  			complex<double> (pCTData[i].real() / nCount, -pCTData[i].imag() / nCount); 	} 	 	// 释放内存 	delete pCWork; 	pCWork = NULL; }
二维快速傅里叶变换与反变换 

/*************************************************************************  *  * \函数名称：  *   FFT_2D()  *  * \输入参数:  *   complex<double> * pCTData	- 图像数据  *   int    nWidth				- 数据宽度  *   int    nHeight				- 数据高度  *   complex<double> * pCFData	- 傅立叶变换后的结果  *  * \返回值:  *   无  *  * \说明:  *   二维快速傅立叶变换。  *  ************************************************************************  */ VOID WINAPI DIBFFT_2D(complex<double> * pCTData, int nWidth, int nHeight, complex<double> * pCFData) { 	 	// 循环控制变量 	int	x; 	int	y; 	 	// 临时变量 	double	dTmpOne; 	double  dTmpTwo; 	 	// 进行傅立叶变换的宽度和高度，（2的整数次幂） 	// 图像的宽度和高度不一定为2的整数次幂 	int		nTransWidth; 	int 	nTransHeight;  	// 计算进行傅立叶变换的宽度	（2的整数次幂） 	dTmpOne = log(nWidth)/log(2); 	dTmpTwo = ceil(dTmpOne)		   ; 	dTmpTwo = pow(2,dTmpTwo)	   ; 	nTransWidth = (int) dTmpTwo	   ; 	 	// 计算进行傅立叶变换的高度 （2的整数次幂） 	dTmpOne = log(nHeight)/log(2); 	dTmpTwo = ceil(dTmpOne)		   ; 	dTmpTwo = pow(2,dTmpTwo)	   ; 	nTransHeight = (int) dTmpTwo	   ;	 	 	// x，y（行列）方向上的迭代次数 	int		nXLev; 	int		nYLev;  	// 计算x，y（行列）方向上的迭代次数 	nXLev = (int) ( log(nTransWidth)/log(2) +  0.5 ); 	nYLev = (int) ( log(nTransHeight)/log(2) + 0.5 ); 	 	for(y = 0; y < nTransHeight; y++) 	{ 		// x方向进行快速傅立叶变换 		FFT_1D(&pCTData[nTransWidth * y], &pCFData[nTransWidth * y], nXLev); 	} 	 	// pCFData中目前存储了pCTData经过行变换的结果 	// 为了直接利用FFT_1D，需要把pCFData的二维数据转置，再一次利用FFT_1D进行 	// 傅立叶行变换（实际上相当于对列进行傅立叶变换） 	for(y = 0; y < nTransHeight; y++) 	{ 		for(x = 0; x < nTransWidth; x++) 		{ 			pCTData[nTransHeight * x + y] = pCFData[nTransWidth * y + x]; 		} 	} 	 	for(x = 0; x < nTransWidth; x++) 	{ 		// 对x方向进行快速傅立叶变换，实际上相当于对原来的图象数据进行列方向的 		// 傅立叶变换 		FFT_1D(&pCTData[x * nTransHeight], &pCFData[x * nTransHeight], nYLev); 	}  	// pCFData中目前存储了pCTData经过二维傅立叶变换的结果，但是为了方便列方向 	// 的傅立叶变换，对其进行了转置，现在把结果转置回来 	for(y = 0; y < nTransHeight; y++) 	{ 		for(x = 0; x < nTransWidth; x++) 		{ 			pCTData[nTransWidth * y + x] = pCFData[nTransHeight * x + y]; 		} 	}  	memcpy(pCTData, pCFData, sizeof(complex<double>) * nTransHeight * nTransWidth ); }  /*************************************************************************  *  * \函数名称：  *   IFFT_2D()  *  * \输入参数:  *   complex<double> * pCFData	- 频域数据  *   complex<double> * pCTData	- 时域数据  *   int    nWidth				- 图像数据宽度  *   int    nHeight				- 图像数据高度  *  * \返回值:  *   无  *  * \说明:  *   二维傅立叶反变换。  *  ************************************************************************  */ VOID WINAPI IFFT_2D(complex<double> * pCFData, complex<double> * pCTData, int nWidth, int nHeight) { 	// 循环控制变量 	int	x; 	int	y; 	 	// 临时变量 	double	dTmpOne; 	double  dTmpTwo; 	 	// 进行傅立叶变换的宽度和高度，（2的整数次幂） 	// 图像的宽度和高度不一定为2的整数次幂 	int		nTransWidth; 	int 	nTransHeight;  	// 计算进行傅立叶变换的宽度	（2的整数次幂） 	dTmpOne = log(nWidth)/log(2); 	dTmpTwo = ceil(dTmpOne)		   ; 	dTmpTwo = pow(2,dTmpTwo)	   ; 	nTransWidth = (int) dTmpTwo	   ; 	 	// 计算进行傅立叶变换的高度 （2的整数次幂） 	dTmpOne = log(nHeight)/log(2); 	dTmpTwo = ceil(dTmpOne)		   ; 	dTmpTwo = pow(2,dTmpTwo)	   ; 	nTransHeight = (int) dTmpTwo	   ; 	 	// 分配工作需要的内存空间 	complex<double> *pCWork= new complex<double>[nTransWidth * nTransHeight];  	//临时变量 	complex<double> *pCTmp ; 	 	// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭 	// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭 	for(y = 0; y < nTransHeight; y++) 	{ 		for(x = 0; x < nTransWidth; x++) 		{ 			pCTmp = &pCFData[nTransWidth * y + x] ; 			pCWork[nTransWidth * y + x] = complex<double>( pCTmp->real() , -pCTmp->imag() ); 		} 	}  	// 调用傅立叶正变换 	::DIBFFT_2D(pCWork, nWidth, nHeight, pCTData) ; 	 	// 求时域点的共轭，求得最终结果 	// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据 	// 相差一个系数 	for(y = 0; y < nTransHeight; y++) 	{ 		for(x = 0; x < nTransWidth; x++) 		{ 			pCTmp = &pCTData[nTransWidth * y + x] ; 			pCTData[nTransWidth * y + x] =  				complex<double>( pCTmp->real()/(nTransWidth*nTransHeight), 								 -pCTmp->imag()/(nTransWidth*nTransHeight) ); 		} 	} 	delete pCWork ; 	pCWork = NULL ; }