【離散コサイン変換】DFTを使ってDCTしてみる

背景

【離散コサイン変換】DCT-IIの式をぼんやり理解するの続きです。

使用したライブラリ

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

DFTを使ってDCT-IIしてみる

$N$点のデータ$x_n$があったとして、$4N$点の$y_n$を

\begin{eqnarray}
y_{2n+1} &=& y_{4N-2n-1}=x_n\;(0\leq n < N)\\
\end{eqnarray}

で定義してやってDFTすればOKです。上の点以外はゼロとします。

std::vector<double> dct(std::vector<double> in)
{
	int N = in.size();
	std::vector<double> out(N);
	for (int k = 0; k < N; k++)
	{
		out[k] = 0;
		for (int n = 0; n < N; n++)
		{
			out[k] += in[n] * cos(pi / N * (n + 0.5) * k);
		}
	}
	return out;
}

int main()
{
	int N = 128;
	std::vector<double> data(N), out(N);
	for (int i = 0; i < N; i++)
	{
		data[i] = cos(2 * pi * i / 10);
	}

	out = dct(data);

	std::vector<std::complex<double>> data2(4 * N);
	for (int i = 0; i < 4 * N; i++) data2[i] = 0;
	for (int i = 0; i < N; i++) data2[2 * i + 1] = data2[4 * N - 2 * i - 1] = data[i] / 2.;

	DFT dft;
	std::vector<std::complex<double>> out2 = dft.forward(data2);

	for (int i = 0; i < N; i++)
	{
		std::cout << out[i] << " " << out2[i].real() << std::endl;
	}

	return 0;
}

定義どおりにDCTした場合と結果が一致することが確認できました。

DFTを使ってDCT-IIIしてみる

$N$点のデータ$x_n$があったとして、$4N$点の$y_n$を

\begin{eqnarray}
0\leq n<N&\to& y_n = x_n\\
n=N&\to& y_n = 0\\
N<n\leq 2N &\to& y_n = -x_{2N-n}\\
2N<n<3N &\to& y_n = -x_{n-2N}\\
n=3N &\to& y_n = 0\\
3N<n<4N &\to& y_n = x_{4N-n}
\end{eqnarray}

で定義してやってDFTすればOKです。

std::vector<double> dct3(std::vector<double> in)
{
	int N = in.size();
	std::vector<double> out(N);
	for (int k = 0; k < N; k++)
	{
		out[k] = in[0] / 2.;
		for (int n = 1; n < N; n++)
		{
			out[k] += in[n] * cos(pi / N * (k + 0.5) * n);
		}
	}
	return out;
}

int main()
{
	int N = 128;
	std::vector<double> data(N), out(N);
	for (int i = 0; i < N; i++)
	{
		data[i] = cos(2 * pi * i / 10);
	}

	out = dct3(data);

	std::vector<std::complex<double>> data2(4 * N);
	for (int i = 0; i < 4 * N; i++) data2[i] = 0;
	for (int n = 0; n < N; n++) data2[n] = data[n];
	for (int n = N + 1; n <= 2 * N; n++) data2[n] = -data[2 * N - n];
	for (int n = 2 * N + 1; n < 3 * N; n++) data2[n] = -data[n - 2 * N];
	for (int n = 3 * N + 1; n < 4 * N; n++) data2[n] = data[4 * N - n];

	DFT dft;
	std::vector<std::complex<double>> out2 = dft.forward(data2);

	for (int i = 0; i < N; i++)
	{
		std::cout << out[i] << " " << out2[2 * i + 1].real() / 4. << std::endl;
	}

	return 0;
}

定義どおりにDCTした場合と結果が一致することが確認できました。

次回は二次元DCTについて少し考えてみるつもりです。