2次元FFTと点像関数

1次元のFFTは時系列データの解析によく使われますが、2次元のFFTは画像処理など空間周波数に関する解析に多く用いられます。

Maple には2次元の FFT は用意されていないのでここでは Maple での2次元 FFT のプログラムを紹介します。

ただしプログラムの解説は省略しその応用例を中心にご紹介します。光学では点像関数の計算に使用されよく知られています。カメラや顕微鏡などは像を結ぶ結像光学系です。結像光学系は点から出た光を像面において点に結びます。点の像は光の回折により点にはなりません。点像関数はこの点の像を表します。点光源からの波面は光学系を出たところで球面波となり1点に結ぼうとします。波面は像点を中心とした球(参照球)に張り付いたようになっています。

参照球上の波面の振幅と位相をフーリエ変換(2次元)したものが点像関数です。射出瞳上(参照球に相当)の振幅と位相を表す関数を瞳関数と呼びます。

>    restart;

2次元 FFT の Maple プログラム

最初にワークシートの初期の準備を行います。

>    with(linalg):with(plots):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

以下は2次元高速フーリエ変換プログラムの本体です。パラメータ isign の値によって

逆フーリエ変換も計算されます。

>    FFT2_body:=proc(n,A,isign)local B,data,ndim,i,j,k,ntot,nprev,idim,nrem,ip1,ip2,ip3,i2rev,
i1,i2,i3,i3rev,tempr,tempi,theta,ibit,ifp1,ifp2,wpr,wpi,wr,wi,k1,k2,wtemp;
ndim:=2;
data:=hfarray(1..n*n*2):
for i to n do
   for j to n do
      k:=j*n+i-n;
      data[k*2-1]:=evalf(evalc(Re(A[i,j])));
      data[k*2]:=evalf(evalc(Im(A[i,j])));
   od;
od:
ntot:=n*n:
nprev:=1;
for idim to ndim do
   nrem:=ntot/(n*nprev);
   ip1:=2*nprev;
   ip2:=ip1*n;   
   ip3:=ip2*nrem;
   i2rev:=1;
   for i2 by ip1 to ip2 do
      if i2<i2rev then
         for i1 from i2 by 2 to i2+ip1-2 do
            for i3 from i1 by ip2 to ip3 do
               i3rev:=i2rev+i3-i2;
               tempr:=data[i3];
               tempi:=data[i3+1];
               data[i3]:=data[i3rev];
               data[i3+1]:=data[i3rev+1];
               data[i3rev]:=tempr;
               data[i3rev+1]:=tempi;
            od;
         od;
      fi:
      ibit:=ip2/2;
      while (ibit>=ip1) and (i2rev>ibit) do
         i2rev:=i2rev-ibit;
         ibit:=ibit/2;
      od;
      i2rev:=i2rev+ibit;
   od;
   ifp1:=ip1;
   while   ifp1<ip2 do
      ifp2:=2*ifp1;
      theta:=isign*6.28318530717959/(ifp2/ip1);
      wpr:=-2.0*sin(0.5*theta)**2;
      wpi:=sin(theta);
      wr:=1.0;
      wi:=0.0;
      for i3 by ip1 to ifp1 do
         for i1 from i3 by 2 to i3+ip1-2 do
            for i2 from i1 by ifp2 to ip3 do
               k1:=i2;     k2:=k1+ifp1;     
               tempr:=evalhf(wr*data[k2]-wi*data[k2+1]);               
               tempi:=evalhf(wr*data[k2+1]+wi*data[k2]);
               data[k2]:=evalhf(data[k1]-tempr);
               data[k2+1]:=evalhf(data[k1+1]-tempi);                     
               data[k1]:=evalhf(data[k1]+tempr);
               data[k1+1]:=evalhf(data[k1+1]+tempi);
         od;od:
         wtemp:=wr;  
         wr:=evalhf(wr*wpr-wi*wpi+wr);
         wi:=evalhf(wi*wpr+wtemp*wpi+wi);
      od;   
      ifp1:=ifp2;
   od;
   nprev:=n*nprev;
od:
B:=matrix(n,n):
for i to n do
   for j to n do
   k:=n*j+i-n;
   B[i,j]:=evalf(data[k*2-1])+I*evalf(data[k*2]):
   od;
od;
RETURN(B):
end:

isign を 1 で呼ぶと通常の2次元の高速フーリエ変換プログラムになります。

>    FFT2:=proc(n,A)FFT2_body(n,A,1)end:

isign を -1 で呼ぶと逆フーリエ変換プログラムになります。

>    invFFT2:=proc(n,A)FFT2_body(n,A,-1)end:

FFTでは点像の中心を行列の真中にシフトする必要があります。このためのシフトを行う関数を用意しておきます。

>    FFTshift:=proc(n,A)local i,j,B:
B:=matrix(n,n):
for i to n do
for j to n do
B[i,j]:=A[(i+n/2-1 mod n) +1,(j+n/2-1 mod n) +1]:
od:od:
RETURN(B)
end:

簡単な例

簡単なフーリエ変換を試してみましょう

>    n:=16:

>    P:=matrix(n,n):

>    for i to n do
x:=2*(i-n/2-1/2)/n;
for j to n do
y:=2*(n/2-j+1/2)/n;
P[i,j]:=0;
if  abs(x)< .6  and abs(y)<0.6  then P[i,j]:=1 fi:
od:od:
p:=NULL:
for th from 0  to 35 do
p:=p,listplot3d(P,orientation=[th/4,th]);
od:
display(p,insequence=true);

[Maple Plot]

この図は振幅が一定で、四角い窓(瞳)のレンズでの点像の例に相当します。位相はゼロで波面は乱れがなくきれいにそろっていると仮定します。これは収差のない理想レンズによる点像関数に相当します。これに2次元 FFT をかけます。

>    Q:=FFT2(n,P):

結果は複素行列なので強度(intensity)を求める関数を用意しておきます。

>    get_intensity:=proc(n,A)local B,i,j:
B:=matrix(n,n):
for i to n do
for j to n do
B[i,j]:=evalf(Re(A[i,j])^2+Im(A[i,j])^2):
od:od:
RETURN(B):
end:

>    R:=get_intensity(n,Q):

シフトして中心を合わせプロットします。

>    S:=FFTshift(n,R):
listplot3d(S,shading=zhue);

[Maple Plot]

これが四角い窓(矩形の瞳)の場合の点像関数です。

視点を動かしてみましょう。

>    p:=NULL:
for th from 0  to 35 do
p:=p,listplot3d(S,shading=ZHUE,orientation=[th/5,th]);
od:
display(p,insequence=true);

[Maple Plot]

レンズの窓(口径)が小さくなったら

レンズのサイズ(システム開口)が小さくなるとどうなるでしょうか?

同じ矩形(正方形)のアパチャーの一様な振幅(位相はゼロ)のデータを用意します。

>    for i to n do
x:=2*(i-n/2-1/2)/n;
for j to n do
y:=2*(n/2-j+1/2)/n;
P[i,j]:=0;
if  abs(x)< .25  and abs(y)<0.25  then P[i,j]:=1 fi:
od:od:

>    listplot3d(P);

[Maple Plot]

これをFFTにかけてみましょう。

>    Q:=FFT2(n,P):

>    R:=get_intensity(n,Q):

>    S:=FFTshift(n,R):
listplot3d(S,shading=zhue);

[Maple Plot]

このように点像関数は前に比べ回折により広がったものになります。この広がりは光の波長によってことなります。

>    listplot3d(S,shading=zhue,orientation=[0,0],style=patchcontour,contours=5);

[Maple Plot]

この点像関数の物理的サイズが重要になります。結果の表示されたデータ行列の大きさは波長(λ)およびレンズのF数(F number)によってきまります。データ行列のグリッド格子間隔サイズGは次のようにきまります。

G=λ*F数

波長が500ナノメートルでF数が2とするとグリッド間隔は1000ナノメートル(=1ミクロン)になります。

点像関数の断面を表示してみましょう。

>    listplot3d(S,shading=zhue,orientation=[0,90]);

[Maple Plot]

このように点像関数は一旦ゼロになりその後も波を打っています。最初にゼロになるところがありますがこれを第一暗輪(この場合は回転対称ではないが)と呼びます。この第一暗輪の円盤(回転対称の場合)はエアリーディスクと呼ばれています。この場合エアリーディスクの半径は4グリッド分となっています。従ってさきほどのGを使うと4xGの半径がエアリーディスクのサイズとなります。

エアリーディスク

このように一様な振幅の場合、円形の瞳で回転対象の場合は点像関数はベッセル関数 (J1) で表されます。

>    J1 := x -> BesselJ(1,x);

J1 := proc (x) options operator, arrow; BesselJ(1,x) end proc

強度分布は次のようになります。

>    Intensity:=r->(2*J1(b*r)/b/r)^2;

Intensity := proc (r) options operator, arrow; 4*J1(b*r)^2/b^2/r^2 end proc

bはF数および波長によってきまる定数です。F数を1波長を0.5ミクロンとしてプロットしてみましょう。

>    F:=1;lambda:=0.5;

F := 1

lambda := .5

>    b:=Pi/lambda/F;

b := 2.000000000*Pi

>    plot(Intensity(r),r=0..2);

[Maple Plot]

b*rを変数にしてプロットしてみましょう。

>    plot(J1(br),br=0..5);

[Maple Plot]

これがゼロになる点を計算すると次のようになります。

>    fsolve(J1(br)=0,br=1..4);

3.831705970

フーリエ変換にかける行列全体の振幅を1にしてみましょう。

>    for i to n do
x:=2*(i-n/2-1/2)/n;
for j to n do
y:=2*(n/2-j+1/2)/n;
P[i,j]:=1;
od:od:
listplot3d(P,view=0..1,axes=box);

[Maple Plot]

これをフーリエ変換して点像関数をプロットしてみます。

>    Q:=FFT2(n,P):

>    R:=get_intensity(n,Q):

>    S:=FFTshift(n,R):
listplot3d(S,shading=zhue,orientation=[0,90]);

[Maple Plot]

このようにグリッドひとつ分がエアリーディスクとなってしまいます。このため瞳だけでフーリエ変換したのでは点像関数の詳細が(特に収差が小さいレンズの場合)わかりません。このため通常瞳の2倍くらいの範囲をフーリエ変換するのが通常行われています。

円形開口の場合

円形の瞳(開口)の人も関数を作成します。

>    n:=32;P:=matrix(n,n):

n := 32

>    for i to n do
x:=2*(i-n/2-1/2)/n;
for j to n do
y:=2*(n/2-j+1/2)/n;
P[i,j]:=0;
if evalf(sqrt(x^2+y^2)) < 0.7^2  then P[i,j]:=1 fi:
od:od:

>    listplot3d(P);

[Maple Plot]

点像関数を求めます。

>    Q:=FFT2(n,P):

>    R:=get_intensity(n,Q):

>    S:=FFTshift(n,R):

>    listplot3d(S,axes=boxed);

[Maple Plot]

対数スケールでプロットします。

>    listplot3d([seq([seq(log(S[i,j]+0.0001),i=2..n)],j=2..n)],shading=zhue,view=-5..15,orientation=[1,1]);

[Maple Plot]

中心部分だけを拡大してプロットしましょう。

>    listplot3d([seq([seq(log(S[i,j]),i=n/2-5..n/2+7)],j=n/2-5..n/2+7)],shading=zhue);

[Maple Plot]

ガウシアン・ビーム

CDやDVDに使われているピックアップ・レンズでは光源にレーザー・ダイオードが使われています。これらのビームは瞳上で強度が一様でなく正規(ガウス)分布とみなされます。また瞳上で回転対象でなく縦方向と横方向で分布の広がり幅がことなります。このような場合の点像関数を求めてみましょう。

>    n:=32;P:=matrix(n,n):

n := 32

>    e:=0.5;

e := .5

>    for i to n do
x:=2*(i-n/2-1/2)/n;
for j to n do
y:=2*(n/2-j+1/2)/n;
P[i,j]:=0;
if evalf(sqrt(x^2+y^2)) < 0.7^2  then P[i,j]:=exp(-(x^2+y^2/e^2)) fi:
od:od:

振幅をプロットしてみましょう。

>    listplot3d(P,shading=zhue);

[Maple Plot]

>    listplot3d(P,shading=zhue,orientation=[0,0],style=patchnogrid);

[Maple Plot]

これをフーリエ変換してみましょう。

>    Q:=FFT2(n,P):

>    R:=get_intensity(n,Q):

>    S:=FFTshift(n,R):

>    listplot3d(S,axes=boxed);

[Maple Plot]

対数スケールでプロットしてみましょう。裾野の部分がよく分かります。

>    listplot3d([seq([seq(log(S[i,j]+0.0001),i=2..n)],j=2..n)],shading=zhue,view=-10..12,orientation=[0,0]);

[Maple Plot]

回転対称でなく横方向と縦方向の広がり方が異なっています。

>    listplot3d([seq([seq(log(S[i,j]),i=n/2-5..n/2+7)],j=n/2-5..n/2+7)],shading=zhue);

[Maple Plot]

2つのスリットの干渉

2つのスリットがあった場合それぞれを通過する波同士が干渉しあいます。まずふたつのスリットを表現する瞳関数を作成します

>    n:=32:P:=matrix(n,n):

>    for i to n do
x:=2*(i-n/2-1/2)/n;
for j to n do
y:=2*(n/2-j+1/2)/n;
P[i,j]:=0;
if  abs(x)< .5  and abs(y-0.5)<0.05  then P[i,j]:=1 fi:
if  abs(x)< .5  and abs(y+0.5)<0.05  then P[i,j]:=1 fi:
od:od:

>    listplot3d(P);

[Maple Plot]

>    listplot3d(P,orientation=[0,0],shading=zgrayscale,style=patchnogrid);

[Maple Plot]

このフーリエ変換を行ってみましょう。

>    Q:=FFT2(n,P):

>    R:=get_intensity(n,Q):

>    S:=FFTshift(n,R):
listplot3d(S,shading=zgrayscale,orientation=[0,0],style=patchnogrid);

[Maple Plot]

このように干渉が観測できます。

>    listplot3d(S,shading=zhue,orientation=[20,20]);

[Maple Plot]

Gerchberg & Saxton の位相再生アルゴリズム

このように瞳関数をフーリエ変換すれば点像関数が得られるのならば、逆フーリエ変換を行えば点像関数が特定のパターンになるように瞳関数を構成できます。振幅は一様として固定し位相だけを再生することができます。これはGerchberg&Saxtonの位相再生アルゴリズムと呼ばれています。この手順を示します。

これにより位相データが収束しそのフーリエ変換が目的のパターンに近づいていきます。

目的のパターンを用意します。

>    n:=32;Summ:=0:

n := 32

>    target:=matrix(n,n);

target := array(1 .. 32,1 .. 32,[])

>    for i to n do
x:=2*(i-n/2-1/2)/n;
for j to n do
y:=2*(n/2-j+1/2)/n;
target[i,j]:=0;
if ((abs(x)< .8 and abs(y) < .17) or (abs(y)<0.8 and abs(x)<.17)) then target[i,j]:=1;Summ:=Summ+2 fi:
od:od:

和 Summ も求めておきます。

>    Summ;

552

目標ののパターンをプロットします。

>    listplot3d(target,shading=zgrayscale,orientation=[5,5]);

[Maple Plot]

位相の行列を用意し乱数(0から2π)をセットします。

>    object:=matrix(n,n):

>    for i to n do
for j to n do
phase:=evalf(2*Pi*rand()/10^12);
object[i,j]:=evalf(cos(phase))+I*evalf(sin(phase));
od:od:

位相を求める関数を用意します。

>    get_phase:=proc(n,A)local B,i,j:
B:=matrix(n,n):
for i to n do
for j to n do
B[i,j]:=evalf(arctan(Re(A[i,j]),Im(A[i,j]))):
od:od:
RETURN(B):
end:

初期データ(object)の位相をプロットしてみましょう。

>    C:=get_phase(n,object):

>    listplot3d(C,orientation=[0,0],shading=zhue);

[Maple Plot]

プロットアニメーション用のデータを初期化しておきます。

>    p1:=NULL:p2:=NULL:

さきほどの手順(Gerchberg&Saxtonアルゴリズム)を10回繰り返します。

>    for k to 10 do
result:=FFTshift(n,FFT2(n,object)):
B:=get_intensity(n,result):
p1:=p1,listplot3d(B,orientation=[10,10],shading=zhue);

for i to n do
for j to n do
amp:=Summ*target[i,j]:
if amp <> 0 then
temp:=sqrt(Re(result[i,j])^2+Im(result[i,j])^2);
if temp > 0 then
object[i,j]:=amp*result[i,j]/temp;
else
object[i,j]:=amp;
fi:fi:
od:od:
result:=invFFT2(n,FFTshift(n,object)):

for i to n do
for j to n do
temp:=sqrt(Re(result[i,j])^2+Im(result[i,j])^2);
if temp >0 then object[i,j]:=result[i,j]/temp;
else
object[i,j]:=1;
fi:
od:od:
C:=get_phase(n,object);
p2:=p2,listplot3d(C,orientation=[0,0],shading=zhue);
od:

フーリエ変換(点像関数)をアニメーションでみてみましょう。

>    display(p1,insequence=true);

[Maple Plot]

位相データの収束をアニメーションでみてみましょう。

>    display(p2,insequence=true);

[Maple Plot]

>    listplot3d(C,axes=boxed,orientation=[10,10]);

[Maple Plot]

この位相データはファイルに出力できます。ここでは弊社光学設計プログラムCODEVの位相データに

出力してみました。

>    fd := fopen("pha.1.int", WRITE):
fprintf(fd,"This is a phase file\n"):
fprintf(fd,"GRD %d %d  WFR  WVL 0.5000  SSZ   1000   NDA -32767\n",n,n):
for i from 1 to n do
for j from 1 to n do
fprintf(fd," %d",trunc(1000*C[i,j]/(2*Pi)));
od:fprintf(fd,"\n"):od:
fclose(fd);

>    B:=get_intensity(n,FFTshift(n,FFT2(n,object)));

B := B

>    listplot3d(B,shading=zgrayscale,orientation=[10,10],style=patchnogrid);

[Maple Plot]

>