【画像処理100本ノックに挑戦】独自の離散フーリエ変換クラスを作る

背景

画像処理100本ノックに挑戦中です。Q.32からしばらくフーリエ変換関連の問題です。見通しよく回答するためには簡単にクラスの形でまとめておいた方がよさそうなのでここでやっちゃいます。ちなみに普段はFFTWを使っています。

簡単な解説

高速フーリエ変換(FFT)ではないので素直に定義式を実装するだけです。1次元離散フーリエ変換(DFT)が実装できれば、2次元では$x$を変数としてDFTした後に、$y$を変数としてDFTすればOKです。

面倒なのはDFTシフトです。定義式通り実装すると、フーリエ空間での添え字を$n$として$n=0$が直流成分なのですが、しばしば$n=\lfloor {N/2} \rfloor$を直流成分とした方が便利です。このような座標の平行移動がDFTシフトです。
離散フーリエ変換では周期性があるため、全体的に右に移動させてはみ出たら左から挿入していけばよいです。結局これは、一次元配列を真ん中で二つに分割し、入れ替える処理に相当します。ただし$N$が奇数の時、真ん中の要素は前半分の配列に含まれると考えます。
逆変換(DFT逆シフト)には注意が必要です。$N$が偶数の場合は順シフトと逆シフトは同じですが、$N$が奇数の場合は異なります。真ん中の要素は後ろ半分の配列に含まれるとみなすことに注意が必要です。2次元の場合は$x$方向のシフト後、$y$方向にシフトするだけです。

コード

冗長な部分が多いので今後修正するかもしれません。

const double pi = 3.14159265358979323846;
const std::complex<double> imag(0, 1);
class DFT
{
public:
	DFT() {}
	std::vector<std::complex<double>> forward(const std::vector<std::complex<double>> &f)
	{
		int N = f.size();
		std::vector<std::complex<double>> F(N);
		for (int t = 0; t < N; t++)
		{
			F[t] = 0;
			for (int x = 0; x < N; x++)
			{
				F[t] += f[x] * exp(-imag * 2. * pi * (double)t * (double)x / (double)N);
			}
		}
		return F;
	}
	std::vector<std::complex<double>> backward(const std::vector<std::complex<double>>& F)
	{
		int N = F.size();
		std::vector<std::complex<double>> f(N);
		for (int t = 0; t < N; t++)
		{
			f[t] = 0;
			for (int x = 0; x < N; x++)
			{
				f[t] += F[x] * exp(imag * 2. * pi * (double)t * (double)x / (double)N) / (double)N;
			}
		}
		return f;
	}
	std::vector<std::complex<double>> shift(const std::vector<std::complex<double>>& F)
	{
		int N = F.size();
		std::vector<std::complex<double>> ret(N);
		if (N % 2 == 0)
		{
			for (int n = 0; n < N / 2; n++)
			{
				ret[n + N / 2] = F[n];
			}
			for (int n = N / 2; n < N; n++)
			{
				ret[n - N / 2] = F[n];
			}
		}
		else
		{
			for (int n = 0; n <= N / 2; n++)
			{
				ret[n + N / 2] = F[n];
			}
			for (int n = N / 2 + 1; n < N; n++)
			{
				ret[n - N / 2 - 1] = F[n];
			}
		}
		return ret;
	}
	std::vector<std::complex<double>> ishift(const std::vector<std::complex<double>>& F)
	{
		int N = F.size();
		std::vector<std::complex<double>> ret(N);
		if (N % 2 == 0)
		{
			ret = shift(F);
		}
		else
		{
			for (int n = N / 2; n < N; n++) ret[n - N / 2] = F[n];
			for (int n = 0; n < N / 2; n++) ret[n + N / 2 + 1] = F[n];
		}
		return ret;
	}
};

class DFT2D
{
private:
	int width, height;
public:
	DFT2D(int width, int height) : width{ width }, height{height}
	{
	}
	std::vector<std::vector<std::complex<double>>> forward(const std::vector<std::vector<std::complex<double>>> &f)
	{
		DFT dft;
		std::vector<std::complex<double>> fx(width);
		std::vector<std::vector<std::complex<double>>> ret(width, std::vector<std::complex<double>>(height));
		for (int j = 0; j < height; j++)
		{
			for (int i = 0; i < width; i++)
			{
				fx[i] = f[i][j];
			}
			fx = dft.forward(fx);//xを変数としてdft
			for (int i = 0; i < width; i++)
			{
				ret[i][j] = fx[i];
			}
		}
		for (int i = 0; i < width; i++)
		{
			ret[i] = dft.forward(ret[i]);//yを変数として
		}
		return ret;
	}
	std::vector<std::vector<std::complex<double>>> backward(const std::vector<std::vector<std::complex<double>>>& f)
	{
		DFT dft;
		std::vector<std::complex<double>> fx(width);
		std::vector<std::vector<std::complex<double>>> ret(width, std::vector<std::complex<double>>(height));
		for (int j = 0; j < height; j++)
		{
			for (int i = 0; i < width; i++)
			{
				fx[i] = f[i][j];
			}
			fx = dft.backward(fx);//xを変数としてdft
			for (int i = 0; i < width; i++)
			{
				ret[i][j] = fx[i];
			}
		}
		for (int i = 0; i < width; i++)
		{
			ret[i] = dft.backward(ret[i]);//yを変数として
		}
		return ret;
	}
	std::vector<std::vector<std::complex<double>>> shift(const std::vector<std::vector<std::complex<double>>>& f)
	{
		DFT dft;
		std::vector<std::complex<double>> fx(width);
		std::vector<std::vector<std::complex<double>>> ret(width, std::vector<std::complex<double>>(height));
		for (int j = 0; j < height; j++)
		{
			for (int i = 0; i < width; i++)
			{
				fx[i] = f[i][j];
			}
			fx = dft.shift(fx);//xを変数として
			for (int i = 0; i < width; i++)
			{
				ret[i][j] = fx[i];
			}
		}
		for (int i = 0; i < width; i++)
		{
			ret[i] = dft.shift(ret[i]);//yを変数として
		}
		return ret;
	}
	std::vector<std::vector<std::complex<double>>> ishift(const std::vector<std::vector<std::complex<double>>>& f)
	{
		DFT dft;
		std::vector<std::complex<double>> fx(width);
		std::vector<std::vector<std::complex<double>>> ret(width, std::vector<std::complex<double>>(height));
		for (int j = 0; j < height; j++)
		{
			for (int i = 0; i < width; i++)
			{
				fx[i] = f[i][j];
			}
			fx = dft.ishift(fx);//xを変数として
			for (int i = 0; i < width; i++)
			{
				ret[i][j] = fx[i];
			}
		}
		for (int i = 0; i < width; i++)
		{
			ret[i] = dft.ishift(ret[i]);//yを変数として
		}
		return ret;
	}

};

使い方 & 動作確認

1DDFT:

int main()
{
	int N = 128;
	std::vector<std::complex<double>> f(N);
	for (int n = 0; n < N; n++)
	{
		f[n] = cos(2 * pi * n * 10 / N);
	}
	DFT dft;

	auto F = dft.forward(f);
	for (int n = 0; n < N; n++)
	{
		std::cout << F[n].real() << " " << F[n].imag() << std::endl;
	}
	return 0;
}

1DDFT + DFT Shift:

int main()
{
	int N = 128;
	std::vector<std::complex<double>> f(N);
	for (int n = 0; n < N; n++)
	{
		f[n] = cos(2 * pi * n * 10 / N);
	}
	DFT dft;

	auto F = dft.forward(f);
	F = dft.shift(F);
	for (int n = 0; n < N; n++)
	{
		std::cout << F[n].real() << " " << F[n].imag() << std::endl;
	}
	return 0;
}