按时间抽选（DIT）的基-2 FFT算法（库利-图基算法）C++程序

基-2 FFT算法的C++程序，按时间抽选、输入倒位序、输出自然顺序，$N=2^L$

#include <complex>

int fft(complex<double> *a,int l)
{
	const double pai = 3.141592653589793;
	complex<double> u,w,t;
	unsigned n=1,nv2,nm1,k,le,lei,ip;
	unsigned i,j,m;
	double tmp;
	
    n=n<<l; // 2^l
    nv2=n>>1; // N/2
    nm1=n-1; // 0,1,2...N-1
    i=0;
    j = 0; //倒位序i,j都从0开始
    
	for (i = 0; i < nm1; i++)
	{
		if (i < j)
		{
			t = a[j];
			a[j] = a[i];
			a[i] = t;
		}
		k = nv2;
        while (j >= k)
        {
            j-=k; //j=j-k
            k>>=1; //k/2
        };
        j+=k; //j=j+k
	}
	le=1;
		for (m = 1; m <= 1; m++)
		{
			lei = le;
			le <<= 1; //2*le
			u = complex<double>(1, 0);
			tmp = pai / lei;
			w = complex<double>(cos(tmp), -sin(tmp));
			for (j = 0; j < lei; j++)
			{
				for (i = j;i<n;i+=le)
				{
					ip = i + lei;
					t = a[ip] * u;
					a[ip] = a[i] - t;
					a[i] += t;
				}
				u*=w;
			}
		}
	return 0;
}

有助于对基-2 FFT算法的理解和掌握，对入门时域/频域处理大有裨益。

.〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓.
用矩阵折半重排理解倒位序算法

a[8]={a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7]}

//设数组a[8]的倒位序为X[8]
//一维矩阵M(8,1)
M=
  0
  1
  2
  3
  4
  5
  6
  7

//经过L级的折半重排，将列向量M(8,1)变换为行向量M_L(1,8)
M_1=     M_2=         M_3=
    0 4      0 4 2 6      0 4 2 6 1 5 3 7
    1 5      1 5 3 7
    2 6
    3 7

//循环结束时得到L=3，与序列点数N=2^L符合
X[8]={a[0],a[4],a[2],a[6],a[1],a[5],a[3],a[7]}

.〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓.
复数$z=a+bi，可用点Z(a，b)表示$

// complex constructor example
#include <iostream>     // std::cout
#include <complex>      // std::complex
using namespace std;

int main()
{
	complex<double> first(2.0, 2.0);
	complex<double> second(first);
	complex<long double> third(second);

	cout << "first= " << first << endl;
	cout << "second= " << second << endl;
	cout << "third= " << third  << endl;
	return 0;
}

.〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓.
对8位序列进行FFT处理

#include <complex>
#define _USE_MATH_DEFINES //要求在 #include <math.h> 前面定义
#include <math.h>

using namespace std;

int FFT(complex<double> array[])
{
	complex<double> u, w, wn, t;
	unsigned n = 1, l = 3, nv2, nm1, k, le, lei;
	unsigned i, j, m;

	//unsigned l = sizeof(array);

	n <<= l; // 2^l
	nv2 = n >> 1; // N/2
	nm1 = n - 1; // 0,1,2...N-1
	i = 0;
	j = 0; //倒位序i,j都从0开始

	for (i = 0; i < nm1; i++)
	{
		if (i < j)
		{
			t = array[j];
			array[j] = array[i];
			array[i] = t;
		}
		k = nv2;
		while (j >= k)
		{
			j -= k; //j=j-k
			k >>= 1; //k/2
		};
		j += k; //j=j+k
	}
	for (int ii = 0; ii < 8; ii++) //下标循环，从0到7
	{
		cout << "I = " << ii << "~" << array[ii] << endl; //输出数组元素值，每行一个
	}
	le = 1;
	for (i = 1; i <= l; i++)
	{
		le <<= 1; //le = 2*le;
		cout << "i= " << i << "\t" << "le= " << le << endl; //与使用<< "\tle= " << le 效果一样，但更清晰
		lei = le / 2;
		//初始化单位复根
		complex<double> wn(cos(2 * M_PI / le), -sin(2 * M_PI / le)); //使用math.h中对pi定义宏
		for (j = 0; j < n; j = j + le)
		{
			complex<double> w(1, 0); // 初始化螺旋因子
			for (m = j; m < j + le / 2; m++)
			{
				cout << "m= " << m << "w= " << w << endl;
				u = array[m];
				t = w*array[m + lei];
				array[m] = u + t;
				array[m + lei] = u - t;
				w = w*wn;
			}
		}
	}
	for (i = 0; i < 8; i++) //下标循环，从0到7
	{
		cout << "I = " << i << "~" << array[i] << endl; //输出数组元素值，每行一个
	}
	cout << endl;
	system("pause"); //暂停，方便查看结果
	return 0;
}