演習6

NumPy(Python+Numeric)を使ってFFTの練習

以下のプログラムを実行してみよう。 ここで、パラメータとして U を変えてみると、実空間での 広がりとスペクトル空間での広がりとの関係を調べることができる。
f1.py 
#!/bin/env python

import Numeric
import FFT
import Gnuplot

gp = Gnuplot.Gnuplot()
gp.title = ('FFT Demo')
gp('set data style lines')

L = 128    # 2のべき乗
U = 0.125  # 0.1から0.5程度で変えてみる
UL = int(L*U)

X = range(L)
f = Numeric.array([1.0]*UL+[0.0]*(L-UL))
d = Gnuplot.Data(X,f,title='Original Function')
gp.plot(d)

raw_input('Press Return')

F = FFT.fft(f)
Fr = Gnuplot.Data(X,F.real,title='FFT Real part')
Fi = Gnuplot.Data(X,F.imag,title='FFT Imaginary part')
gp.plot(Fr,Fi)

raw_input('Press Return')

2次元FFT

2次元のFFTも試してみよう。
f2.py 
#!/bin/env python

import Numeric
import FFT
import Gnuplot

gp = Gnuplot.Gnuplot()
gp.title = ('FFT Demo')
gp('set data style lines')

L = 32     # 2のべき乗
U = 0.125  # 0.1から0.5程度で変えてみる
UL = int(L*U)

X = range(L)
f = Numeric.array([1.0]*UL+[0.0]*(L-UL))
f = Numeric.outerproduct(f,f)
d = Gnuplot.GridData(f,X,X,title='Original Function')
gp.splot(d)
raw_input('Press Return')

F = FFT.fft2d(f)
Fr = Gnuplot.GridData(F.real,X,X,title='FFT Real part')
Fi = Gnuplot.GridData(F.imag,X,X,title='FFT Imaginary part')
gp.splot(Fr,Fi)

raw_input('Press Return')

実際の画像に適用

FFTの計算は要素数が2のべき乗のときに最適化されている。そこで、 画像を以下のように2のべき乗サイズにしておく。 
pnmscale -width=128 -height=128 cloud.pgm > cl.pgm
(画像は何を使っても良いが、PGM形式の白黒画像を用いる)

このプログラムでは、画像をデータとしてプロットし、 次に FFTを通したスペクトルを表示し、さらに、 波数空間でのフィルターを掛けてから逆変換して 画像にし表示する。フィルターのパラメータを変えると 画像がどのように変化するかを調べなさい。
fc.py 
#!/bin/env python

import Image
import Numeric,NumTut
import FFT
import Gnuplot

L = 128                   # 画像のサイズ(L x L)
U = 0.25                  # フィルターのパラメータ
im = Image.open('cl.pgm') # 画像ファイル

f  = Numeric.array(im.getdata())
f.shape = (L,L)
X = range(L)

gp = Gnuplot.Gnuplot()
gp('set data style lines')

d = Gnuplot.GridData(f,X,X,title='Original Data')
gp.splot(d)
raw_input('Press Return')
F = FFT.fft2d(f)
Fr = Gnuplot.GridData(F.real,X,X,title='FFT Real part')
Fi = Gnuplot.GridData(F.imag,X,X,title='FFT Imaginary part')
gp.splot(Fr,Fi)
raw_input('Press Return')

UL = int(L*U)
G = Numeric.array([1.0]*(UL/2)+[0.0]*(L-UL)+[1.0]*(UL/2))
G = Numeric.outerproduct(G,G)
F = F * G
f = FFT.inverse_fft2d(F)
NumTut.view(Numeric.transpose(f.real))
最後に、xvのsave機能を使って表示された画像を保存してレポート作成に用いること。

畳み込みによるフィルター

FFTの試験に用いた画像(この例ではcl.pgm)に5x5の低域通過フィルター (先週の授業を参照)をかけてその効果を調べレポートにしなさい。
本日の演習のレポートは ~/public_html/image/ex6.html に作成すること。

参考：
Numerical Python Manual
Gnuplot.py
Python Imaging Library