一维快速傅里叶变换与反变换
/************************************************************************* * * \函数名称： * 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 ; } 
