痛すぎる卒論指導も終盤。
ひょんなことから一人のゼミ生には
超解像の原理についてまとめてもらっている。
原理は関数解析的な話になるがそれはさておき、
実際にそれができるものなのか試したくなった。
10進BASICで簡単なシミュレーションをやってみたので
そのプログラムを公開。
ただ、かなり汚い書き方なのでご容赦を。
【使い方】
とにかく実行すると超解像を順次繰り返します。
画像をファイルとして保存したいときは
!GSAVE
となっている行の ! を外してください。
こんなかんじ↓
10進BASICプログラム→[superResolution1D.bas]
REM [superResolution1D.bas] REM 超解像のシミュレーション REM 帯域制限された1D波形を回復するシミュレーション REM REM FFTおよびIFFTサブルーチンを使用するために必要な宣言 OPTION ARITHMETIC complex ! 複素数モードに設定 DECLARE EXTERNAL SUB FFT ! FFT関数 DECLARE EXTERNAL SUB IFFT ! FFT関数 REM FFTのための配列および定数宣言 LET dn = 1024*16 ! データ総数 DIM d(dn) ! 入出力バッファ用配列 DIM g(dn) ! 帯域制限領域周波数用配列 REM 【元の関数の定義】 DEF f(x)=MIN(MAX(2-(x/.0032)^2,0),1.8) LET dt = .001 ! サンプリング間隔 LET XB=.4 !空間幅(割合) LET FB=.4 !周波数制限(割合) LET ht=2.2 !元の表示windowの高さ LET wd=dt*dn/64 !元の表示wondowの幅 LET fht=250 !周波数領域の表示windowの高さ LET win=.01 LET boxx=2 LET boxy=3 LET bmp=300 SET bitmap SIZE bmp*(1+win*2)*boxx,bmp*(1+win*2)*boxy LET boxz=1/MAX(boxx,boxy) SET AREA COLOR 0 REM 【元の関数の表示】==================================================== SET VIEWPORT win,boxz-win , 1-boxz+win,1-win SET WINDOW -wd*1.1, wd*1.1, -ht, ht DRAW axes0 LET t = 0 FOR i = 1 TO dn LET d(i) = f(i/dn-.5) PLOT LINES: (i-dn/2)*dt, d(i); LET t = t + dt NEXT i PLOT TEXT ,AT -wd,ht*.9: "元の信号" REM 【周波数リストの作成】(後半半分のデータは無視する) dim freq(dn/2+1) let maxFreq = 1/dt / 2 ! 最大周波数 for i=1 to dn/2+1 let freq(i) = maxFreq * (i-1) / (dn/2) NEXT i REM 【帯域制限定数】 LET tauX=wd*XB/dt !空間制限 LET tau=INT(maxFreq*FB) !周波数制限 REM 【元信号の周波数スペクトラム】================================================ SET VIEWPORT boxz+win,boxx*boxz-win , 1-boxz+win,1-win SET WINDOW -maxFreq/100, maxFreq/10, -fht, fht DRAW axes(maxFreq/20, fht/2) LET t = 0 call FFT(d) PLOT LINES FOR i = 1 TO dn/2+1 PLOT LINES: freq(i), re(d(i)); IF i>tau THEN LET d(i)=0 LET g(i)=d(i) NEXT i PLOT TEXT ,AT maxFreq*.02,fht*.9: "元の信号のフーリエ変換" LET c=0 !超解像の適用回数カウンタ DO WHILE(c<1000) REM 【周波数帯域制限】==================================================== SET VIEWPORT boxz+win,boxx*boxz-win , 1-boxz*2+win,1-boxz-win SET WINDOW -maxFreq/100, maxFreq/10, -fht, fht PLOT AREA : -maxFreq,-fht;-maxFreq,fht; maxFreq,fht; maxFreq,-fht DRAW axes(maxFreq/20, fht/2) PLOT LINES FOR i = 1 TO dn/2+1 IF i<=tau THEN LET d(i)=g(i) PLOT LINES: freq(i), re(d(i)); NEXT i PLOT LINES SET LINE COLOR 2 PLOT LINES: freq(tau),-fht*.9; freq(tau),fht*.9 PLOT TEXT ,AT freq(tau), -fht*.4: "τ" PLOT TEXT ,AT maxFreq*.02,fht*.9: "帯域制限信号との合成" REM 【フーリエ逆変換して表示】============================================ SET VIEWPORT win,boxz-win , 1-boxz*2+win,1-boxz-win SET WINDOW -wd*1.1, wd*1.1, -ht, ht PLOT AREA : -wd*1.1,-ht;-wd*1.1,ht; wd*1.1,ht; wd*1.1,-ht DRAW axes0 LET t = 0 PLOT lines CALL IFFT(d) SET LINE COLOR 1 FOR i = 1 TO dn PLOT LINES: (i-dn/2)*dt, d(i); LET t = t + dt NEXT i PLOT TEXT ,AT -wd,ht*.9: "超解像"&STR$(c)&"回目" REM 【空間帯域制限】====================================================== SET VIEWPORT win,boxz-win , win,1-boxz*2-win SET WINDOW -wd*1.1, wd*1.1, -ht, ht PLOT AREA : -wd*1.1,-ht;-wd*1.1,ht; wd*1.1,ht; wd*1.1,-ht DRAW axes0 LET t = 0 PLOT LINES FOR i = 1 TO dn IF ABS(i-dn/2)>tauX THEN LET d(i)=0 PLOT LINES: (i-dn/2)*dt, d(i); LET t = t + dt NEXT i PLOT LINES SET LINE COLOR 2 PLOT LINES: tauX*dt,-ht*.9; tauX*dt,ht*.9 PLOT LINES: -tauX*dt,-ht*.9; -tauX*dt,ht*.9 PLOT TEXT ,AT -wd,ht*.9: "空間制限" PLOT TEXT ,AT -tauX*dt, -ht*.2: "-T" PLOT TEXT ,AT tauX*dt*1.1, -ht*.2: "T" REM 【超解像周波数スペクトラム】================================================ SET VIEWPORT boxz+win,boxx*boxz-win , win,1-boxz*2-win SET WINDOW -maxFreq/100, maxFreq/10, -fht, fht PLOT AREA : -maxFreq,-fht;-maxFreq,fht; maxFreq,fht; maxFreq,-fht DRAW axes(maxFreq/20, fht/2) LET t = 0 call FFT(d) PLOT LINES SET LINE COLOR 1 FOR i = 1 TO dn/2+1 PLOT LINES: freq(i), re(d(i)); NEXT i PLOT TEXT ,AT maxFreq*.02,fht*.9: "空間制限信号のフーリエ変換" REM 【画像ファイルとして出力】========================================== LET sw=0 IF MOD(c,50)=0 THEN LET sw=1 ELSEIF (c<31 AND c>9) AND MOD(c,10)=0 THEN LET sw=1 ELSEIF c<10 THEN LET sw=1 END IF IF sw=1 THEN LET f$="spr"&right$("000"&STR$(c),3)&".bmp" !gsave f$, "4bit" END IF LET c=c+1 LOOP END REM ====================================================================== REM [1D 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 【1D FFT】 REM θ = -2π/N REM W = exp(-��(-1)θ) EXTERNAL SUB FFT(x()) OPTION ARITHMETIC complex DECLARE EXTERNAL SUB FFTMAIN LET sx = SIZE(x) LET th = -2 * PI / sx CALL FFTMAIN(x, th) END SUB REM 【1D IFFT】 REM θ = 2π/N REM W = exp(-��(-1)θ) EXTERNAL SUB IFFT(x()) OPTION ARITHMETIC complex DECLARE EXTERNAL SUB FFTMAIN LET sx = SIZE(x) LET th = 2 * PI / sx CALL FFTMAIN(x, th) MAT x = (1/sx) * x 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