快速傅里叶变换(FFT)的C#实现及详细注释
快速傅里葉變換(FFT)的C#實現及詳細注釋
-------------------------------------------------------------------------------------------------------------------
作者:隨煜而安
時間:2015/7/21
注:本文為作者原創文章,所有參考內容均在參考文獻中列出,轉載請注明出處。
參考文獻:
數字信號處理第三版(華中科技大學出版社)
C#數字圖像處理算法典型實例(人民郵電出版社)
-------------------------------------------------------------------------------------------------------------------
在閱讀下面的快速傅立葉變換實現代碼前,需要先了解FFT的基本原理及過程??梢詤⒖嘉业纳弦黄┛涂焖俑道锶~變換原理解析
下面給出的代碼實現的是一維頻率抽取的基2FFT算法,上面鏈接的博客中已經說過,對于快速傅立葉變換FFT,按時間抽取和按頻率抽取的快速算法的計算量是相同的,所以這里只給出了按頻率抽取的代碼。
首先,我們知道,傅里葉變換的序列和變換結果都應該是復數序列,所以我們要封裝一個復數Complex類。詳細代碼我已經在我的博客C#復數類的封裝中給出(我今天發現其中有個ToString的重載方法寫的有點問題,如果要使用注意一下,但無傷大雅也不影響本博客的任何內容)。
封裝好了復數Complex類,我們就可以開始實現FFT算法了。首先給出兩個私有方法,也是真正實現FFT的核心方法。
一維頻率抽取基2的FFT算法:
/// <summary>/// 一維頻率抽取基2快速傅里葉變換/// 頻率抽取:輸入為自然順序,輸出為碼位倒置順序/// 基2:待變換的序列長度必須為2的整數次冪/// </summary>/// <param name="sourceData">待變換的序列(復數數組)</param>/// <param name="countN">序列長度,可以指定[0,sourceData.Length-1]區間內的任意數值</param>/// <returns>返回變換后的序列(復數數組)</returns>private Complex[] fft_frequency(Complex[] sourceData, int countN){//2的r次冪為N,求出r.r能代表fft算法的迭代次數int r = Convert.ToInt32(Math.Log(countN, 2));//分別存儲蝶形運算過程中左右兩列的結果Complex[] interVar1 = new Complex[countN];Complex[] interVar2 = new Complex[countN];interVar1 = (Complex[])sourceData.Clone();//w代表旋轉因子Complex[] w = new Complex[countN / 2];//為旋轉因子賦值。(在蝶形運算中使用的旋轉因子是已經確定的,提前求出以便調用)//旋轉因子公式 \ /\ /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));}//蝶形運算for (int i = 0; i < r; i++){//i代表當前的迭代次數,r代表總共的迭代次數.//i記錄著迭代的重要信息.通過i可以算出當前迭代共有幾個分組,每個分組的長度//interval記錄當前有幾個組// <<是左移操作符,左移一位相當于*2//多使用位運算符可以人為提高算法速率^_^int interval = 1 << i;//halfN記錄當前循環每個組的長度Nint halfN = 1 << (r - i);//循環,依次對每個組進行蝶形運算for (int j = 0; j < interval; j++){//j代表第j個組//gap=j*每組長度,代表著當前第j組的首元素的下標索引int gap = j * halfN;//進行蝶形運算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];}}//將結果拷貝到輸入端,為下次迭代做好準備interVar1 = (Complex[])interVar2.Clone();}//將輸出碼位倒置for (uint j = 0; j < countN; j++){//j代表自然順序的數組元素的下標索引//用rev記錄j碼位倒置后的結果uint rev = 0;//num作為中間變量uint num = j;//碼位倒置(通過將j的最右端一位最先放入rev右端,然后左移,然后將j的次右端一位放入rev右端,然后左移...)//由于2的r次冪=N,所以任何j可由r位二進制數組表示,循環r次即可for (int i = 0; i < r; i++){rev <<= 1;rev |= num & 1;num >>= 1;}interVar2[rev] = interVar1[j];}return interVar2;}
??一維頻率抽取基2的IFFT算法:
實現IFFT即逆變換的方法有多種,我采用的是能夠重復調用寫好的FFT方法的方式。這種方式需要在進行逆變換前對輸入序列稍作處理取每個元素的共軛。然后調用FFT方法,最終同意對結果做除N處理即可。具體實現代碼如下
/// <summary>/// 一維頻率抽取基2快速傅里葉逆變換/// </summary>/// <param name="sourceData">待反變換的序列(復數數組)</param>/// <param name="countN">序列長度,可以指定[0,sourceData.Length-1]區間內的任意數值</param>/// <returns>返回逆變換后的序列(復數數組)</returns>private Complex[] ifft_frequency(Complex[] sourceData, int countN){//將待逆變換序列取共軛,再調用正變換得到結果,對結果統一再除以變換序列的長度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;}
這樣我們有了核心的兩個方法,當然我們實現的是基2的FFT,對于其他情況,我打算在考完研后補充一個普通DFT的算法,一個針對N為合數的FFT算法。這樣我們就可以封裝一個供用戶調用的公共方法,針對N的類型,智能的選擇合適的算法。當然,目前就先實現到這里了,下面是我封裝的公共方法,N為合數的FFT算法用注釋代替了,虛位以待,考完研后補上!
/// <summary>/// 對給定的序列進行指定長度的離散傅里葉變換DFT/// 內部將使用快速傅里葉變換FFT/// </summary>/// <param name="sourceData">待變換的序列</param>/// <param name="countN">變換的長度N</param>/// <returns>返回變換后的結果(復數數組)</returns>public Complex[] DFT(Complex[] sourceData, int countN){if (countN > sourceData.Length || countN < 0)throw new Exception("指定的傅立葉變換長度越界!");//求出r,2的r次冪為Ndouble dr = Math.Log(countN, 2);int r = Convert.ToInt32(dr);//獲取整數部分//初始化存儲變換結果的數組Complex[] result = new Complex[countN];//判斷選擇合適的算法進行快速傅里葉變換FFTif ((dr - r) != 0){//待變換序列長度不是基2的}else{//待變換序列長度是基2的//使用一維頻率抽取基2快速傅里葉變換result = fft_frequency(sourceData, countN);}return result;}
總結
以上是生活随笔為你收集整理的快速傅里叶变换(FFT)的C#实现及详细注释的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 包管理工具conda极简教程
- 下一篇: 实现Modbus ASCII多主站应用