快速傅立叶变换(FFT)源程序

1、/实验用的头文件 MYFFT.H /作用：为帮助小虎子做实验，这个头文件提供了完整的一维与二维FFT算法，我想应改是够你折腾了吧！ #include / complex using namespace std; typedef complex Comp; / 复数类型定义 const float _2PI_ = 2.0f * 3.14159265f; / 常数2PI定义 const int MAX_N = 256; / 最大DFT点数 /*-*-*-*-*-*-*-*-*-*-*-*-* FFT算法模块接口定义 *-*-*-*-*-*-*-*-*-*-*-*-*/ / / Function n

2、ame : BitReverse / Description : 二进制倒序操作 / Return type : int / Argument : int src 待倒读的数 / Argument : int size 二进制位数 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 : 数据二进制整序 / Return

3、type : void / Argument : Comp xMAX_N 待整序数组 / Argument : int N FFT点数 / Argument : int M 点数的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); / 重新存入原数据中(已经是二进制整序过了的数据) for (i=0; iN; i+) x = new_x; / / Function name : CalcW / Description : 计算

4、旋转因子 / Return type : void / Argument : Comp WMAX_N 存放因子的数组 / Argument : int N FFT的点数 / Argument : int flag 正逆变换标记 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 数据 / Argument : int M 幂次 / Argument : int flag 正逆变换标记 以下本应由自己完成。 void FFT_1D(Comp* x, int M, int flag) int N = (1 M); / 二进制整序 Reorder(x, N, M); / 旋转因子计算 Comp WMAX_N; CalcW(W, N, flag); / 级内群数 int GroupNum = N/2; / 第一级的群

6、数为N/2 / 群内蝶形单元数 int CellNum = 1; / 第一级的群内蝶形单元数为1 / 处理各级 for (int i=0; iM; i+) / 处理各群 for (int j=0; jGroupNum; j+) / 处理各蝶形单元 for (int k=0; k= 1; / 级别增加, 则相应的群数减少一半 CellNum = 1; / 级别增加, 则相应的群内单元数增加一倍 / 如果是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 二维数据 / Argument : int M 幂次 / Argument : int flag 正逆变换标记 以下本应由自己完成。 void FFT_2D(Comp xMAX_NMAX_N, int M, int flag) int N = (1 M); / 先逐行进行 1D-FFT for (int i=0; iN; i+) FFT_1D(x, M, flag); / - 计算结果再存入矩阵x中 / 再逐列进行 1D-FFT Comp colMAX_

8、N; for (int j=0; jN; j+) / 取得第j列的数据 for (int i=0; iN; i+) col = xj; / 对第j列数据进行 1D-FFT FFT_1D(col, M, flag); / - 计算结果在数组col中 / 将结果放回矩阵第j列中 for (i=0; iN; i+) xj = col; / 0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数/ pr: l=0时，存放N点采样数据的实部/ l=1时, 存放傅立叶变换的N个实部/ pi: l=0时，存放N点采样数据的虚部 / l=1时, 存放傅立叶变换的N个虚部/ 出口参数：/ fr: l=0,

9、返回傅立叶变换的实部/ l=1, 返回逆傅立叶变换的实部/ fi: l=0, 返回傅立叶变换的虚部/ l=1, 返回逆傅立叶变换的虚部/ pr: il = 1,i = 0 时，返回傅立叶变换的模/ il = 1,i = 1 时，返回逆傅立叶变换的模/ pi: il = 1,i = 0 时，返回傅立叶变换的辐角/ il = 1,i = 1 时，返回逆傅立叶变换的辐角/ 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