scipy.fftpackのコンベンション(2次元DFT)

1次元の場合に引き続き、2次元配列の場合の動作を確認する。

● 2次元配列の1次元DFT
2次元DFTに入る前に、1次元DFT、つまりscipy.fftpack.fft()の2次元配列に対する動作を見よう。

>>> from scipy.fftpack import fft
>>> a = [[1]]
>>> fft(a)
array([[ 1.+0.j]])
>>> a = [[1, 0], [0, 1]]
>>> fft(a)
array([[ 1.+0.j, 1.+0.j],
[ 1.+0.j, -1.+0.j]])

これは、各成分に対して1次元DFTしている。

長さ3にすれば、

>>> a = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
>>> fft(a)
array([[ 1.0+0.j , 1.0+0.j , 1.0-0.j ],
[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j],
[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j]])


長さ4にすれば、

>>> a = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
>>> fft(a)
array([[ 1.+0.j, 1.+0.j, 1.+0.j, 1.-0.j],
[ 1.+0.j, 0.-1.j, -1.+0.j, 0.+1.j],
[ 1.+0.j, -1.+0.j, 1.+0.j, -1.-0.j],
[ 1.+0.j, 0.+1.j, -1.+0.j, 0.-1.j]])



この動作はscipy.fftpack.ifft()も同じ。

>>> a = [[1]]
>>> ifft(a)
array([[ 1.+0.j]])
>>> a = [[1, 0], [0, 1]]
>>> ifft(a)
array([[ 0.5+0.j, 0.5+0.j],
[ 0.5+0.j, -0.5+0.j]])
>>> a = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
>>> ifft(a)
array([[ 0.33333333+0.j , 0.33333333-0.j , 0.33333333+0.j ],
[ 0.33333333+0.j , -0.16666667+0.28867513j,
-0.16666667-0.28867513j],
[ 0.33333333+0.j , -0.16666667-0.28867513j,
-0.16666667+0.28867513j]])
>>> a = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
>>> ifft(a)
array([[ 0.25+0.j , 0.25-0.j , 0.25+0.j , 0.25+0.j ],
[ 0.25+0.j , 0.00+0.25j, -0.25+0.j , 0.00-0.25j],
[ 0.25+0.j , -0.25-0.j , 0.25+0.j , -0.25+0.j ],
[ 0.25+0.j , 0.00-0.25j, -0.25+0.j , 0.00+0.25j]])


1次元配列の場合に見た通り、、fft()とifft()の合成は恒等写像。

>>> a = [[1]]
>>> fft(ifft(a))
array([[ 1.+0.j]])
>>> ifft(fft(a))
array([[ 1.+0.j]])
>>> a = [[1, 0], [0, 1]]
>>> fft(ifft(a))
array([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, 1.+0.j]])
>>> ifft(fft(a))
array([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, 1.+0.j]])
>>> a = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
>>> fft(ifft(a))
array([[ 1.00000000e+00+0.j, 0.00000000e+00+0.j, 0.00000000e+00+0.j],
[ 0.00000000e+00+0.j, 1.00000000e+00+0.j, 5.55111512e-17+0.j],
[ 0.00000000e+00+0.j, 5.55111512e-17+0.j, 1.00000000e+00+0.j]])
>>> ifft(fft(a))
array([[ 1.00000000e+00+0.j, 0.00000000e+00+0.j, 0.00000000e+00+0.j],
[ 0.00000000e+00+0.j, 1.00000000e+00+0.j, 7.40148683e-17+0.j],
[ 0.00000000e+00+0.j, 7.40148683e-17+0.j, 1.00000000e+00+0.j]])
>>> a = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
>>> fft(ifft(a))
array([[ 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[ 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
[ 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
[ 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j]])
>>> ifft(fft(a))
array([[ 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[ 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
[ 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
[ 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j]])


だが、各成分の長さが異なると怒られる。

>>> a = [[1], [0, 1], [0, 0, 1], [0, 0, 0, 1]]
>>> fft(a)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/Library/Python/2.7/site-packages/scipy-0.11.0.dev_1983db6_20120208-py2.7-macosx-10.7-x86_64.egg/scipy/fftpack/basic.py", line 212, in fft
tmp = _asfarray(x)
File "/Library/Python/2.7/site-packages/scipy-0.11.0.dev_1983db6_20120208-py2.7-macosx-10.7-x86_64.egg/scipy/fftpack/basic.py", line 118, in _asfarray
return numpy.asfarray(x)
File "/Library/Python/2.7/site-packages/numpy-2.0.0.dev_4c0576f_20120208-py2.7-macosx-10.7-x86_64.egg/numpy/lib/type_check.py", line 103, in asfarray
return asarray(a,dtype=dtype)
File "/Library/Python/2.7/site-packages/numpy-2.0.0.dev_4c0576f_20120208-py2.7-macosx-10.7-x86_64.egg/numpy/core/numeric.py", line 333, in asarray
maskna=maskna, ownmaskna=ownmaskna)
ValueError: setting an array element with a sequence.


● 2次元配列の2次元DFT
これから本番の2次元DFT。

1×1の場合、

>>> from scipy.fftpack import fft2
>>> a = [[1]]
>>> fft2(a)
array([[ 1.+0.j]])

恒等写像だ。


>>> from scipy.fftpack import ifft2
>>> a = [[1]]
>>> ifft2(a)
array([[ 1.+0.j]])

逆も然り。
正確には、逆変換も然り、だけど。


2×2の場合、

>>> a = [[1, 0], [0, 1]]
>>> fft2(a)
array([[ 2.+0.j, 0.+0.j],
[ 0.+0.j, 2.+0.j]])


これは次と等しいはず。

>>> a = [[1, 0], [0, 1]]
>>> fft(fft(a).transpose()).transpose()
array([[ 2.+0.j, 0.+0.j],
[ 0.+0.j, 2.+0.j]])

確かにそうなっている。

ちゃんと確かめよう。
対称行列じゃ転置の具合が怪しいし。

>>> a = [[1, 0], [0, 0]]
>>> fft2(a)
array([[ 1.+0.j, 1.+0.j],
[ 1.+0.j, 1.+0.j]])
>>> fft(fft(a).transpose()).transpose()
array([[ 1.+0.j, 1.+0.j],
[ 1.+0.j, 1.+0.j]])
>>>
>>> a = [[0, 1], [0, 0]]
>>> fft2(a)
array([[ 1.+0.j, -1.+0.j],
[ 1.+0.j, -1.+0.j]])
>>> fft(fft(a).transpose()).transpose()
array([[ 1.+0.j, -1.+0.j],
[ 1.+0.j, -1.+0.j]])
>>>
>>> a = [[0, 0], [1, 0]]
>>> fft2(a)
array([[ 1.+0.j, 1.+0.j],
[-1.+0.j, -1.+0.j]])
>>> fft(fft(a).transpose()).transpose()
array([[ 1.+0.j, 1.+0.j],
[-1.+0.j, -1.+0.j]])
>>>
>>> a = [[0, 0], [0, 1]]
>>> fft2(a)
array([[ 1.+0.j, -1.+0.j],
[-1.+0.j, 1.+0.j]])
>>> fft(fft(a).transpose()).transpose()
array([[ 1.+0.j, -1.+0.j],
[-1.+0.j, 1.+0.j]])

確かにそうなっている。

逆変換ifft2(a)。

>>> a = [[1, 0], [0, 0]]
>>> ifft2(a)
array([[ 0.25+0.j, 0.25+0.j],
[ 0.25+0.j, 0.25+0.j]])

今度は(1/2)の2乗が係数に掛かっている。

逆変換ifft2(a)も、ifft(ifft(a).transpose()).transpose()と等しいはずだ。

>>> a = [[1, 0], [0, 0]]
>>> ifft2(a)
array([[ 0.25+0.j, 0.25+0.j],
[ 0.25+0.j, 0.25+0.j]])
>>> ifft(ifft(a).transpose()).transpose()
array([[ 0.25+0.j, 0.25+0.j],
[ 0.25+0.j, 0.25+0.j]])
>>>
>>> a = [[0, 1], [0, 0]]
>>> ifft2(a)
array([[ 0.25+0.j, -0.25+0.j],
[ 0.25+0.j, -0.25+0.j]])
>>> ifft(ifft(a).transpose()).transpose()
array([[ 0.25+0.j, -0.25+0.j],
[ 0.25+0.j, -0.25+0.j]])
>>>
>>> a = [[0, 0], [1, 0]]
>>> ifft2(a)
array([[ 0.25+0.j, 0.25+0.j],
[-0.25+0.j, -0.25+0.j]])
>>> ifft(ifft(a).transpose()).transpose()
array([[ 0.25+0.j, 0.25+0.j],
[-0.25+0.j, -0.25+0.j]])
>>>
>>> a = [[0, 0], [0, 1]]
>>> ifft2(a)
array([[ 0.25+0.j, -0.25+0.j],
[-0.25+0.j, 0.25+0.j]])
>>> ifft(ifft(a).transpose()).transpose()
array([[ 0.25+0.j, -0.25+0.j],
[-0.25+0.j, 0.25+0.j]])

確かに。

ifft2(a)が本当にfft(a)の逆変換か。

>>> a = [[1, 0], [0, 0]]
>>> ifft2(fft2(a))
array([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, 0.+0.j]])
>>> fft2(ifft2(a))
array([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, 0.+0.j]])
>>>
>>> a = [[0, 1], [0, 0]]
>>> ifft2(fft2(a))
array([[ 0.+0.j, 1.+0.j],
[ 0.+0.j, 0.+0.j]])
>>> fft2(ifft2(a))
array([[ 0.+0.j, 1.+0.j],
[ 0.+0.j, 0.+0.j]])
>>>
>>> a = [[0, 0], [1, 0]]
>>> ifft2(fft2(a))
array([[ 0.+0.j, 0.+0.j],
[ 1.+0.j, 0.+0.j]])
>>> fft2(ifft2(a))
array([[ 0.+0.j, 0.+0.j],
[ 1.+0.j, 0.+0.j]])
>>>
>>> a = [[0, 0], [0, 1]]
>>> ifft2(fft2(a))
array([[ 0.+0.j, 0.+0.j],
[ 0.+0.j, 1.+0.j]])
>>> fft2(ifft2(a))
array([[ 0.+0.j, 0.+0.j],
[ 0.+0.j, 1.+0.j]])

ちゃんと元に戻っている。

以上で十分な情報は得られたけれど、ついでに3×3の場合の動作を明示してしまおう。


>>> a = [[1, 0, 0], [0, 0, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.+0.j, 1.+0.j, 1.+0.j],
[ 1.+0.j, 1.+0.j, 1.+0.j],
[ 1.+0.j, 1.+0.j, 1.+0.j]])
>>> a = [[0, 1, 0], [0, 0, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j],
[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j],
[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j]])
>>> a = [[0, 0, 1], [0, 0, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j],
[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j],
[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j]])
>>> a = [[0, 0, 0], [1, 0, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.0+0.j , 1.0+0.j , 1.0+0.j ],
[-0.5-0.8660254j, -0.5-0.8660254j, -0.5-0.8660254j],
[-0.5+0.8660254j, -0.5+0.8660254j, -0.5+0.8660254j]])
>>> a = [[0, 0, 0], [0, 1, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j],
[-0.5-0.8660254j, -0.5+0.8660254j, 1.0+0.j ],
[-0.5+0.8660254j, 1.0+0.j , -0.5-0.8660254j]])
>>> a = [[0, 0, 0], [0, 0, 1], [0, 0, 0]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j],
[-0.5-0.8660254j, 1.0+0.j , -0.5+0.8660254j],
[-0.5+0.8660254j, -0.5-0.8660254j, 1.0+0.j ]])
>>> a = [[0, 0, 0], [0, 0, 0], [1, 0, 0]]
>>> fft2(a)
array([[ 1.0+0.j , 1.0+0.j , 1.0+0.j ],
[-0.5+0.8660254j, -0.5+0.8660254j, -0.5+0.8660254j],
[-0.5-0.8660254j, -0.5-0.8660254j, -0.5-0.8660254j]])
>>> a = [[0, 0, 0], [0, 0, 0], [0, 1, 0]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j],
[-0.5+0.8660254j, 1.0+0.j , -0.5-0.8660254j],
[-0.5-0.8660254j, -0.5+0.8660254j, 1.0+0.j ]])
>>> a = [[0, 0, 0], [0, 0, 0], [0, 0, 1]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j],
[-0.5+0.8660254j, -0.5-0.8660254j, 1.0+0.j ],
[-0.5-0.8660254j, 1.0+0.j , -0.5+0.8660254j]])


あと、縦横の長さは異なってOK。

>>> a = [[1, 0, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.+0.j, 1.+0.j, 1.+0.j],
[ 1.+0.j, 1.+0.j, 1.+0.j]])


やはり、fft(fft(a).transpose()).transpose()と一致する。

>>> a = [[1, 0, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.+0.j, 1.+0.j, 1.+0.j],
[ 1.+0.j, 1.+0.j, 1.+0.j]])
>>> fft(fft(a).transpose()).transpose()
array([[ 1.+0.j, 1.+0.j, 1.-0.j],
[ 1.+0.j, 1.+0.j, 1.+0.j]])
>>>
>>> a = [[0, 1, 0], [0, 0, 0]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j],
[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j]])
>>> fft(fft(a).transpose()).transpose()
array([[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j],
[ 1.0+0.j , -0.5-0.8660254j, -0.5+0.8660254j]])
>>>
>>> a = [[0, 0, 0], [0, 0, 1]]
>>> fft2(a)
array([[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j],
[-1.0+0.j , 0.5-0.8660254j, 0.5+0.8660254j]])
>>> fft(fft(a).transpose()).transpose()
array([[ 1.0+0.j , -0.5+0.8660254j, -0.5-0.8660254j],
[-1.0+0.j , 0.5-0.8660254j, 0.5+0.8660254j]])

めでたし、めでたし。