Discrete Fourier Transform

이산 푸리에 변환이란 주기가 \(N\)인 복소수열 \(x_{0},x_{1},\dots,x_{N-1}\)에 대하여 다음과 같이 정의되는 주기가 \(N\)인 복소수열 \(X_{0},X_{1},\dots,X_{N-1}\)로 바꾸는 변환을 말한다.

\[X_{k}=\sum_{n=0}^{N-1}e^{-2\pi ikn/N}x_{n}\]

여기서 \(e^{iz}=\cos z+i\sin z\)이며, 이로부터 \(e^{2\pi i}=e^{-2\pi i}=1\)임을 쉽게 알 수 있다.

~

이산 푸리에 역변환

또한 이렇게 변환된 복소수열 \(X_{k}\)는 다음과 같은 과정을 통하여 복소수열 \(x_{n}\)으로 역변환될 수 있다.

\[x_{n}=\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}\]

~

이산 푸리에 역변환의 증명

\(X_{k}\)의 정의대로 식을 풀어보면 다음과 같다.

\[\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}=\sum_{k=0}^{N-1}\sum_{j=0}^{N-1}e^{2\pi ik(n-j)/N}x_{j}=\sum_{j=0}^{N-1}x_{j}\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}\]

\(j=n\)인 경우에는 다음처럼 된다.

\[\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}=\sum_{k=0}^{N-1}1=N\]

\(j\neq n\)인 경우에는 다음처럼 된다.

\[\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}=\frac{e^{2\pi i(n-j)}-1}{e^{2\pi i(n-j)/N}-1}=0\]

따라서 다음이 성립한다.

\[\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}=\frac{1}{N}Nx_{n}=x_{n}\]

~

이산 푸리에 변환과 Convolution

주기가 \(N\)인 수열 \(a_{n}\)과 \(b_{n}\)을 각각 이산 푸리에 변환한 수열을 \(A_{n}\)과 \(B_{n}\)이라 하자.

또한 새로운 수열 \(C_{n}=A_{n}B_{n}\)을 정의하고, 이를 이산 푸리에 역변환한 수열을 \(c_{n}\)이라 하면 다음과 같은 식이 성립한다.

\[c_{n}=\sum_{k=0}^{N-1}a_{k}b_{n-k}\]

~

이산 푸리에 변환과 Convolution의 증명

\(C_{n}\)의 정의대로 식을 풀어보면 다음과 같다.

\[C_{n}=A_{n}B_{n}=\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}e^{-2\pi i(j+k)n/N}a_{k}b_{j}=\sum_{j=0}^{N-1}e^{-2\pi ijn/N}\sum_{k=0}^{N-1}a_{k}b_{j-k}\]

따라서 이를 역변환한 수열 \(c_{n}\)에 대한 식이 성립한다.

~

고속 푸리에 변환

이산 푸리에 변환을 단순한 반복문으로 계산하면 \(O(N^{2})\)의 시간복잡도를 가지지만, 이를 \(O(N\log N)\)의 시간복잡도에 행할 수 있는 방법이 고속 푸리에 변환이다.

고속 푸리에 변환에는 다양한 알고리즘이 있지만, 여기서는 \(N=2^{M}\) 꼴의 푸리에 변환에 최적화된 Cooley-Tukey 알고리즘에 대해 설명한다.

\(X_{k}\)를 다음처럼 고쳐 적을 수 있다.

\[X_{k}=\sum_{n=0}^{N-1}e^{-2\pi ikn/N}x_{n}=\sum_{n=0}^{N/2-1}e^{-2\pi ik(2n)/N}x_{2n}+\sum_{n=0}^{N/2-1}e^{-2\pi ik(2n+1)/N}x_{2n+1}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n}+e^{-2\pi ik/N}\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n+1}\]

이제 수열 \(Y_{k}\)와 \(Z_{k}\)를 다음처럼 정의하자.

\[Y_{k}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n}\]

\[Z_{k}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n+1}\]

수열 \(Y_{k}\)와 \(Z_{k}\)를 계산할 수 있다면, 수열 \(Z_{k}\)는 다음과 같은 방법으로 \(O(N)\)의 시간복잡도로 계산할 수 있다.

\[X_{k}=Y_{k}+e^{-2\pi ik/N}Z_{k}\]

\[X_{k+N/2}=Y_{k}-e^{-2\pi ik/N}Z_{k}\]

그런데, 수열 \(Y_{k}\)와 \(Z_{k}\)는 수열 \(x_{2n}\)과 \(x_{2n+1}\)에 대한 이산 푸리에 변환이므로, 분할정복을 통해 총 \(O(N\log N)\)의 시간복잡도로 이산 푸리에 변환을 구할 수 있다.

#include <cstdio>
#include <cmath>
#include <complex>

const int SZ = 20, N = 1 << SZ;
const double PI = acos(-1);

using namespace std;

int Rev(int x) {
	int i, r = 0;
	for (i = 0; i < SZ; i++) {
		r = r << 1 | x & 1;
		x >>= 1;
	}
	return r;
}

void FFT(complex<double> *a, bool f) {
	complex<double> x, y, z;
	double w;
	int i, j, k;
	for (i = 0; i < N; i++) {
		j = Rev(i);
		if (i < j) {
			z = a[i];
			a[i] = a[j];
			a[j] = z;
		}
	}
	for (i = 1; i < N; i <<= 1) {
		w = PI / i;
		if (f) w = -w;
		x = complex<double>(cos(w), sin(w));
		for (j = 0; j < N; j += i << 1) {
			y = complex<double>(1);
			for (k = 0; k < i; k++) {
				z = a[i + j + k] * y;
				a[i + j + k] = a[j + k] - z;
				a[j + k] += z;
				y *= x;
			}
		}
	}
	if (f) for (i = 0; i < N; i++) a[i] /= N;
}

complex<double> X[N];

int main() {
	double x;
	int i, n;
	scanf("%d", &n);
	for (i = 0; i <= n; i++) {
		scanf("%lf", &x);
		X[i] = complex<double>(x, 0);
	}
	FFT(X, false);
	for (i = 0; i < N; i++) X[i] *= X[i];
	FFT(X, true);
	for (i = 0; i <= n + n; i++) printf("%.0f ", X[i].real());
}

~

일반화된 이산 푸리에 변환

지금까지는 함수 \(\mathbb{Z}_{N}\rightarrow\mathbb{C}\)에 대한 이산 푸리에 변환에 대해서만 살펴보았다.

한편 이산 푸리에 변환에 필요한 연산은 덧셈과 복소수곱셈뿐이므로 이 둘이 정의된 모든 집합에 대해서 이산 푸리에 변환이 가능하다.

예를 들어, 함수 \(\mathbb{Z}_{M}\rightarrow\mathbb{C}\)는 덧셈과 복소수곱셈이 가능하므로 \(\mathbb{Z}_{N}\rightarrow\mathbb{Z}_{M}\rightarrow\mathbb{C}\)에 대한 이산 푸리에 변환이 가능하다.

이를 통해 두 \(N\times M\) 행렬 \(A_{i,j}\)와 \(B_{i,j}\)에 대하여 다음을 만족하는 행렬 \(C_{i,j}\)를 구할 수 있다.

\[C_{i,j}=\sum_{k=0}^{N-1}\sum_{l=0}^{M-1}A_{k,l}B_{i-k,j-l}\]

template을 사용하여 구현하면 다음과 같다.

#include <cstdio>
#include <cmath>
#include <complex>

const int SZ = 10, N = 1 << SZ;
const double PI = acos(-1);

using namespace std;

int Rev(int x) {
	int i, r = 0;
	for (i = 0; i < SZ; i++) {
		r = r << 1 | x & 1;
		x >>= 1;
	}
	return r;
}

template <typename T>
struct Func {
	T a[N];

	T &operator [](int x) {
		return a[x];
	}

	Func<T> operator +(Func<T> &x) {
		Func<T> r;
		int i;
		for (i = 0; i < N; i++) r.a[i] = a[i] + x[i];
		return r;
	}

	Func<T> operator -(Func<T> &x) {
		Func<T> r;
		int i;
		for (i = 0; i < N; i++) r.a[i] = a[i] - x[i];
		return r;
	}

	Func<T> operator *(complex<double> &x) {
		Func<T> r;
		int i;
		for (i = 0; i < N; i++) r.a[i] = a[i] * x;
		return r;
	}

	Func<T> operator *=(complex<double> &x) {
		int i;
		for (i = 0; i < N; i++) a[i] *= x;
		return *this;
	}

	Func<T> operator *=(Func<T> &x) {
		int i;
		for (i = 0; i < N; i++) a[i] *= x[i];
		return *this;
	}

	void FFT(bool f) {
		complex<double> x, y;
		T z;
		double w;
		int i, j, k;
		for (i = 0; i < N; i++) {
			j = Rev(i);
			if (i < j) {
				z = a[i];
				a[i] = a[j];
				a[j] = z;
			}
		}
		for (i = 1; i < N; i <<= 1) {
			w = PI / i;
			if (f) w = -w;
			x = complex<double>(cos(w), sin(w));
			for (j = 0; j < N; j += i << 1) {
				y = complex<double>(1);
				for (k = 0; k < i; k++) {
					z = a[i + j + k] * y;
					a[i + j + k] = a[j + k] - z;
					a[j + k] = a[j + k] + z;
					y *= x;
				}
			}
		}
		if (f) {
			x = complex<double>((double)1 / N);
			for (i = 0; i < N; i++) a[i] *= x;
		}
	}
};

Func<Func<complex<double> > > X;

int main() {
	double x;
	int i, j, n, m;
	scanf("%d%d", &n, &m);
	for (i = 0; i <= n; i++) for (j = 0; j <= m; j++) {
		scanf("%lf", &x);
		X[i][j] = complex<double>(x, 0);
	}
	for (j = 0; j < N; j++) X[j].FFT(false);
	X.FFT(false);
	X *= X;
	X.FFT(true);
	for (j = 0; j < N; j++) X[j].FFT(true);
	for (i = 0; i <= n + n; i++) {
		for (j = 0; j <= m + m; j++) printf("%.0f ", X[i][j].real());
		printf("\n");
	}
}

~

XOR Convolution

\(2^{M}\) 크기의 배열 \(A_{i}\)와 \(B_{i}\)에 대하여 다음처럼 정의되는 배열 \(C_{i}\)를 구하는 문제를 생각해 보자. (\(\oplus\)는 Bitwise XOR 연산을 뜻함)

\[C_{i}=\sum_{k=0}^{2^{M}-1}A_{k}B_{i\oplus k}\]

XOR 연산은 각 자리별로 더한 뒤 2로 나눈 나머지를 구하는 것과 같으므로, \(\mathbb{Z}_{2^{M}}\rightarrow\mathbb{C}\)를 다음처럼 펼쳐서 생각해 보자.

\[\mathbb{Z}_{2}\rightarrow\mathbb{Z}_{2}\rightarrow\cdots\rightarrow\mathbb{Z}_{2}\rightarrow\mathbb{C}\]

각각이 이산 푸리에 변환 가능하므로, \(\mathbb{Z}_{2}\)에 대한 이산 푸리에 변환을 \(M\)번 행하는 것임을 알 수 있다.

#include <cstdio>
#include <complex>

const int SZ = 20, N = 1 << SZ;

using namespace std;

int Rev(int x) {
	int i, r = 0;
	for (i = 0; i < SZ; i++) {
		r = r << 1 | x & 1;
		x >>= 1;
	}
	return r;
}

void FFT(int *a, bool f) {
	int i, j, k, z;
	for (i = 0; i < N; i++) {
		j = Rev(i);
		if (i < j) {
			z = a[i];
			a[i] = a[j];
			a[j] = z;
		}
	}
	for (i = 1; i < N; i <<= 1) for (j = 0; j < N; j += i << 1) for (k = 0; k < i; k++) {
		z = a[i + j + k];
		a[i + j + k] = a[j + k] - z;
		a[j + k] += z;
	}
	if (f) for (i = 0; i < N; i++) a[i] /= N;
}

int X[N];

int main() {
	int i, n;
	scanf("%d", &n);
	for (i = 0; i < 1 << n; i++) scanf("%d", &X[i]);
	FFT(X, false);
	for (i = 0; i < N; i++) X[i] *= X[i];
	FFT(X, true);
	for (i = 0; i < 1 << n; i++) printf("%d ", X[i]);
}