フーリエ級数と周期性ノイズについてまとめている学生のネタ。
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
