目的

画像処理100本ノックに挑戦中です。Q.27はバイキュービック補間ですが、どうもサンプリング定理の応用らしいです。ここらへんをちゃんと理解するためにサンプリング定理を復習します。

サンプリング定理

観測対象は連続実関数$f(x)$です。離散サンプリングは、コム関数で表現されます。

s(x) = \sum_n \delta(x-nL)

$L$はサンプリング間隔です。サンプリング後の信号は

\tilde f(x) = f(x)\sum_n \delta(x-nL)

とあらわされます。

これをフーリエ変換します。畳み込みの定理より

\mathcal F[\tilde f](\nu) = \mathcal F[f] \otimes \mathcal F[\sum_n\delta(x-nL)]

となります。コム関数のフーリエ変換は

\mathcal F[\sum_n\delta(x-nL)] = \frac{1}{L}\sum_n\delta(\nu- n/L)

です。これから、

\mathcal F[\tilde f](\nu) = \frac{1}{L}\sum_n F(\nu- n/L)

となります。元の連続信号$f$のフーリエ変換$F(\nu)$を$1/L$周期で並べたものになるわけです。サンプリング間隔を小さくすれば、この周期は長くなります。つまり、十分に小さなサンプリング間隔をとれば、隣り合った$F(\nu)$同士が重なり合わなくなります。実関数のフーリエ変換が偶関数となることに注意すれば、

\nu_{\rm max} < \frac{1}{2L}

を満たすようにサンプリングすればよいわけです。これがサンプリング定理です。信号が含む最大周波数の2倍以上のサンプリング周波数が必要ということです。2倍なのは実関数のフーリエ変換が偶関数であることに由来するわけです。

さて、この条件を満たしている場合、隣り合う$F$同士は重なり合っていないわけですから、周波数空間上で幅$1/L$の矩形関数($\times L$)をかけてやれば元の$F(\nu)$が得られるわけです。実空間で言えば、

\mathcal F[L\times {\rm rect}_{1/L}] = L\times \frac{\sin(\pi x/L)}{\pi x}

より、sinc関数を離散サンプリングされた信号$\tilde f$に畳み込んでやれば元の連続信号$f$が復元出来ることになります。このsinc関数の周期は$2L$ですから、畳み込んでも、サンプリングした点において寄与するのはその点のみであり値が変わらないことがわかります。うまく出来ています。

電気信号などであれば、サンプリングされた離散パルス列$\tilde f$に理想ローパスフィルタをかけてやれば連続信号$f$が生成出来るわけです。
画像処理のように計算機上でデジタル的な計算をする場合は、実空間での畳み込みによって任意の位置$x$におけるアナログ信号値が得られます。

離散信号から連続信号を復元してみる

元の信号を$f(x) = \exp(-\frac{x^2}{2\sigma^2})$とします。サンプリング間隔は$L=1$とします。

# include <iostream>
# include <vector>
const double pi = 3.14159265358979323846;

int main()
{
	int N = 1024;
	auto gauss = [&](double x)
	{
		double m = N / 2;
		double s = 100;
		return exp(-(x - m) * (x - m)/2./s/s);
	};
	std::vector<double> data;
	for (int i = 0; i < N; i++)
	{
		data.push_back(gauss(i));
	}

	auto f = [&](double x)
	{
		double ret = 0;
		for (int i = 0; i < N; i++)
		{
			double d = x - i;
			if (std::abs(d) > 1e-5)ret += data[i] * sin(pi * d) / pi / d;
			else ret += data[i];
		}
		return ret;
	};

	for (double x = N/2-10; x < N/2+10; x += 0.2)
	{
		std::cout << x << " " << f(x) << " " << gauss(x) << std::endl;
	}


	return 0;
}

実線が真の値、点が復元した値です。整数点しか測定していないのにちゃんと復元出来ています。

バイキュービック補間との関係

バイキュービック補間はsinc関数を3次近似した式を用いるらしいです。補間というよりは"回復"と言った方がしっくりきます。もちろん、単一画素内での急激な変化がない前提でないと正しく回復出来るわけではないのですが。
詳しい計算は今後書く予定の【画像処理100本ノックに挑戦】記事中でしてみるつもりです。