快速傅立葉變換(FFT)源程序

1、/實(shí)驗(yàn)用的頭文件 MYFFT.H /作用：為幫助小虎子做實(shí)驗(yàn)，這個(gè)頭文件提供了完整的一維與二維FFT算法，我想應(yīng)改是夠你折騰了吧！ #include / complex using namespace std; typedef complex Comp; / 復(fù)數(shù)類型定義 const float _2PI_ = 2.0f * 3.14159265f; / 常數(shù)2PI定義 const int MAX_N = 256; / 最大DFT點(diǎn)數(shù) /*-*-*-*-*-*-*-*-*-*-*-*-* FFT算法模塊接口定義 *-*-*-*-*-*-*-*-*-*-*-*-*/ / / Function n

2、ame : BitReverse / Description : 二進(jìn)制倒序操作 / Return type : int / Argument : int src 待倒讀的數(shù) / Argument : int size 二進(jìn)制位數(shù) int BitReverse(int src, int size) int tmp = src; int des = 0; for (int i=size-1; i=0; i-) des = (tmp & 0x1) 1; return des; / / Function name : Reorder / Description : 數(shù)據(jù)二進(jìn)制整序 / Return

3、type : void / Argument : Comp xMAX_N 待整序數(shù)組 / Argument : int N FFT點(diǎn)數(shù) / Argument : int M 點(diǎn)數(shù)的2的冪次 void Reorder(Comp xMAX_N, int N, int M) Comp new_xMAX_N; for (int i=0; iN; i+) new_x = xBitReverse(i, M); / 重新存入原數(shù)據(jù)中(已經(jīng)是二進(jìn)制整序過了的數(shù)據(jù)) for (i=0; iN; i+) x = new_x; / / Function name : CalcW / Description : 計(jì)算

4、旋轉(zhuǎn)因子 / Return type : void / Argument : Comp WMAX_N 存放因子的數(shù)組 / Argument : int N FFT的點(diǎn)數(shù) / Argument : int flag 正逆變換標(biāo)記 void CalcW(Comp WMAX_N, int N, int flag) for (int i=0; iN/2; i+) if (flag = 1) W = exp(Comp(0, -_2PI_ * i / N); / 正FFT變換 else W = exp(Comp(0, _2PI_ * i / N); / 逆FFT變換 / / Function name :

5、 FFT_1D_Kernel / Description : FFT算法核心 / Return type : void / Argument : Comp* x 數(shù)據(jù) / Argument : int M 冪次 / Argument : int flag 正逆變換標(biāo)記 以下本應(yīng)由自己完成。 void FFT_1D(Comp* x, int M, int flag) int N = (1 M); / 二進(jìn)制整序 Reorder(x, N, M); / 旋轉(zhuǎn)因子計(jì)算 Comp WMAX_N; CalcW(W, N, flag); / 級(jí)內(nèi)群數(shù) int GroupNum = N/2; / 第一級(jí)的群

6、數(shù)為N/2 / 群內(nèi)蝶形單元數(shù) int CellNum = 1; / 第一級(jí)的群內(nèi)蝶形單元數(shù)為1 / 處理各級(jí) for (int i=0; iM; i+) / 處理各群 for (int j=0; jGroupNum; j+) / 處理各蝶形單元 for (int k=0; k= 1; / 級(jí)別增加, 則相應(yīng)的群數(shù)減少一半 CellNum = 1; / 級(jí)別增加, 則相應(yīng)的群內(nèi)單元數(shù)增加一倍 / 如果是IFFT，各元素還要再除以N if (flag != 1) for (i=0; iN; i+) x /= N; / / Function name : FFT_2D_Kernel / Descr

7、iption : 2D FFT核心算法 / Return type : void / Argument : Comp xMAX_NMAX_N 二維數(shù)據(jù) / Argument : int M 冪次 / Argument : int flag 正逆變換標(biāo)記 以下本應(yīng)由自己完成。 void FFT_2D(Comp xMAX_NMAX_N, int M, int flag) int N = (1 M); / 先逐行進(jìn)行 1D-FFT for (int i=0; iN; i+) FFT_1D(x, M, flag); / - 計(jì)算結(jié)果再存入矩陣x中 / 再逐列進(jìn)行 1D-FFT Comp colMAX_

8、N; for (int j=0; jN; j+) / 取得第j列的數(shù)據(jù) for (int i=0; iN; i+) col = xj; / 對(duì)第j列數(shù)據(jù)進(jìn)行 1D-FFT FFT_1D(col, M, flag); / - 計(jì)算結(jié)果在數(shù)組col中 / 將結(jié)果放回矩陣第j列中 for (i=0; iN; i+) xj = col; / 0),實(shí)質(zhì)上k是n個(gè)采樣數(shù)據(jù)可以分解為偶次冪和奇次冪的次數(shù)/ pr: l=0時(shí)，存放N點(diǎn)采樣數(shù)據(jù)的實(shí)部/ l=1時(shí), 存放傅立葉變換的N個(gè)實(shí)部/ pi: l=0時(shí)，存放N點(diǎn)采樣數(shù)據(jù)的虛部 / l=1時(shí), 存放傅立葉變換的N個(gè)虛部/ 出口參數(shù)：/ fr: l=0,

9、返回傅立葉變換的實(shí)部/ l=1, 返回逆傅立葉變換的實(shí)部/ fi: l=0, 返回傅立葉變換的虛部/ l=1, 返回逆傅立葉變換的虛部/ pr: il = 1,i = 0 時(shí)，返回傅立葉變換的模/ il = 1,i = 1 時(shí)，返回逆傅立葉變換的模/ pi: il = 1,i = 0 時(shí)，返回傅立葉變換的輻角/ il = 1,i = 1 時(shí)，返回逆傅立葉變換的輻角/ data: 2005.8.15,Mend Xin Dongvoid kkfft(double pr, double pi, int n, int k, double fr, double fi, int l, int il)int

10、 it,m,is,i,j,nv,l0;double p,q,s,vr,vi,poddr,poddi;for (it=0; it=n-1; it+) m = it; is = 0;for(i=0; i=k-1; i+) j = m/2; is = 2*is+(m-2*j); m = j;frit = pris; fiit = piis;/-pr0 = 1.0; pi0 = 0.0;p = 6.283185306/(1.0*n);pr1 = cos(p); pi1 = -sin(p);if (l!=0) pi1=-pi1;for (i=2; i=n-1; i+) p = pri-1*pr1; q

11、= pii-1*pi1;s = (pri-1+pii-1)*(pr1+pi1);pri = p-q; pii = s-p-q;for (it=0; it=0; l0-) m = m/2; nv = 2*nv;for(it=0; it=(m-1)*nv; it=it+nv)for (j=0; j=(nv/2)-1; j+) p = prm*j*frit+j+nv/2;q = pim*j*fiit+j+nv/2;s = prm*j+pim*j;s = s*(frit+j+nv/2+fiit+j+nv/2);poddr = p-q; poddi = s-p-q;frit+j+nv/2 = frit+

12、j-poddr;fiit+j+nv/2 = fiit+j-poddi;frit+j = frit+j+poddr;fiit+j = fiit+j+poddi;if(l!=0)for(i=0; i=n-1; i+) fri = fri/(1.0*n);fii = fii/(1.0*n);if(il!=0)for(i=0; i=n-1; i+) pri = sqrt(fri*fri+fii*fii);if(fabs(fri)0) pii = 90.0;else pii = -90.0;elsepii = atan(fii/fri)*360.0/6.283185306;return;/快速傅立葉變換

13、的源代碼/-/ fft.cpp/ The implementation of the/ Fast Fourier Transform algorithm/ (c) Reliable Software, 1996/-#include fft.h#include recorder.h/ log (1) = 0, log(2) = 1, log(3) = 2, log(4) = 2 .#define PI (2.0 * asin(1.0)/ Points must be a power of 2Fft:Fft (int Points, long sampleRate): _Points (Point

14、s), _sampleRate (sampleRate)_aTape = new double _Points;#if 0/ 1 kHz calibration wavefor (int i = 0; i _Points; i+)_aTapei = 1600 * sin (2 * PI * 1000. * i / _sampleRate);#elsefor (int i = 0; i = 1;_logPoints+;_aBitRev = new int _Points;_X = new Complex_Points;_W = new Complex* _logPoints+1;/ Precom

15、pute complex exponentialsint _2_l = 2;for (int l = 1; l = _logPoints; l+)_Wl = new Complex _Points;for ( int i = 0; i _Points; i+ )double re = cos (2. * PI * i / _2_l);double im = -sin (2. * PI * i / _2_l);_Wli = Complex (re, im);_2_l *= 2;/ set up bit reverse mappingint rev = 0;int halfPoints = _Po

16、ints/2;for (i = 0; i = mask)rev -= mask; / turn off this bitmask = 1;rev += mask;_aBitRev _Points-1 = _Points-1;Fft:Fft()delete _aTape;delete _aBitRev;for (int l = 1; l _Points)return;/ make space for cSample samples at the end of tape/ shifting previous samples towards the beginningmemmove (_aTape,

17、 &_aTapecSample,(_Points - cSample) * sizeof(double);/ copy samples from iterator to tail end of tapeint iTail = _Points - cSample;for (int i = 0; i cSample; i+, iter.Advance()_aTape i + iTail = (double) iter.GetSample();/ Initialize the FFT bufferfor (i = 0; i _Points; i+)PutAt (i, _aTapei);/ 0 1 2 3 4 5 6 7/ level 1/ step 1 0/ increm 2 W/ j = 0 1/ level 2/ step 2/ increm 4 0/ j = 0 W 1/ j = 1 2 W/ level 3 2/ step 4/ increm 8 0/ j = 0 W 1/ j = 1 3 W 2/ j = 2 3 W 3/ j = 3 3 W/ 3/void Fft:Transform ()/ step = 2 (level-1)/ increm = 2 level;int step = 1