時間抽出(DIT)のベース-2 FFTアルゴリズム(クリー-グラフベースアルゴリズム)C++プログラム

ベース-2 FFTアルゴリズムのC++プログラムは,時間抽出,入力逆位順,出力自然順,N=2 L N=2^L N=2 L

#include 

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+b i、利用可能点Z(a,b)はz=a+bi、利用可能点Z(a,b)はz=a+bi、利用可能点Z(a,b)はz=a+bi

// complex constructor example
#include      // std::cout
#include       // 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 
#define _USE_MATH_DEFINES //    #include      
#include 

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;
}