● 波動光学

Youngの干渉実験でお馴染みのように光の波動性は干渉として物理的に明確と なる。干渉現象も含めて光の伝搬の様子を光の波動性から数学的に導く手法が 波動光学である。後で解説するように光の干渉とFourier変換には密接な関連 がある。本節では光の干渉を計算機で再現することを通じて、Fourier変換の 意味を理解することを目指そう。

■ 波動光学の理論

一本の光路を走る波長λの光の波動関数 f を光路上の距離 l と 時刻 t によって次のように表す。

(1)

(2)

(3)

(4)

(5)

フィルムが非常に小さいならば、(5)式中の x^2+y^2 による位相の変化も無視でき簡単な次式となる。

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(12)式は次の通りとなる。

(13)

(13)式は各 k について N個の格子点での計算が 必要なので全部のkについて合計 N^2回の計算が必要に思われるが、実は 高速Fourier変換(Fast Fourier Transform(FFT))と呼ばれるアルゴリズム があり、この方法では N log_2{N} のオーダーの計算量で済む。 単純な比較は危険だが、N=256でのN^2とN log_2{N}の差は歴然である。

このFFTのアルゴリズムを簡単に紹介しよう。 FFTでは系のサイズ、つまり格子点の個数 N が例えば256=2^8のように 2の巾乗となる場合にもっとも威力を発揮する。なのでNをそのような 値にとることにしよう。

(13)式のjの和について、jが偶数の項と奇数の 項で別々に和をとるように変形する。

(14)
(15)

この論法をさらに繰り返すと最終的にはサイズ1のN本のFourier変換に 位相をかけた和になる。サイズ1のFourier変換とは恒等変換であるから、 ここから逆に分解を遡って位相をかけて足し合わせる作業を進めていくと Fourier変換になるのである。 N=8の簡単な場合についてそのF_0とF_1のこの分解による表現を記しておく。

(16)
(17)
(18)
(19)

class RevTable
{
  int reve[N];
public:
  RevTable( void ){
    for( int i=0 ; i<N ; i++ ){
      int rev=0;
      for( u_int mask1=1, mask2=N/2 ; mask2 ; mask1<<=1, mask2>>=1 ){
        if( i&mask1 ) rev|=mask2;
      }
      reve[i] = rev;
    }
  }
  inline int operator [] ( int i ){
    return( reve[i] );
  }
} Rev;

class ExpTable
{
  Complex expo[2*N];
public:
  ExpTable( void ){
    for( int p=-N+1 ; p<=N-1 ; p++ ){
      expo[p+N] = Complex( cos(2*M_PI/N*p), sin(2*M_PI/N*p) );
    }
  }
  inline Complex operator [] ( int p ){
    return( expo[p+N] );
  }
} Exp;

//---- 1-Dimensional Fast Fourier Transform
void FFT1( int sign, Complex src[N], Complex dst[N] )
{
  int k;
    // copy from source array at reversed index to destination array
  for( int j=0 ; j<N ; j++ ) dst[Rev[j]] = src[j];

    // Uncomment below lines, and x=0, k=0 are treated as the center of the arrays.
  /**********************
  for( k=N/2 ; k<N ; k++ ) dst[k] *= -1.0;
  **********************/

    // Danielson-Lanczos's algorithm
  int k0, k1, k2, dk, p, dp;
  Complex ex, work;
  for( dk=0x01, dp=HALF_N ; dp ; dk<<=1, dp>>=1 ){
    for( k0=0, p=0 ; k0<dk ; k0++, p+=sign*dp ){
      for( k1=k0, ex=Exp[p] ; k1<N ; k1=k2+dk ){
        k2=k1+dk;
        work     = dst[k2] * ex;
        dst[k2]  = dst[k1] - work;
        dst[k1] += work;
      }
    }
  }
    // Uncomment below lines, and x=0, k=0 are treated as the center of the arrays.
  /**********************
  for( k=1 ; k<N ; k+=2 ) dst[k] *= -1.0;
  **********************/
}

(20)

波数の配列への格納の仕方に注意がいる。波数の配列のindex と (13)式の波数kの対応関係を以下に記す。

index	0	1	…	N/2	…	N-2	N-1
k	0	+1	…	+N/2,,-N/2	…	-2	-1

この対応関係は奇妙に思えるかも知れないが、これが世に出回っている Fourier変換関数の標準仕様である。 しかしこの仕様では光の干渉を扱うには不自然である。 フィルムとスクリーンの中心が共に配列の中心に配置させるには FFT1() 関数の2つのコメントを解除する。 こうすると以下の対応関係となる。

index	0	1	…	N/2	…	N-2	N-1
k	-N/2	1-N/2	…	0	…	N/2-2	N/2-1

▲ 2次元の高速Fourier変換

サイズ N times Nの2次元のFFTは x軸方向のN本の1次元FFTを行ってから、 y軸方向のN本の1次元FFTを行うので、複素数の演算回数は合計

(21)

//---- 2-Dimensional Fast Fourier Transform
void FFT2( int sign, Complex src[N][N], Complex dst[N][N] )
{
  int kx, ky;
    // copy from source array at reversed index to destination array
  for( int jx=0 ; jx<N ; jx++ ){
    for( int jy=0 ; jy<N ; jy++ ){
      dst[Rev[jx]][Rev[jy]] = (1.0/N)*src[jx][jy];
    }
  }
    // Uncomment below lines, and x=0, k=0 are treated as the center of the arrays.
  /**********************
  for( kx=N/2 ; kx<N ; kx++ ){
    for( ky=0 ; ky<N/2 ; ky++ ){
      dst[kx][ky] *= -1.0;      dst[ky][kx] *= -1.0;
    }
  }
  **********************/

    // Danielson-Lanczos's algorithm
  int k0, k1, k2, dk, p, dp;
  Complex ex, work;
  for( dk=0x01, dp=HALF_N ; dp ; dk<<=1, dp>>=1 ){
    for( k0=0, p=0 ; k0<dk ; k0++, p+=sign*dp ){
      for( k1=k0, ex=Exp[p] ; k1<N ; k1=k2+dk ){
        for( ky=0, k2=k1+dk ; ky<N ; ky++ ){      // transform columns
          work         = dst[k2][ky] * ex;
          dst[k2][ky]  = dst[k1][ky] - work;
          dst[k1][ky] += work;
        }
        for( kx=0, k2=k1+dk ; kx<N ; kx++ ){      // transform rows
          work         = dst[kx][k2] * ex;
          dst[kx][k2]  = dst[kx][k1] - work;
          dst[kx][k1] += work;
        }
      }
    }
  }
    // Uncomment below lines, and x=0, k=0 are treated as the center of the arrays.
  /**********************
  for( kx=0 ; kx<N ; kx++ ){
    for( ky=kx+1 ; ky<N ; ky+=2 ){
      dst[kx][ky] *= -1.0;      dst[ky][kx] *= -1.0;
    }
  }
  **********************/
}

▲ 畳み込み演算

畳み込み演算(convolution)と呼ばれるFourier変換特有の演算がある。 これは大事な法則を含むので学んでおこう。畳み込みの演算子 * の 演算規則は次式で定義される。

(22)

(23)
(24)
(25)

(26)

(27)

小さな穴による光の回折をFraunhofer回折の式 5()式 で計算する。プログラムソースは fraunhofer.cc である。

下図のように四角形、円形の開口があるフィルムに よる干渉模様を計算すると、それぞれ下図になる。

開口の部分で f(Q)=1 で、他でf(Q)=0として、Fourier変換を 代数的に計算すると、四角形の開口は

(28)

(29)

もっとも単純な原理による簡単な像のholographyを載せる。 propagatorには単純なFourier変換 (12)式を用いる。 ソースファイル名は holography.cc である。

元の像

holographyに撮影された像

再現された像

▲ 他の応用例

Fourier変換のデータ処理における重要性は色々あるが、ひとつに、繰り返し 構造をもつ像の特徴が極めて明確になり精密な分析が可能になることである。 ふたつに、繰り返し構造をもたないノイズを識別することが容易になり、 データからノイズを除去できることである。 このため、電子顕微鏡によるタンパク質分析や、X線,電子線による結晶表面解 析などで実際に活用されている。近年では電子線を用いたholographyによる 結晶解析の研究も盛んで、結晶試料の干渉像から結晶の構造が直ちに判明でき るようになりつつある。

また、偏微分方程式の解析にもFourier変換は使われる。 実空間で方程式を解くより、波数空間で解く方が安定性や精度の面で 優れる場合がある。高Reynolds numberでの流体 footnote{系の対称性が必要で、工学の設計ではあまり使われない。} や物質内部の電子状態 footnote{最近は実空間でより良く解く研究が盛んである。} の計算などがある。

本節で示したように、Fourier変換の本質は目に見えて楽しい光の干渉現象 から理解することができるのである。例に載せた像以外にも各自さまざまな 像の干渉を計算して、Fourier変換の性質の理解を深めてほしい。