読者です 読者をやめる 読者になる 読者になる

遊び tokidoki 仕事

数学と音楽と教育と遊び

Top | おしごと | ゼミ | がくせい | すうがく | かがく | きょういく | おんがく | おきにー | Tips | Photo | イベント | ものもう | あれこれ | About

フーリエでモアレ

フーリエ級数と周期性ノイズについてまとめている学生のネタ。
2次元のうなりとしてのモアレを
わざわざフーリエ変換したデータを傷つけることによって
モアレを擬似的に起こすシミュレーションをした。
以下がそのプログラム。
やっつけで作っているので、相当汚い。


上段が元画像とそのフーリエ変換、
下段がフーリエ変換データに傷を付けたものと
その逆フーリエ変換。

モアレノイズがどう消せるのかの
一つの考察。


こんな感じ↓





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

REM [FFT2D_moare.bas]
REM 2D FFT/IFFT でモアレを合成する
REM FFT/IFFTは外部関数として使う

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

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

! 配列サイズ dn は必ず2^nで
LET dn = 256            ! データ総数
LET  dn2 = dn/2
LET  win=.02            ! windowの余白

LET  dt = 1   ! サンプリング間隔
LET  ct= 10^0  !元画像のコントラスト
LET  mct= 10^(.5)  !モアレ画像のコントラスト
LET  fct= 10^(-.5)  !フーリエスペクトルのコントラスト

!極座標で(r,θ)=±(mr,mth)の周波数のところにmstの強さで傷を付ける
LET  mst=100  !モアレ傷の強さ
LET  mth=30   !モアレの角度
LET  mr=.5 !モアレの細かさ

DEF f(x,y)=(1-SGN(x^2+y^2-.5)-.1)+10*(1-SGN(2*(x-.2)^2+(y-.2)^2-.01)-.1)
 +10*(1-SGN(2*(x+.2)^2+(y-.2)^2-.01)-.1)
 +10*(MAX(-SGN(x^2-y-.3)+SGN(x^2/2-y-.2)-.1,0))

SET POINT COLOR 10
SET POINT STYLE 1
SET bitmap SIZE dn*2*(1+win*2),dn*2*(1+win*2)

REM 【元の画像の表示】

SET VIEWPORT win,.5-win , .5+win,1-win
DIM bmp(dn,dn)          ! 画像用配列 dn×dn

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

CALL FFT(bmp)

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

SET VIEWPORT .5+win,1-win , .5+win,1-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)

REM 【モアレ傷を付ける】

SET VIEWPORT .5+win,1-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2
LET  mx=mr*COS(mth/180*PI)
LET  my=mr*SIN(mth/180*PI)

FOR y= -dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      IF (ABS(x-dn2*mx)<1 AND ABS(y-dn2*my)<1) THEN
         LET  bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1)=mst
      END IF
      IF (ABS(x+dn2*mx)<1 AND ABS(y+dn2*my)<1) THEN
         LET  bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1)=mst
      END IF
      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)

CALL IFFT(bmp)

REM 【モアレ画像の表示】
SET VIEWPORT win,.5-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2

FOR y= 1 TO dn
   FOR x = 1 TO dn
      LET g=ATN(re(bmp(x,y)+.1)*mct)/PI/2+.5
      SET COLOR mix(10) g,g,g
      PLOT POINTS: x-dn2,y-dn2
   NEXT x
NEXT y

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 NATIVE
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 NATIVE
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 NATIVE
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