２次元FFTで画像のフーリエ変換

引き続き、卒論用のプログラムシリーズ。
今回は、画像のフーリエ変換を表示する10進BASICプログラム。
たとえば円環領域のフーリエ変換は

といった感じに出力。


10進BASICプログラム→[FFT2D.bas]

REM [FFT2D.bas]
REM 2D FFT/IFFT
REM FFT/IFFTは外部関数として使う

OPTION ARITHMETIC complex
SET TEXT FONT "",6
SET AXIS COLOR 2

DECLARE EXTERNAL SUB FFT    ! FFT
DECLARE EXTERNAL SUB IFFT   ! IFFT

! 配列サイズ dn は必ず2^nで
LET dn = 256            ! データ総数
LET  dn2 = dn/2
LET  win=.02            ! windowの余白
DIM bmp(dn,dn)          ! 画像用配列 dn×dn
SET bitmap SIZE dn*2*(1+win*2),dn*(1+win*2)

REM 【元の関数の定義】
DEF f(x,y)=SGN(x^2+y^2-.45)-SGN(x^2+y^2-.5)-.1

LET  dt = 1   ! サンプリング間隔
LET  ct= 10^2  !元画像のコントラスト
LET  fct= 10^1  !フーリエスペクトルのコントラスト

SET POINT COLOR 10
SET POINT STYLE 1

REM 【元の画像の表示】

SET VIEWPORT win,.5-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2
FOR y=-dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      LET  g=f(x/dn2,y/dn2)
      LET bmp(x+dn2+1,y+dn2+1) = g
      LET  g=ATN(g*ct)/PI+.5
      SET COLOR MIX(10) g,g,g 
      PLOT POINTS: x,y
   NEXT x
NEXT y
DRAW axes(dn/8,dn/8)

CALL FFT(bmp)

REM 【周波数スペクトラム】

SET VIEWPORT .5+win,1-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2
FOR y= -dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      LET g=ATN(ABS(bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1))*fct)/PI*2
      SET COLOR MIX(10) g,g,g 
      PLOT POINTS: x,y
   NEXT x
NEXT y
DRAW axes(dn/8,dn/8)

END

REM  [2D FFT/IFFT]
REM   元データ x(,) → FFT/IFFT → 変換後データ x(,)
REM   ユニタリ変換になるようFFT/IFFT両方を1/N
REM   Fk = 1/N * Σ[n=1...N] (fn * W^( (k-1)*(n-1)) )   :FFT
REM   fk = 1/N * Σ[n=1...N] (Fn * W^(-(k-1)*(n-1)) )  :IFFT
REM   W = exp(-��(-1) * 2π/N)
REM   配列添字は 1 から出発

REM 【2D FFT】
REM θ = -2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB FFT(x(,))
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x,1)
DIM d(sx)
LET th = -2 * PI / sx
FOR t=1 TO sx
   FOR s=1 TO sx
      LET  d(s)=x(s,t)
   NEXT s
   CALL FFTMAIN(d, th)
   FOR s=1 TO sx
      LET  x(s,t)=d(s)
   NEXT s
NEXT t
MAT x = (1/sx) * x
FOR s=1 TO sx
   FOR t=1 TO sx
      LET  d(t)=x(s,t)
   NEXT t
   CALL FFTMAIN(d, th)
   FOR t=1 TO sx
      LET  x(s,t)=d(t)
   NEXT t
NEXT s
END SUB

REM 【2D IFFT】
REM θ = 2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB IFFT(x(,))
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x,1)
DIM d(sx)
LET th = 2 * PI / sx
FOR t=1 TO sx
   FOR s=1 TO sx
      LET  d(s)=x(s,t)
   NEXT s
   CALL FFTMAIN(d, th)
   FOR s=1 TO sx
      LET  x(s,t)=d(s)
   NEXT s
NEXT t
MAT x = (1/sx) * x
FOR s=1 TO sx
   FOR t=1 TO sx
      LET  d(t)=x(s,t)
   NEXT t
   CALL FFTMAIN(d, th)
   FOR t=1 TO sx
      LET  x(s,t)=d(t)
   NEXT t
NEXT s
END SUB

REM [バタフライ演算]
REM 配列サイズを半分にしながら再帰呼び出しする
REM wk= exp( (k-1)θ )
REM x(k)    --- ×1  -- x0(k)
REM         × 
REM x(k+N/2)--- ×wk -- x1(k)

EXTERNAL SUB FFTMAIN(x(), th)
OPTION ARITHMETIC complex
LET sx = SIZE(x)
IF sx>1 THEN
   LET sx2 = sx/2
   DIM x0(sx2), x1(sx2)
   FOR k=1 TO sx2
      LET x0(k) = x(k) + x(k+sx2)
      LET wk = EXP( complex(0, th * (k-1)) )
      LET x1(k) = wk * ( x(k) - x(k+sx2) )
   NEXT k
   LET th0 = 2 * th
   CALL FFTMAIN(x0, th0)
   LET th1 = 2 * th
   CALL FFTMAIN(x1, th1)
   FOR k=1 TO sx2
      LET x(2*k-1) = x0(k)
      LET x(2*k) = x1(k)
   NEXT k
END IF
END SUB