作為學習過程，我們盡量少用 庫函數(shù)。由于我們需要進行復數(shù)乘法梭纹、加法運算，下面先粗糙地定義一下復數(shù)運算：
template<typename T>
struct Complex {
T r;
T i;
Complex(): r((T)0), i((T)0) { }
Complex(T real, T image): r(real), i(image) { }
Complex(const Complex<T>& other) {
r = other.r;
i = other.i;
}
Complex<T>& operator = (const Complex<T>& c) {
if (&c == this) return *this;
r = c.r;
i = c.i;
return *this;
}
};
template<typename T>
Complex<T> operator + (const Complex<T>& c1, const Complex<T>& c2) {
Complex<T> x;
x.r = c1.r + c2.r;
x.i = c1.i + c2.i;
return x;
}
template<typename T>
Complex<T> operator - (const Complex<T>& c1, const Complex<T>& c2) {
Complex<T> x;
x.r = c1.r - c2.r;
x.i = c1.i - c2.i;
return x;
}
template<typename T>
Complex<T> operator * (const Complex<T>& c1, const Complex<T>& c2) {
Complex<T> x;
x.r = c1.r * c2.r - c1.i * c2.i;
x.i = c1.r * c2.i + c1.i * c2.r;
return x;
}
我們考慮使用靜態(tài)大小的 FFT 樣本點數(shù) N致份，也就是在編譯期变抽，我們已經知道 FFT 運算的樣本點數(shù)了。因此我們將該參數(shù)直接作為 C++ 模板參數(shù)進行傳遞氮块。
template <size_t N, typename T>
class naivefft {
public:
typedef Complex<T> cpx_t;
public:
naivefft() ;
void cdft(cpx_t* a) ;
private:
int order_of_2(int n);
void bit_reverse(cpx_t* a);
void butterfly_op(cpx_t* a, int l, int r, cpx_t tw) ;
private:
cpx_t _twiddle_factors[N/2+1];
int _stages;
int _N;
};
我們先來實現(xiàn)一下 蝶形運算绍载，為了減少復數(shù)乘法，我們用一個中間變量保存乘法結果：
template<size_t N, typename T>
void naivefft<N, T>::butterfly_op(cpx_t* a, int l, int r, cpx_t tw) {
// we get two FFT coefficients with only one multiplication and two additions.
cpx_t x = a[r] * tw;
cpx_t y = a[l] + x;
a[r] = a[l] - x;
a[l] = y;
}
進行 FFT 主體循環(huán)時滔蝉，最外層的循環(huán)是遍歷層數(shù)击儡，我們稱為 階段 stages；為此蝠引，我們需要一個計算 2 的整數(shù)次冪 的函數(shù)：
template<size_t N, typename T>
int naivefft<N, T>::order_of_2(int n) {
int ord = 0;
// counting bits followed by the 1
while (n > 1) {
ord += 1;
n >>= 1;
}
return ord;
}
不過在我們計算 order_of_2 時阳谍，我們要先確認參數(shù) n 是 2 的整數(shù)次冪，為此我們使用一個 編譯期間進行檢查 的元函數(shù)：
template <size_t N>
struct is_order_2 {
enum { yes = (N > 0) && ((N & (N - 1)) == 0) };
};
在構造函數(shù)中螃概，我們初始化 旋轉因子 需要一個函數(shù)矫夯，即 復數(shù)版 exp 函數(shù)，這是我們唯一調用 庫函數(shù) 的地方：
template<typename T>
Complex<T> exp(Complex<T> c) {
return Complex<T>(std::exp(c.r) * std::cos(c.i), std::exp(c.r) * std::sin(c.i));
}
現(xiàn)在完成我們的構造函數(shù)吊洼，初始化必要的 旋轉因子训貌；
template<size_t N, typename T>
naivefft<N,T>::naivefft() {
// checking the validity of FFT length at compile stage.
static_assert(is_order_2<N>::yes, "N must be power of 2.");
_N = N;
// the primary root of _Nth order of 1.
T phi = -2 * std::acos((T)-1) / _N;
// we only calculate the necessary roots of 1.
for (size_t i = 0; i < _N / 2 + 1; ++i)
_twiddle_factors[i] = exp(cpx_t(0, i * phi));
// calculate the order of 2 of N, and this indicates our fft stages.
_stages = order_of_2(_N);
}
到此，我們迎來了第一個重要的函數(shù)：按位逆序 (bit reversing)冒窍。
template<size_t N, typename T>
void naivefft<N,T>::bit_reverse(cpx_t* a) {
cpx_t x;
// traverse each index, since we know the
// 1st and last indices keep unchanged,
// we scan from 1 to _N-2.
for (int i = 1; i < _N - 1; ++i) {
int o = _stages; // valid bits length
int c = i; // temp value for index i
int ri = 0; // reversed index of i
// traverse each bit of i and
// construct the reversed index of i
// bit by bit.
while (o--) {
ri = (ri << 1) | (c & 1);
c = (c >> 1);
}
// swap the two numbers when ri < i
// to prevent swapping back.
if (ri < i) {
x = a[ri];
a[ri] = a[i];
a[i] = x;
}
}
}
最后递沪，完成我們的 Radix2-DIT-FFT 的主循環(huán)，具體過程可以對照著 我的 “基2-快速傅立葉變換” 的樸素實現(xiàn)（推導部分） 中的 圖 3）來看综液；事實上我寫代碼的時候畫了長度為 16 的流圖款慨。
template<size_t N, typename T>
void naivefft<N, T>::cdft(cpx_t* a) {
// first, permutating the numbers with bit reversed ordering.
bit_reverse(a);
// at each stage, we have several independent blocks.
int block_size = 2;
// the butterfly stride within the blocks.
int stride = 1;
// temp val for calculating the stride
// scanning at the twiddle factors table.
int tw_stride_temp = 1;
// taversing each stage.
for (int stage = 0; stage < _stages; stage++) {
// at new stage, we need to shrinking the twiddle stride,
// so we double the denominator.
tw_stride_temp *= 2;
int tw_stride = _N / tw_stride_temp;
// traversing each block of current stage.
for (int b = 0; b < _N; b += block_size) {
// calculating all the butterfly operations within one block.
for (int s = 0; s + stride < block_size; s++) {
butterfly_op(a, b + s, b + s + stride, _twiddle_factors[s * tw_stride]);
}
}
// according the flow graph, we need to double
// the butterfly stride and block size at the next stage.
stride *= 2;
block_size *= 2;
}
}
