快速傅里葉變換（FFT）的C#實(shí)現(xiàn)及詳細(xì)注釋

-------------------------------------------------------------------------------------------------------------------

作者：隨煜而安

時(shí)間：2015/7/21

注：本文為作者原創(chuàng)文章，所有參考內(nèi)容均在參考文獻(xiàn)中列出，轉(zhuǎn)載請(qǐng)注明出處。

參考文獻(xiàn)：

數(shù)字信號(hào)處理第三版（華中科技大學(xué)出版社）

C#數(shù)字圖像處理算法典型實(shí)例（人民郵電出版社）

-------------------------------------------------------------------------------------------------------------------

在閱讀下面的快速傅立葉變換實(shí)現(xiàn)代碼前，需要先了解FFT的基本原理及過(guò)程。可以參考我的上一篇博客快速傅里葉變換原理解析

下面給出的代碼實(shí)現(xiàn)的是一維頻率抽取的基2FFT算法，上面鏈接的博客中已經(jīng)說(shuō)過(guò)，對(duì)于快速傅立葉變換FFT，按時(shí)間抽取和按頻率抽取的快速算法的計(jì)算量是相同的，所以這里只給出了按頻率抽取的代碼。

首先，我們知道，傅里葉變換的序列和變換結(jié)果都應(yīng)該是復(fù)數(shù)序列，所以我們要封裝一個(gè)復(fù)數(shù)Complex類。詳細(xì)代碼我已經(jīng)在我的博客C#復(fù)數(shù)類的封裝中給出（我今天發(fā)現(xiàn)其中有個(gè)ToString的重載方法寫的有點(diǎn)問(wèn)題，如果要使用注意一下，但無(wú)傷大雅也不影響本博客的任何內(nèi)容）。

封裝好了復(fù)數(shù)Complex類，我們就可以開始實(shí)現(xiàn)FFT算法了。首先給出兩個(gè)私有方法，也是真正實(shí)現(xiàn)FFT的核心方法。

一維頻率抽取基2的FFT算法：

/// <summary>/// 一維頻率抽取基2快速傅里葉變換/// 頻率抽取：輸入為自然順序，輸出為碼位倒置順序/// 基2：待變換的序列長(zhǎng)度必須為2的整數(shù)次冪/// </summary>/// <param name="sourceData">待變換的序列(復(fù)數(shù)數(shù)組)</param>/// <param name="countN">序列長(zhǎng)度,可以指定[0,sourceData.Length-1]區(qū)間內(nèi)的任意數(shù)值</param>/// <returns>返回變換后的序列（復(fù)數(shù)數(shù)組）</returns>private Complex[] fft_frequency(Complex[] sourceData, int countN){//2的r次冪為N，求出r.r能代表fft算法的迭代次數(shù)int r = Convert.ToInt32(Math.Log(countN, 2));//分別存儲(chǔ)蝶形運(yùn)算過(guò)程中左右兩列的結(jié)果Complex[] interVar1 = new Complex[countN];Complex[] interVar2 = new Complex[countN];interVar1 = (Complex[])sourceData.Clone();//w代表旋轉(zhuǎn)因子Complex[] w = new Complex[countN / 2];//為旋轉(zhuǎn)因子賦值。（在蝶形運(yùn)算中使用的旋轉(zhuǎn)因子是已經(jīng)確定的，提前求出以便調(diào)用）//旋轉(zhuǎn)因子公式 \ /\ /k __// \/ \/N -- exp(-j*2πk/N)//這里還用到了歐拉公式for (int i = 0; i < countN / 2; i++){double angle = -i * Math.PI * 2 / countN;w[i] = new Complex(Math.Cos(angle), Math.Sin(angle));}//蝶形運(yùn)算for (int i = 0; i < r; i++){//i代表當(dāng)前的迭代次數(shù)，r代表總共的迭代次數(shù).//i記錄著迭代的重要信息.通過(guò)i可以算出當(dāng)前迭代共有幾個(gè)分組，每個(gè)分組的長(zhǎng)度//interval記錄當(dāng)前有幾個(gè)組// <<是左移操作符，左移一位相當(dāng)于*2//多使用位運(yùn)算符可以人為提高算法速率^_^int interval = 1 << i;//halfN記錄當(dāng)前循環(huán)每個(gè)組的長(zhǎng)度Nint halfN = 1 << (r - i);//循環(huán)，依次對(duì)每個(gè)組進(jìn)行蝶形運(yùn)算for (int j = 0; j < interval; j++){//j代表第j個(gè)組//gap=j*每組長(zhǎng)度，代表著當(dāng)前第j組的首元素的下標(biāo)索引int gap = j * halfN;//進(jìn)行蝶形運(yùn)算for (int k = 0; k < halfN / 2; k++){interVar2[k + gap] = interVar1[k + gap] + interVar1[k + gap + halfN / 2];interVar2[k + halfN / 2 + gap] = (interVar1[k + gap] - interVar1[k + gap + halfN / 2]) * w[k * interval];}}//將結(jié)果拷貝到輸入端，為下次迭代做好準(zhǔn)備interVar1 = (Complex[])interVar2.Clone();}//將輸出碼位倒置for (uint j = 0; j < countN; j++){//j代表自然順序的數(shù)組元素的下標(biāo)索引//用rev記錄j碼位倒置后的結(jié)果uint rev = 0;//num作為中間變量uint num = j;//碼位倒置（通過(guò)將j的最右端一位最先放入rev右端，然后左移，然后將j的次右端一位放入rev右端，然后左移...）//由于2的r次冪=N，所以任何j可由r位二進(jìn)制數(shù)組表示，循環(huán)r次即可for (int i = 0; i < r; i++){rev <<= 1;rev |= num & 1;num >>= 1;}interVar2[rev] = interVar1[j];}return interVar2;}

??一維頻率抽取基2的IFFT算法：

實(shí)現(xiàn)IFFT即逆變換的方法有多種，我采用的是能夠重復(fù)調(diào)用寫好的FFT方法的方式。這種方式需要在進(jìn)行逆變換前對(duì)輸入序列稍作處理取每個(gè)元素的共軛。然后調(diào)用FFT方法，最終同意對(duì)結(jié)果做除N處理即可。具體實(shí)現(xiàn)代碼如下

/// <summary>/// 一維頻率抽取基2快速傅里葉逆變換/// </summary>/// <param name="sourceData">待反變換的序列（復(fù)數(shù)數(shù)組）</param>/// <param name="countN">序列長(zhǎng)度,可以指定[0,sourceData.Length-1]區(qū)間內(nèi)的任意數(shù)值</param>/// <returns>返回逆變換后的序列（復(fù)數(shù)數(shù)組）</returns>private Complex[] ifft_frequency(Complex[] sourceData, int countN){//將待逆變換序列取共軛，再調(diào)用正變換得到結(jié)果，對(duì)結(jié)果統(tǒng)一再除以變換序列的長(zhǎng)度Nfor (int i = 0; i < countN; i++){sourceData[i] = sourceData[i].Conjugate();}Complex[] interVar = new Complex[countN];interVar = fft_frequency(sourceData, countN);for (int i = 0; i < countN; i++){interVar[i] = new Complex(interVar[i].Real / countN, -interVar[i].Imaginary / countN);}return interVar;}

這樣我們有了核心的兩個(gè)方法，當(dāng)然我們實(shí)現(xiàn)的是基2的FFT，對(duì)于其他情況，我打算在考完研后補(bǔ)充一個(gè)普通DFT的算法，一個(gè)針對(duì)N為合數(shù)的FFT算法。這樣我們就可以封裝一個(gè)供用戶調(diào)用的公共方法，針對(duì)N的類型，智能的選擇合適的算法。當(dāng)然，目前就先實(shí)現(xiàn)到這里了，下面是我封裝的公共方法，N為合數(shù)的FFT算法用注釋代替了，虛位以待，考完研后補(bǔ)上！

/// <summary>/// 對(duì)給定的序列進(jìn)行指定長(zhǎng)度的離散傅里葉變換DFT/// 內(nèi)部將使用快速傅里葉變換FFT/// </summary>/// <param name="sourceData">待變換的序列</param>/// <param name="countN">變換的長(zhǎng)度N</param>/// <returns>返回變換后的結(jié)果（復(fù)數(shù)數(shù)組）</returns>public Complex[] DFT(Complex[] sourceData, int countN){if (countN > sourceData.Length || countN < 0)throw new Exception("指定的傅立葉變換長(zhǎng)度越界！");//求出r,2的r次冪為Ndouble dr = Math.Log(countN, 2);int r = Convert.ToInt32(dr);//獲取整數(shù)部分//初始化存儲(chǔ)變換結(jié)果的數(shù)組Complex[] result = new Complex[countN];//判斷選擇合適的算法進(jìn)行快速傅里葉變換FFTif ((dr - r) != 0){//待變換序列長(zhǎng)度不是基2的}else{//待變換序列長(zhǎng)度是基2的//使用一維頻率抽取基2快速傅里葉變換result = fft_frequency(sourceData, countN);}return result;}

總結(jié)

以上是生活随笔為你收集整理的快速傅里叶变换（FFT）的C#实现及详细注释的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。