快速傅里葉變換與逆變換——CFFT2D類的實現(xiàn) 前言在快速傅里葉變換（FFT）未出之前，傅里葉變換并未大為流行，主要的原因就是變換需要的計算時間過多。在FFT出現(xiàn)后，傅里葉變換才算是真正的運用到了實際生產(chǎn)生活當中,。所以，可以這樣說是FFT推廣了傅里葉變換。FFT變換是頻率域圖像處理的基石,，所以對于數(shù)字圖像處理學**來說，學**FFT算法是非常必須的,。雖然,，現(xiàn)在Internet有很多的商業(yè)、開源的代碼庫提供了FFT算法實現(xiàn),，更有甚者,，在DSP的配套硬件中，有專門的FFT實現(xiàn)硬件,，運行速度得到更高的提升,。但是，對于學**來說確實不利的,，就好像一個“黑盒子”一樣,，我把圖像放入黑盒子，它就輸出了傅里葉變換結(jié)果,，提供后續(xù)處理步驟的數(shù)據(jù),。今天,，我們就來拆開這個“黑盒子”，看個明白,，并且使用C++代碼把它寫出來,。 一、簡介目前流行的大多數(shù)成熟的FFT算法的基本思路可以分為兩大類：一類是按時間抽取的快速傅里葉算法（DIT-FFT）,，一類是按照頻率抽取的快速傅里葉算法（DIF-FFT）,。 按時間抽取的FFT 按時間抽取的FFT算法是基于輸入序列 按頻率抽取的FFT 按頻率抽取的FFT算法是基于將輸出序列 這次教程我們采用實現(xiàn)較為簡單的按時間抽取的FFT,，因為針對的序列長度是2的整數(shù)次冪,，所以，這些算法又稱基-2FFT,。 （1）,、一維按時間抽取的基-2FFT算法對于基-2的FFT，可以設(shè)序列長度為 則 根據(jù)DFT的周期性和
上式一個遞推公式，是FFT的蝶形運算的理論基礎(chǔ),。由于我們只關(guān)注基-2的FFT算法,，所以序列分解到最后都會以2點的DFT進行計算，所以我們先看看2點的FFT蝶形圖,。
voiddip::CFFT2D::FFT_(std::complex<double> * TD, std::complex<double> *FD, int r) { // 傅立葉變換點數(shù) long count; // 循環(huán)變量 int i,j,k; // 中間變量 int bfsize,p; // 角度 double angle; std::complex<double> *W,*X1,*X2,*X; // 計算傅立葉變換點數(shù) count = 1 << r; // 分配運算所需存儲器 W = new std::complex<double>[count / 2]; X1 = newstd::complex<double>[count]; X2 = newstd::complex<double>[count]; // 計算加權(quán)系數(shù) 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]; } // 釋放內(nèi)存 delete W; delete X1; delete X2; } （2）、一維按時間抽取的基-2IFFT算法 在得到一維福利葉變換后,，我進行一維逆傅里葉變換就變得非常簡單了，只要將傅里葉變換結(jié)果進行取共軛,，再次進行FFT變換,，將結(jié)果縮放 voiddip::CFFT2D::IFFT_(std::complex<double> * FD, std::complex<double>* TD, int r) { // 傅立葉變換點數(shù) long count; // 循環(huán)變量 int i; std::complex<double> *X; // 計算傅立葉變換點數(shù) 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()); } // 調(diào)用快速傅立葉變換 FFT_(X, TD, r); // 求時域點的共軛 for(i = 0; i < count; i++) { TD[i]= std::complex<double> (TD[i].real() / count, -TD[i].imag() / count); } // 釋放內(nèi)存 delete X; } （3）,、二維FFT與IFFT的實現(xiàn)根據(jù)維度拓展的理論,，在一維FFT的基礎(chǔ)上，我們只要將二維數(shù)據(jù)在X方向上進行一維FFT變換,，然后在X方向變化的結(jié)果的基礎(chǔ)上,，進行Y方向的FFT變換，就可以完成二維FFT,。對于二維IFFT,，進行類似的操作就可以實現(xiàn)。 在編程實現(xiàn)時,，x方向的變換,，就是將每一行當作一個序列，由于內(nèi)存的連續(xù),，所以我們只要提供每一行的首地址,，就可以訪問整行的數(shù)據(jù)。對于Y方向的變換,，我們可以將X方向的變換結(jié)果進行矩陣轉(zhuǎn)置,，這樣每一行的首地址就變成了每一列的首地址，就可以連續(xù)的訪問每一列了,。 二,、規(guī)劃從簡介中，我們都是討論序列尺寸都是2的整數(shù)次冪的FFT和IFFT,，對于圖像而言就是圖像的寬度和高度都必須是2的整數(shù)次冪,。為了使我們的代碼更加具有適用性，我們將FFT和IFFT封裝成一個C++類,，提供可以填充圖像,、裁剪圖像、FFT,、IFFT并且可以輸出幅度譜圖像等接口,。 那么我們的類規(guī)劃如下： 命名空間：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: //構(gòu)造函數(shù) CFFT2D(void); //析構(gòu)函數(shù) ~CFFT2D(void); //最優(yōu)尺寸求解 voidGetOptimalSize( int nWidth, int nHeight, int &nOptimalWidth, int&OptimalHeight ); //復數(shù)矩陣創(chuàng)建 voidCreateComplexMatrix(unsigned char *pImageData, int nWidth, int nHeight,std::complex<double> *pTimeDomain ); //圖像數(shù)據(jù)提取 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); //計算一個數(shù)的與2的n次冪最接近的值 intMinPower2_(int nOrg ); //判斷一個數(shù)是不是2的n次冪 boolIsPower2_( int nOrg, int &nMaxPower ); }; } #endif 四、測試（1）,、原始圖像
（2）,、擴充圖像
（3）、幅度譜圖像
（4）,、相位譜圖像 （5）,、逆變換圖像
（6）、裁剪圖像
五,、函數(shù)實體#include"StdAfx.h" #include"FFT2D.h" #include<assert.h> #include<limits> /****************************************************** 函數(shù)名：CFFT2D 屬性：公有 功能：構(gòu)造函數(shù) 參數(shù)： 無 返回值： 無 ******************************************************/ dip::CFFT2D::CFFT2D(void) { } /****************************************************** 函數(shù)名：~CFFT2D 屬性：公有 功能：析構(gòu)函數(shù) 參數(shù)： 無 返回值： 無 ******************************************************/ dip::CFFT2D::~CFFT2D(void) { } /****************************************************** 函數(shù)名：GetOptimalSize 屬性：公有 功能：最優(yōu)的傅里葉尺寸求解 參數(shù)： nWidth -in- 原始寬度 nHeight -in- 原始高度 nOptimalWidth -in- 最優(yōu)的寬度 OptimalHeight -in- 最優(yōu)的高度 返回值： 無 ******************************************************/ voiddip::CFFT2D::GetOptimalSize(int nWidth, int nHeight, int &nOptimalWidth,int &OptimalHeight) { nOptimalWidth = MinPower2_( nWidth ); OptimalHeight = MinPower2_( nHeight ); } /****************************************************** 函數(shù)名：CreateComplexMatrix 屬性：公有 功能：根據(jù)圖像創(chuàng)建對應的復數(shù)矩陣 參數(shù)： pImageData -in- 圖像數(shù)據(jù)指針 nWidth -in- 圖像寬度 nHeight -in- 圖像高度 pTimeDomain -in- 時域復數(shù)矩陣 返回值： 無 ******************************************************/ 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); } } } /****************************************************** 函數(shù)名：CreateImageMatrix 屬性：公有 功能：根據(jù)復數(shù)矩陣創(chuàng)建對應的圖像 參數(shù)： pInverseResult -in- 傅里葉逆變換結(jié)果 nWidth -in- 圖像寬度 nHeight -in- 圖像高度 pImageData -in- 圖像數(shù)據(jù)指針 返回值： 無 ******************************************************/ 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; } } } /****************************************************** 函數(shù)名：ForwardTransate 屬性：公有 功能：2維離散傅里葉變換 參數(shù)： pTimeDomain -in- 時域復數(shù)矩陣 nWidth -in- 矩陣寬度 nHeight -in- 矩陣高度 pFrequencyDomain -in- 頻域復數(shù)矩陣 返回值： 無 ******************************************************/ voiddip::CFFT2D::ForwardTransate(std::complex<double> *pTimeDomain, intnWidth, int nHeight, std::complex<double> *pFrequencyDomain) { //參數(shù)檢查,，高度和寬度必須為2的n次冪 int nW = 1, nH = 1; assert(IsPower2_( nWidth, nW) &&IsPower2_(nWidth, nH)); //申請轉(zhuǎn)置緩存 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); } //轉(zhuǎn)置 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); } //轉(zhuǎn)置 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; } 下轉(zhuǎn)《快速傅里葉變換與逆變換（二）》 |
|
|
分解成較短序列然從這些較短的DFT求得輸入序列的
的方法。由于抽取后的較短序列可以繼續(xù)分解下去,，繼續(xù)得到更短的序列,，從而可以更簡單的進行運算。
,。由于N是偶數(shù),，按照序列項的奇偶分成2組。
重復周期性（具體推導流程請參考《數(shù)字信號處理》）,，上式可以轉(zhuǎn)為：
