らんだむな記憶

blogというものを体験してみようか!的なー

Fourier変換~理想と現実

$d$次元Fourier変換を
\begin{equation}
\mathcal{F}u(\xi) = \frac{1}{(2\pi)^{d/2}}\int_{\R^d}\exp(- i x \cdot \xi) u(x) dx
\end{equation}

で定義する。
ここでは、2次元Fourier変換を考え、$u(x) = \chi_{[-1,1]^2}(x)$ つまり、
\begin{equation}
\chi_{[-1,1]^2}(x) =
\begin{cases}
1, \hspace{1em} x_1, x_2 \in [-1,1], \\
0, \hspace{1em} \text{otherwise}
\end{cases}
\end{equation}

とする。$u$ のFourier変換はFubiniの定理により、
\begin{equation}
\fbox{$\mathcal{F}\chi_{[-1,1]^2}(\xi) = \frac{2}{\pi} \left( \frac{\sin \xi_1}{\xi_1} \right)\left( \frac{\sin \xi_2}{\xi_2} \right)$}
\end{equation}

と求まる。理論上は以下のような性質になっている:
$\chi_{[-1,1]^2} \in L^1(\R^2) \cap L^2(\R^2)$ なので $\mathcal{F}\chi_{[-1,1]^2} \in L^\infty(\R^2) \cap L^2(\R^2)$


見えにくいが、以下の画像において左の画像は原点付近に3x3くらいの大きさの白い箱を置いている。Fourier変換の積分核(三角函数)の周期くらいになってほしいなーという気持ちの大きさだが、どうだろうな。要するに、原点付近では、値255をとり、それ以外では0をとるような定義函数を用意した形だ。これをFFTをかますと理論上は上記のような結果に近いものになるはず。但し、周波数空間で係数は正と負に散らばっている可能性があって、画像に落とし込むには値の範囲を[0,255]に収め込まないとならない。ということで、周波数の絶対値をとったマグニチュードスペクトルにしている。パワースペクトルだと2乗の作用で小さいところがとことん小さくなるかもしれないのでやめた。
セッティングの妥当性はよく分からないが、縦軸方向と横軸方向の$\left|\frac{\sin \xi}{\xi}\right|$が掛け合わされたような形状になっているのではないだろうか?減衰が一番少ない原点付近が値が一番大きいので白味が強くなっている。減衰の影響で分かりにくいが、5x5 = 25個の白の長方形領域が見えるはずである。
f:id:derwind:20160130210657p:plain

これには以下のようなコードを使った。対数をとらなくても十分だったのでとっていない。

f = np.fft.fft2(img)
f = np.fft.fftshift(f)
raw_mag = np.abs(f)

magnitude_spectrum = raw_mag

Amazon.co.jp: ディジタル画像処理[改訂新版]: 画像情報教育振興協会: 本という本は入門書(まったくの1冊目に適当かは分からないんだが)にしては話題も広く浅くという感じで、ラベリングやらFourier変換やらパターン認識までまぁ色々載ってはいる。Fourier変換の話題は「6-1 画像のフーリエ変換」が該当し、周期は上記の設定とは異なるのだが、ほぼ同等の結果が p.127 に乗っている。よくありがち(?)な CHAR_MAX に等しい番号のページと言っても良いかもしれない。

さて、ついでなので、プロットもしよう。これもpythonを使っておく。

#! /usr/bin/env python
# -*- coding: utf-8 -*-

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x = np.arange(-20, 20, 0.5)
y = np.arange(-20, 20, 0.5)
xx, yy = np.meshgrid(x, y)
zz = (2/math.pi) * (np.sin(xx)/(xx)) * (np.sin(yy)/(yy))

figure = plt.figure()
ax = Axes3D(figure)
ax.plot_wireframe(xx, yy, zz)
plt.show()

結果は以下。マグニチュードスペクトルのやつと組み合わせて理解するためには多少はイメージ力が必要かもしれない。
f:id:derwind:20160130215137p:plain