遊び tokidoki 仕事

数学と音楽と教育と遊び

| おしごと - きょういく - がくせい - ゼミ - イベント | すうがく - おんがく - 数理音楽 - DTM - かがく - scratch
| Art - photo - おきにー - Tips - ものもう - あれこれ | About - Top

超解像シミュレーション

すうがく ゼミ Tips

痛すぎる卒論指導も終盤。
ひょんなことから一人のゼミ生には
超解像の原理についてまとめてもらっている。
原理は関数解析的な話になるがそれはさておき、
実際にそれができるものなのか試したくなった。
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