快速傅里叶变换与逆变换——CFFT2D类的实现 前言在快速傅里叶变换（FFT）未出之前，傅里叶变换并未大为流行，主要的原因就是变换需要的计算时间过多。在FFT出现后，傅里叶变换才算是真正的运用到了实际生产生活当中。所以，可以这样说是FFT推广了傅里叶变换。FFT变换是频率域图像处理的基石，所以对于数字图像处理学**来说，学**FFT算法是非常必须的。虽然，现在Internet有很多的商业、开源的代码库提供了FFT算法实现，更有甚者，在DSP的配套硬件中，有专门的FFT实现硬件，运行速度得到更高的提升。但是，对于学**来说确实不利的，就好像一个“黑盒子”一样，我把图像放入黑盒子，它就输出了傅里叶变换结果，提供后续处理步骤的数据。今天，我们就来拆开这个“黑盒子”，看个明白，并且使用C++代码把它写出来。 一、简介目前流行的大多数成熟的FFT算法的基本思路可以分为两大类：一类是按时间抽取的快速傅里叶算法（DIT-FFT），一类是按照频率抽取的快速傅里叶算法（DIF-FFT）。 按时间抽取的FFT 按时间抽取的FFT算法是基于输入序列 按频率抽取的FFT 按频率抽取的FFT算法是基于将输出序列 这次教程我们采用实现较为简单的按时间抽取的FFT，因为针对的序列长度是2的整数次幂，所以，这些算法又称基-2FFT。 （1）、一维按时间抽取的基-2FFT算法对于基-2的FFT，可以设序列长度为 则 根据DFT的周期性和
上式一个递推公式，是FFT的蝶形运算的理论基础。由于我们只关注基-2的FFT算法，所以序列分解到最后都会以2点的DFT进行计算，所以我们先看看2点的FFT蝶形图。
voiddip::CFFT2D::FFT_(std::complex<double> * TD, std::complex<double> *FD, int r) { // 傅立叶变换点数 long count; // 循环变量 int i,j,k; // 中间变量 int bfsize,p; // 角度 double angle; std::complex<double> *W,*X1,*X2,*X; // 计算傅立叶变换点数 count = 1 << r; // 分配运算所需存储器 W = new std::complex<double>[count / 2]; X1 = newstd::complex<double>[count]; X2 = newstd::complex<double>[count]; // 计算加权系数 for(i = 0; i < count / 2; i++) { angle= -i * 3.141****26 * 2 / count; W[i]= std::complex<double> (cos(angle), sin(angle)); } // 将时域点写入X1 memcpy(X1, TD,sizeof(std::complex<double>) * count); // 采用蝶形算法进行快速傅立叶变换 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] = X1[i + p] + X1[i + p + bfsize / 2]; X2[i+ p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)]; } } //输入输出交换，即原位运算 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]=X1[p]; } // 释放内存 delete W; delete X1; delete X2; } （2）、一维按时间抽取的基-2IFFT算法 在得到一维福利叶变换后，我进行一维逆傅里叶变换就变得非常简单了，只要将傅里叶变换结果进行取共轭，再次进行FFT变换，将结果缩放 voiddip::CFFT2D::IFFT_(std::complex<double> * FD, std::complex<double>* TD, int r) { // 傅立叶变换点数 long count; // 循环变量 int i; std::complex<double> *X; // 计算傅立叶变换点数 count = 1 << r; // 分配运算所需存储器 X = newstd::complex<double>[count]; // 将频域点写入X memcpy(X, FD,sizeof(std::complex<double>) * count); // 求共轭 for(i = 0; i < count; i++) { X[i]= std::complex<double> (X[i].real(), -X[i].imag()); } // 调用快速傅立叶变换 FFT_(X, TD, r); // 求时域点的共轭 for(i = 0; i < count; i++) { TD[i]= std::complex<double> (TD[i].real() / count, -TD[i].imag() / count); } // 释放内存 delete X; } （3）、二维FFT与IFFT的实现根据维度拓展的理论，在一维FFT的基础上，我们只要将二维数据在X方向上进行一维FFT变换，然后在X方向变化的结果的基础上，进行Y方向的FFT变换，就可以完成二维FFT。对于二维IFFT，进行类似的操作就可以实现。 在编程实现时，x方向的变换，就是将每一行当作一个序列，由于内存的连续，所以我们只要提供每一行的首地址，就可以访问整行的数据。对于Y方向的变换，我们可以将X方向的变换结果进行矩阵转置，这样每一行的首地址就变成了每一列的首地址，就可以连续的访问每一列了。 二、规划从简介中，我们都是讨论序列尺寸都是2的整数次幂的FFT和IFFT，对于图像而言就是图像的宽度和高度都必须是2的整数次幂。为了使我们的代码更加具有适用性，我们将FFT和IFFT封装成一个C++类，提供可以填充图像、裁剪图像、FFT、IFFT并且可以输出幅度谱图像等接口。 那么我们的类规划如下： 命名空间：dip 类名：CFFT2D 公共接口：
头文件： #ifndefCFFT2D_H_H_H_H_H_H_H_H #defineCFFT2D_H_H_H_H_H_H_H_H #include<complex> namespace dip { class CFFT2D { public: //构造函数 CFFT2D(void); //析构函数 ~CFFT2D(void); //最优尺寸求解 voidGetOptimalSize( int nWidth, int nHeight, int &nOptimalWidth, int&OptimalHeight ); //复数矩阵创建 voidCreateComplexMatrix(unsigned char *pImageData, int nWidth, int nHeight,std::complex<double> *pTimeDomain ); //图像数据提取 voidCreateImageMatrix(std::complex<double> *pInverseResult, int nWidth, intnHeight, unsigned char *pImageData ); //2维离散傅里叶变换 voidForwardTransate(std::complex<double> *pTimeDomain, int nWidth, intnHeight, std::complex<double> *pFrequencyDomain ); //2维离散傅里叶逆变换 voidInverseTransate(std::complex<double> *pFrequencyDomain, int nWidth, intnHeight, std::complex<double> *pTimeDomain); //原点移动 voidShift(std::complex<double> *pSrc, int nWidth, int nHeight); //幅度谱图像求解 voidGetAmplitudeSpectrum( std::complex<double> *pFrequencyDomain, int nWidth,int nHeight, unsigned char *pAmplitudeSpectrum); //相位谱图像求解 voidGetPhaseSpectrum(std::complex<double> *pFrequencyDomain, int nWidth, intnHeight, unsigned char *pPhaseSpectrum); //功率谱图像求解 voidGetPowerSpectrum(std::complex<double> *pFrequencyDomain, int nWidth, intnHeight, unsigned char *pPowerSpectrum); //图像扩充 voidExpand(unsigned char *pImageData, int nWidth, int nHeight, int nExpandWidth, int nExpandHeight,unsigned char *pExpandData ); //图像裁剪 voidCroppe(unsigned char *pSrc, int nWidth, int nHeight, unsigned char *pCroppeDst, intnCropWidth, int nCropHeight); private: //一维傅里叶变换 voidFFT_(std::complex<double> * TD, std::complex<double> * FD, int r); //一维逆傅里叶变换 voidIFFT_(std::complex<double> * FD, std::complex<double> * TD, int r); //计算一个数的与2的n次幂最接近的值 intMinPower2_(int nOrg ); //判断一个数是不是2的n次幂 boolIsPower2_( int nOrg, int &nMaxPower ); }; } #endif 四、测试（1）、原始图像
（2）、扩充图像
（3）、幅度谱图像
（4）、相位谱图像 （5）、逆变换图像
（6）、裁剪图像
五、函数实体#include"StdAfx.h" #include"FFT2D.h" #include<assert.h> #include<limits> /****************************************************** 函数名：CFFT2D 属性：公有 功能：构造函数 参数： 无 返回值： 无 ******************************************************/ dip::CFFT2D::CFFT2D(void) { } /****************************************************** 函数名：~CFFT2D 属性：公有 功能：析构函数 参数： 无 返回值： 无 ******************************************************/ dip::CFFT2D::~CFFT2D(void) { } /****************************************************** 函数名：GetOptimalSize 属性：公有 功能：最优的傅里叶尺寸求解 参数： nWidth -in- 原始宽度 nHeight -in- 原始高度 nOptimalWidth -in- 最优的宽度 OptimalHeight -in- 最优的高度 返回值： 无 ******************************************************/ voiddip::CFFT2D::GetOptimalSize(int nWidth, int nHeight, int &nOptimalWidth,int &OptimalHeight) { nOptimalWidth = MinPower2_( nWidth ); OptimalHeight = MinPower2_( nHeight ); } /****************************************************** 函数名：CreateComplexMatrix 属性：公有 功能：根据图像创建对应的复数矩阵 参数： pImageData -in- 图像数据指针 nWidth -in- 图像宽度 nHeight -in- 图像高度 pTimeDomain -in- 时域复数矩阵 返回值： 无 ******************************************************/ voiddip::CFFT2D::CreateComplexMatrix(unsigned char *pImageData, int nWidth, intnHeight, std::complex<double> *pTimeDomain) { for (int i = 0; i < nHeight; i++) { unsignedchar *pImageRow = pImageData + i * nWidth; std::complex<double>*pTimeRow = pTimeDomain + i * nWidth; for(int j = 0; j < nWidth; j++) { double dReal = pImageRow[j]; pTimeRow[j] =std::complex<double>(dReal, 0); } } } /****************************************************** 函数名：CreateImageMatrix 属性：公有 功能：根据复数矩阵创建对应的图像 参数： pInverseResult -in- 傅里叶逆变换结果 nWidth -in- 图像宽度 nHeight -in- 图像高度 pImageData -in- 图像数据指针 返回值： 无 ******************************************************/ voiddip::CFFT2D::CreateImageMatrix(std::complex<double> *pInverseResult, intnWidth, int nHeight, unsigned char *pImageData) { //初始化最大值最小值 double dMin = (std::numeric_limits<double>::max)(); double dMax =(std::numeric_limits<double>::min)(); //寻找最大值最小值 for (int i = 0; i < nHeight; i++) { std::complex<double>*pResultRow = pInverseResult + i * nWidth; for(int j = 0; j < nWidth; j++) { double dReal = pResultRow[j].real(); if (dMax <dReal) { dMax= dReal; } if (dMin > dReal) { dMin= dReal; } } } //将灰度范围拉伸到0-255 for (int i = 0; i < nHeight; i++) { unsignedchar *pImageRow = pImageData + i * nWidth; std::complex<double>*pResultRow = pInverseResult + i * nWidth; for(int j = 0; j < nWidth; j++) { double dReal = pResultRow[j].real(); pImageRow[j] = (dReal - dMin)/(dMax -dMin) * 255.0; } } } /****************************************************** 函数名：ForwardTransate 属性：公有 功能：2维离散傅里叶变换 参数： pTimeDomain -in- 时域复数矩阵 nWidth -in- 矩阵宽度 nHeight -in- 矩阵高度 pFrequencyDomain -in- 频域复数矩阵 返回值： 无 ******************************************************/ voiddip::CFFT2D::ForwardTransate(std::complex<double> *pTimeDomain, intnWidth, int nHeight, std::complex<double> *pFrequencyDomain) { //参数检查，高度和宽度必须为2的n次幂 int nW = 1, nH = 1; assert(IsPower2_( nWidth, nW) &&IsPower2_(nWidth, nH)); //申请转置缓存 std::complex<double> *pTD = newstd::complex<double>[nWidth * nHeight]; std::complex<double> *pFD = newstd::complex<double>[nWidth * nHeight]; // 对y方向进行快速傅立叶变换 for(int i = 0; i < nHeight; i++) { FFT_(&pTimeDomain[nWidth* i], &pFD[nWidth * i], nW); } //转置 for (int i = 0; i < nHeight; i++) { for(int j = 0; j < nWidth; j++) { pTD[j * nHeight + i] = pFD[ i * nWidth +j]; } } // 对x方向进行快速傅立叶变换 for(int i = 0; i < nWidth; i++) { FFT_(&pTD[i* nHeight], &pFD[i * nHeight], nH); } //转置 for (int i = 0; i < nHeight; i++) { for(int j = 0; j < nWidth; j++) { pFrequencyDomain[i * nWidth + j] =pFD[j * nHeight + i]; } } delete pFD; delete pTD; } 下转《快速傅里叶变换与逆变换（二）》 |
|
|
分解成较短序列然从这些较短的DFT求得输入序列的
的方法。由于抽取后的较短序列可以继续分解下去，继续得到更短的序列，从而可以更简单的进行运算。
。由于N是偶数，按照序列项的奇偶分成2组。
重复周期性（具体推导流程请参考《数字信号处理》），上式可以转为：
