1. 程式人生 > >快速傅立葉變換FFT(模板)

快速傅立葉變換FFT(模板)

轉載出處 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;
		}
	}
}