#include "FFT_code.h"
//比特翻转
void bitrp (vector<complex<double>> &Data_In, int n)
{
// 位反转置换 Bit-reversal Permutation
int i, j, a, b, p;
complex<double> temp;
for (i = 1, p = 0; i < n; i *= 2)
{
p ++;
}
for (i = 0; i < n; i ++)
{
a = i;
b = 0;
for (j = 0; j < p; j ++)
{
b = (b << 1) + (a & 1); // b = b * 2 + a % 2;
a >>= 1; // a = a / 2;
} //令b为i的p位二进制数中心翻转后的数
if ( b > i)
{
temp=Data_In[b];
Data_In[b]=Data_In[i];
Data_In[i]=temp;
} //交换原数据第i位和b位的值
}
}
//FFT
int FFT(vector<complex<double>> &Data_in, int n)
{
int N=n,i,p;
complex<double> temp=(0.0,0.0);
if( N& (N-1) ) //判断N是否为2的幂
{
for (i = 1, p = 0; i < N; i *= 2)
{
p ++;
}
N=pow(2.0,p);
} //找到比n大的N，N是离n最近的一个2的幂
for(int len=0;len<N-n;len++)
Data_in.push_back(temp); //补0使数据长度为N
double* wreal=new double[N/2];
double* wimag=new double[N/2];
double treal, timag, ureal, uimag, arg;
int m, k, j, t, index1, index2;
bitrp (Data_in, N); //比特翻转
// 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
arg = - 2 * 3.14159 / N;
treal = cos (arg);
timag = sin (arg);
wreal [0] = 1.0;
wimag [0] = 0.0;
for (j = 1; j < N / 2; j ++)
{
wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
}
for (m = 2; m <= N; m *= 2)
{
for (k = 0; k < N; k += m)
{
for (j = 0; j < m / 2; j ++)
{
index1 = k + j;
index2 = index1 + m / 2;
t = N * j / m; // 旋转因子w的实部在wreal[]中的下标为t
treal = wreal [t] * Data_in[index2].real() - wimag [t] * Data_in[index2].imag();
timag = wreal [t] * Data_in[index2].imag() + wimag [t] * Data_in[index2].real();
ureal = Data_in[index1].real();
uimag = Data_in[index1].imag();
Data_in[index1].real(ureal + treal);
Data_in[index1].imag(uimag + timag);
Data_in[index2].real(ureal - treal);
Data_in[index2].imag(uimag - timag);
}
}
}
delete[] wreal;
delete[] wimag;
return N;
}
评论0