C++实现DIT-FFT

代码思路

序列倒序

蝶型运算

代码技巧

原位计算

旋转因子的变化规律

运算中更新

查找表方式

蝶型运算规律

减小运算量：多类蝶形单元运算

代码实现

#include"fft.h"

int main()
{
input_xn();
FFT(xn, N);
print(xk, N);
delete[]xk;
return 0;
}

void print(complex<double> *x, int num)
{
for (int i = 0; i < num; i++)
{
//cout << (*(x + i)) << " ";
cout << (*(x + i)) << endl;
}
//cout << endl;
}
void print(int *x, int num)
{
for (int i = 0; i < num; i++)
{
cout << *(x + i) << " ";
}
cout << endl;
}
void input_xn()
{
cout << "请输入序列个数（2的整数幂次，例如：1024）：";
cin >> n;
xn = new complex<double>
;
cout << "请输入序列：";
for (int i = 0; i < n; i++)
{
cin >> *(xn + i);
}
cout << "您要做几点FFT？";
cin >> N;
}
void cal_input(complex<double> *x)
{
int seq_index, inv_index, cur_index_weight;
complex<double> temp;
inv_index = 0;

for (seq_index = 0; seq_index < N - 1; seq_index++)
{
//判断倒序后序列号大小，只需交换一半数值，且相等时不必交换
if (seq_index < inv_index)
{
temp = *(x + seq_index);
*(x + seq_index) = *(x + inv_index);
*(x + inv_index) = temp;
}

cur_index_weight = N >> 1;  //当前权值为最高位权值

while (cur_index_weight <= inv_index) //清空次高位，产生低位进位权值
{
inv_index = inv_index - cur_index_weight;
cur_index_weight >>= 1;
}

inv_index = inv_index + cur_index_weight;	//低位进1
}

}

void butterfly(complex<double> *x)
{
for (int L = 1; L <= M; ++L)
{
int B = (int)pow(2, L - 1);
for (int J = 0; J <= B - 1; ++J)
{
int P = (int)pow(2, M - L)*J;
complex<double> Wn_P(cos(2 * pi *P / N), sin(-2 * pi * P / N));
for (int k = J; k <= N - 1; k += (int)pow(2, L))
{
complex<double> tmp = *(x + k) + (*(x + k + B)*Wn_P);
*(x + k + B) = *(x + k) - (*(x + k + B)*Wn_P);
*(x + k) = tmp;
}
}
}
}

void FFT(complex<double> *x, int N)
{
xk = new complex<double>
;
for (int i = 0; i < N; i++)
{
*(xk + i) = *(x + i);
if (i >= n)
*(xk + i) = 0;
}
delete[]xn;//释放xn空间
M = int(log2(N));//计算级数
cal_input(xk);
butterfly(xk);

}