快速傅立葉變換FFT(模板)
阿新 • • 發佈:2018-11-08
轉載出處 https://blog.csdn.net/f_zyj/article/details/76037583
摘自大佬的部落格 FFT(最詳細最通俗的入門手冊)
const double PI=acos(-1.0); // 複數結構體 struct Complex { double x,y; // 實部和虛部 x + yi Complex(double _x=0.0, double _y = 0.0) { x = _x; y = _y; } Complex operator - (const Complex &b) const { return Complex(x - b.x, y - b.y); } Complex operator + (const Complex &b) const { return Complex(x + b.x, y + b.y); } Complex operator * (const Complex &b) const { return Complex(x * b.x - y * b.y, x * b.y + y * b.x); } }; // 進行FFT和IFFT前的反轉變換 // 位置i和(i二進位制反轉後的位置)互換 // len必須去2的冪 void change(Complex y[], int len) { int i, j, k; for (i = 1, j = len / 2; i < len - 1; i++) { if (i < j) { swap(y[i], y[j]); } // 交換護衛小標反轉的元素,i < j保證交換一次 // i做正常的+1,j左反轉型別的+1,始終保持i和j是反轉的 k = len / 2; while (j >= k) { j -= k; k /= 2; } if (j < k) { j += k; } } return ; } // FFT // len必須為2 ^ k形式 // on == 1時是DFT,on == -1時是IDFT void fft(Complex y[], int len, int on) { change(y, len); for (int h = 2; h <= len; h <<= 1) { Complex wn(cos(-on * 2 * PI / h), sin(-on * 2 * PI / h)); for (int j = 0; j < len; j += h) { Complex w(1, 0); for (int k = j; k < j + h / 2; k++) { Complex u = y[k]; Complex t = w * y[k + h / 2]; y[k] = u + t; y[k + h / 2] = u - t; w = w * wn; } } } if (on == -1) { for (int i = 0; i < len; i++) { y[i].x /= len; } } }