フーリエ光学と4f系 — レンズは何をフーリエ変換しているのか

顕微鏡で無色透明な細胞を覗いても、ふつうは何も見えません。細胞は光をほとんど吸収しないので、通り抜けた光の「強さ」は背景とほぼ同じだからです。ところが位相差顕微鏡という装置を使うと、同じ細胞が突然くっきりと浮かび上がります。試料には一切触れていません。変えたのは、レンズの後ろの焦点面にほんの小さな「位相板」を1枚置いただけです。

なぜ焦点面に部品を1枚置くだけで、見えなかったものが見えるようになるのでしょうか。答えは、レンズの後側焦点面には、入力像の「フーリエ変換」がそのまま物理的に並んでいるからです。焦点面の位置は空間周波数に対応しており、そこに置いたマスクは、そのまま画像処理でいうフィルタのカーネルとして働きます。デジタル画像処理でやっている「FFT → 周波数マスク → 逆FFT」を、光は電子回路もCPUも使わず、レンズ2枚とマスク1枚で、しかも光速で実行してしまうわけです。

この「レンズ＝フーリエ変換器」という視点はフーリエ光学と呼ばれ、応用は非常に広い分野に及びます。レーザー加工装置では、ピンホール1つを焦点面に置くだけでビームの空間ノイズを除去します（スペイシャルフィルタ／ビームクリーンアップ）。半導体露光装置では、投影レンズの瞳を伝達関数として設計し、位相シフトマスクや輪帯照明で解像限界を押し広げます。流体実験のシュリーレン装置は、焦点にナイフエッジを差し込むだけで空気の密度勾配を可視化します。さらに天体望遠鏡の補償光学、光相関器によるパターン認識、位相回復（ptychography）など、現代の光計測はほぼすべてこの枠組みの上に立っています。

本記事では、「なぜレンズがフーリエ変換をするのか」を回折積分から一行ずつ計算して証明し、その帰結として4f系の空間フィルタリングを導きます。最後にPythonで、遠方場の解析解との照合、ローパス／ハイパス／ナイフエッジ／位相コントラストの再現までを行います。

これが本記事のゴールである4f系の全体像です。入力面から数えて $f$ ごとに、レンズ・フーリエ面・レンズ・出力面が並び、全長がちょうど $4f$ になっています。前半の $2f$ でレンズ $L_1$ が入力のフーリエ変換を中央の面に作り、後半の $2f$ でレンズ $L_2$ がそれを像に戻す。したがって中央のフーリエ面に置いた黒いマスク1枚が、そのまま「周波数領域での掛け算」として働きます。以降の議論は、この絵の各段階を数式で裏付けていく作業だと思ってください。

本記事の内容

• レンズがフーリエ変換器になる理由 — フレネル回折積分からの厳密な導出
• 焦点面の座標と空間周波数の対応（スケーリング則）
• 4f系の伝達関数と「出力＝入力とマスクのインパルス応答の畳み込み」の導出
• Pythonによる遠方場の検証と、4種類の空間フィルタの再現

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 波動光学とは？光を波として扱う考え方の基礎
• 単スリットのフラウンホーファー回折
• フーリエ変換とは？定義から性質まで
• 畳み込み定理をわかりやすく解説

複素振幅（フェーザ）の表記と、フーリエ変換の平行移動・畳み込みの性質を使います。回折の基本的な考え方（ホイヘンス・フレネルの原理）を押さえていれば十分です。

レンズは何をしているのか — まず直感から

数式に入る前に、幾何光学の知識だけで「レンズがフーリエ変換をする」ことの気配を掴んでおきましょう。

高校物理で習ったレンズの性質を思い出してください。平行光線は、後側焦点面の1点に集まります。しかもその点の位置は、光線が「どの方向から来たか」だけで決まり、レンズのどこを通ったかには依存しません。真っ直ぐ光軸に平行に入った光線は焦点の中心へ、光軸から角度 $\theta$ 傾いて入った平行光線束は、焦点面上の $x = f\tan\theta \approx f\theta$ の位置へ集まります。

つまりレンズは、「方向」を「位置」に変換する装置です。入射側では方向という角度情報だったものが、焦点面では位置という空間情報に化けている。これが第一のポイントです。

3色の平行光線束が、それぞれ焦点面の別々の1点に集まっていることに注目してください。同じ色の光線はレンズの上端も下端も通っていますが、行き先は完全に一致します。つまり集光点の位置を決めているのは「レンズのどこを通ったか」ではなく「どの方向から来たか」だけです。角度が2倍になれば焦点面での高さも2倍になり、$x = f\theta$ という比例関係が読み取れます。この「角度 → 位置」の一対一対応こそが、あとで空間周波数とフーリエ面の座標を結ぶ土台になります。

第二のポイントは波動光学から来ます。ある平面上の複素振幅分布 $U(\xi,\eta)$ は、さまざまな方向に進む平面波の重ね合わせとして書けます。これを角スペクトル分解と呼びます。方向余弦 $(\alpha,\beta)$ の平面波の成分の強さは、ちょうど $U$ のフーリエ変換 $\hat U(f_X,f_Y)$ を $f_X=\alpha/\lambda$, $f_Y=\beta/\lambda$ で評価したものになります。フーリエ変換とは要するに「この分布はどの方向へ進む平面波をどれだけ含んでいるか」の目録なのです。

この2つを繋ぐと、答えが見えてきます。

1. 入力面の複素振幅は、方向ごとの平面波成分に分解できる。その重みがフーリエ変換である。
2. レンズは、方向 $\theta$ の平面波を焦点面の位置 $x=f\theta$ に集める。
3. したがって焦点面には、方向ごとの重み＝フーリエ変換が、位置の関数として並ぶ。

言い換えると、入力面の細かい模様（高い空間周波数）は大きく回折して大きな角度で進むので、焦点面では外側に落ちる。粗い模様（低い空間周波数）は真っ直ぐ進むので中心付近に落ちる。焦点面の中心は「べったり一様な成分（直流）」、外側へ行くほど「細かい縞模様の成分」に対応しているわけです。

ここまでは幾何光学的な当たり付けにすぎません。本当に嬉しいのは、この対応が振幅だけでなく位相まで含めて、厳密なフーリエ変換になっているという点です。そしてそれは、レンズを置く位置を「入力面から距離 $f$」に選んだときにだけ成立します。この「なぜ $f$ でなければならないのか」を明らかにするために、次節でフレネル回折積分という道具を準備します。

準備1：スカラー回折とフレネル回折積分

光の伝搬を扱う土台を作ります。単色光（角周波数 $\omega$）を仮定し、電場の1成分を複素振幅 $U(x,y,z)$ で表します。実際の場は $\mathrm{Re}\{U e^{-j\omega t}\}$ です。この記事では光学分野の標準（Goodman）に合わせ、時間依存を $e^{-j\omega t}$ にとります。この約束のもとでは、原点から広がっていく球面波は $e^{jkr}/r$ と書けます（$k=2\pi/\lambda$）。符号の約束は最後まで一貫させないと、レンズの位相の符号が逆になってしまうので注意してください。

真空中（あるいは一様媒質中）で $U$ はヘルムホルツ方程式

$$ \begin{equation} (\nabla^2 + k^2)U = 0 \end{equation} $$

を満たします。これを平面 $z=0$ の境界値 $U(\xi,\eta,0)$ から解いたのがレイリー・ゾンマーフェルトの回折公式です。

$$ \begin{equation} U(x,y,z) = \frac{1}{j\lambda}\iint U(\xi,\eta,0)\,\frac{e^{jkr}}{r}\cos\theta\;d\xi\,d\eta, \qquad r = \sqrt{z^2 + (x-\xi)^2 + (y-\eta)^2} \end{equation} $$

この式は「境界面の各点が球面波 $e^{jkr}/r$ を出し、それを全部足し合わせる」というホイヘンス・フレネルの原理を、そのまま数式にしたものです。$\cos\theta = z/r$ は傾斜係数で、真横方向に出る波を抑える因子、$1/(j\lambda)$ は振幅と位相の規格化です。

近軸近似でフレネル積分に落とす

このままでは平方根が指数の中にあって扱えません。そこで近軸近似（観測点が光軸の近くにあり、$z$ が横方向のずれよりずっと大きい）を入れます。

まず $r$ を展開します。$\rho^2 = (x-\xi)^2+(y-\eta)^2$ とおくと

$$ r = z\sqrt{1 + \frac{\rho^2}{z^2}} \approx z\left(1 + \frac{\rho^2}{2z^2} – \frac{\rho^4}{8z^4}+\cdots\right) $$

ここで扱いが2段階に分かれます。分母の $r$ は振幅を決めるだけなので、いちばん粗い近似 $r\approx z$ で十分です。一方指数の中の $kr$ は位相を決めるので、$\rho^2$ の項まで残さなければいけません。位相は $2\pi$ で一周するので、$k$ 倍された微小量が無視できるとは限らないからです。同じ理由で $\cos\theta\approx 1$ としてよいのは振幅だからです。

第3項 $-k\rho^4/(8z^3)$ を落としてよい条件が、フレネル近似の妥当性条件

$$ z^3 \gg \frac{k}{8}\left[(x-\xi)^2+(y-\eta)^2\right]^2_{\max} $$

です。以上をまとめると、フレネル回折積分が得られます。

$$ \begin{equation} U_2(x,y) = \frac{e^{jkz}}{j\lambda z}\iint U_1(\xi,\eta)\,\exp\!\left\{\frac{jk}{2z}\left[(x-\xi)^2 + (y-\eta)^2\right]\right\}d\xi\,d\eta \end{equation} $$

この式の構造を見てください。積分の中身は $(x-\xi)$ の関数だけになっています。つまりフレネル伝搬は、チャープ関数（二次位相因子）$\exp[jk\rho^2/(2z)]$ との畳み込みです。自由空間を距離 $z$ 進むことは、線形不変系を1つ通すことに等しい。この見方は後で効いてきます。

フラウンホーファー近似 — 遠方場がフーリエ変換になる

指数の中の二乗を展開してみましょう。

$$ (x-\xi)^2+(y-\eta)^2 = (x^2+y^2) – 2(x\xi+y\eta) + (\xi^2+\eta^2) $$

これを代入すると、$x,y$ だけの項は積分の外に出せます。

$$ \begin{equation} U_2(x,y) = \frac{e^{jkz}}{j\lambda z}e^{\frac{jk}{2z}(x^2+y^2)}\iint U_1(\xi,\eta)\,e^{\frac{jk}{2z}(\xi^2+\eta^2)}\,e^{-j\frac{2\pi}{\lambda z}(x\xi+y\eta)}\,d\xi\,d\eta \end{equation} $$

積分の中に $e^{jk(\xi^2+\eta^2)/2z}$ という余計な因子が残っているせいで、これはまだフーリエ変換ではありません。しかし距離 $z$ を十分大きくして

$$ \begin{equation} z \gg \frac{k(\xi^2+\eta^2)_{\max}}{2} = \frac{\pi D^2}{4\lambda} \end{equation} $$

（$D$ は入力の広がり）とすると、この因子の位相は入力の全域で $1$ ラジアンにも満たなくなり、$e^{jk(\xi^2+\eta^2)/2z}\approx 1$ と置けます。これがフラウンホーファー近似（遠方場近似）で、このとき

$$ \begin{equation} U_2(x,y) = \frac{e^{jkz}}{j\lambda z}\,e^{\frac{jk}{2z}(x^2+y^2)}\;\hat U_1\!\left(\frac{x}{\lambda z},\,\frac{y}{\lambda z}\right) \end{equation} $$

となり、確かに入力のフーリエ変換が現れます。ただし前に $e^{jk(x^2+y^2)/2z}$ という残留二次位相が付いていることに注意してください。強度 $|U_2|^2$ だけを測るなら気になりませんが、複素振幅としては純粋なフーリエ変換ではありません。

そして何より、この近似には実用上の致命傷があります。開口 $D=1$ mm、波長 $\lambda=633$ nm で必要な距離を見積もると

$$ z \gg \frac{\pi (10^{-3})^2}{4\times 633\times 10^{-9}} \approx 1.2\ \text{m} $$

つまり実際には十数メートル離れないと遠方場になりません。しかも必要距離は $D^2$ に比例するので、開口を10 mmにすれば100倍の約124 m。これでは実験室に置けません。

左のグラフは、積分の中に残る邪魔な二次位相が距離とともにどう小さくなるかを示しています。$z=0.05$ m では開口の端で24.8 rad もあり、$e^{jk\xi^2/2z}\approx1$ とは到底言えません。$z=1.2$ m でようやく端が 1.03 rad まで下がり、破線の目安に届きます。右のグラフはその必要距離を開口 $D$ の関数で描いたもので、$D^2$ に比例して急激に伸びるのが分かります。$D=1$ mm なら1.2 mですが、$D=10$ mm では124 m。光学定盤の長さ（数メートル）をあっさり超えてしまいます。遠方場は「遠ざかれば得られる」が「遠ざかること自体が不可能」なのです。

ここで自然な問いが立ちます。「遠くまで行かないと得られないフーリエ変換を、手元で作れないか」。答えがレンズです。レンズは「二次位相を掛ける素子」であり、フレネル伝搬が生む二次位相をちょうど打ち消してくれます。次節で、レンズが本当に二次位相因子として書けることを確かめましょう。

準備2：薄肉レンズは二次位相板である

レンズを「光線を曲げる部品」ではなく、「波面に位相をつける板」として捉え直します。

イメージとしてはこうです。ガラスの中では光の速度が $c/n$ に落ちるので、ガラスが厚い場所ほど波は「遅れます」。凸レンズは中心が厚く周辺が薄いので、中心を通った波ほど大きく遅れる。結果として、入ってきた平面波の中心部が引き止められ、周辺部が先に進み、波面が凹んで（＝焦点に向かって収束する球面波になって）出ていく。これが結像の波動光学的な描像です。

これを式にします。レンズを、光線が横方向にずれない程度に薄いとみなす近似（薄肉近似）のもとで、位置 $(x,y)$ での厚さを $\Delta(x,y)$、レンズ全体の最大厚さを $\Delta_0$ とします。この点を通る光が受け取る位相は、「ガラスの中を $\Delta$ だけ進んだ分」と「空気の中を $\Delta_0-\Delta$ だけ進んだ分」の和です。

$$ \phi(x,y) = kn\Delta(x,y) + k\left[\Delta_0 – \Delta(x,y)\right] $$

したがってレンズの複素透過率は

$$ \begin{equation} t_l(x,y) = e^{jk\Delta_0}\,e^{jk(n-1)\Delta(x,y)} \end{equation} $$

次に厚さ $\Delta(x,y)$ を、2つの球面の半径 $R_1,R_2$（符号は光の進行方向を正とする通常の約束）で書きます。第1面のサグ（球面の落ち込み）は $R_1 – \sqrt{R_1^2-(x^2+y^2)}$ ですから

$$ \Delta(x,y) = \Delta_0 – R_1\!\left(1-\sqrt{1-\frac{x^2+y^2}{R_1^2}}\right) + R_2\!\left(1-\sqrt{1-\frac{x^2+y^2}{R_2^2}}\right) $$

ここで近軸近似 $\sqrt{1-a}\approx 1-a/2$ を使います。すると各サグは

$$ R\left(1-\sqrt{1-\frac{x^2+y^2}{R^2}}\right) \approx \frac{x^2+y^2}{2R} $$

と、きれいに二次式になります。これを代入すると

$$ \Delta(x,y) \approx \Delta_0 – \frac{x^2+y^2}{2}\left(\frac{1}{R_1}-\frac{1}{R_2}\right) $$

これを $t_l$ の式に戻します。

$$ t_l(x,y) = e^{jkn\Delta_0}\exp\!\left[-jk(n-1)\frac{x^2+y^2}{2}\left(\frac{1}{R_1}-\frac{1}{R_2}\right)\right] $$

最後にレンズメーカーの公式

$$ \frac{1}{f} = (n-1)\left(\frac{1}{R_1}-\frac{1}{R_2}\right) $$

を代入すれば、定数位相 $e^{jkn\Delta_0}$（以降は無視します）を除いて

$$ \begin{equation} \boxed{\;t_l(x,y) = \exp\!\left[-\frac{jk}{2f}\left(x^2+y^2\right)\right]\;} \end{equation} $$

を得ます。薄肉レンズとは、負の符号を持つ二次位相因子そのものです。

符号の意味を確認しておきましょう。平面波 $U=1$ を入れると、出た直後の場は $e^{-jk(x^2+y^2)/2f}$ です。これは焦点 $z=f$ に向かって収束する球面波 $\propto e^{-jk\rho^2/2f}$ の近軸表現になっています（発散球面波なら $+$ 符号）。凸レンズ（$f>0$）が平行光を焦点に集める、という当たり前の事実がちゃんと再現されています。

そして重要なのは、フレネル伝搬が持つチャープ $\exp[+jk\rho^2/2z]$ と、レンズの位相 $\exp[-jk\rho^2/2f]$ が符号が逆だという点です。両者は打ち消し合える。この打ち消しが完璧に起きる配置を見つけるのが、次節の仕事です。

左のパネルでは、球面で作った実際のレンズ断面（青の塗り）に近軸近似の放物面（赤破線）が重なっており、$\Delta(x)$ を二次式で置き換える近似が中心付近で十分よいことが確認できます。中央のパネルが本節の要点で、レンズの位相（下向きの放物線）と距離 $f$ の伝搬チャープ（上向きの放物線）は $\pm2000$ rad を超える巨大な量なのに、足すと厳密にゼロの直線になります。右のパネルはその物理的な帰結で、平面だった等位相面がレンズを境に焦点を中心とする球面へ変わっています。レンズは光線を曲げているのではなく、波面に二次位相を書き込んでいるという見方が、この3枚で腑に落ちるはずです。

本題：レンズがフーリエ変換器になることの証明

いよいよ主定理です。証明する内容を先に宣言しておきます。

定理　入力の複素振幅 $U_o(\xi,\eta)$ を焦点距離 $f$ のレンズの前側焦点面（レンズから距離 $f$ 手前）に置く。このとき後側焦点面（レンズから距離 $f$ 後ろ）の複素振幅は、定数因子を除いて $U_o$ の厳密な2次元フーリエ変換になる。残留二次位相は現れない。

計算は3ステップです。

• ステップ1：入力面からレンズ直前まで、距離 $f$ のフレネル伝搬
• ステップ2：レンズの位相因子を掛ける
• ステップ3：レンズ直後から後側焦点面まで、距離 $f$ のフレネル伝搬

ステップ1：入力面 → レンズ直前

フレネル回折積分に $z=f$ を入れるだけです。レンズ面の座標を $(u,v)$ とします。

$$ U_l^-(u,v) = \frac{e^{jkf}}{j\lambda f}\iint U_o(\xi,\eta)\exp\!\left\{\frac{jk}{2f}\left[(u-\xi)^2+(v-\eta)^2\right]\right\}d\xi\,d\eta $$

ステップ2：レンズを通す

前節で求めた透過率を掛けます。

$$ U_l^+(u,v) = U_l^-(u,v)\cdot\exp\!\left[-\frac{jk}{2f}(u^2+v^2)\right] $$

ステップ3：レンズ直後 → 後側焦点面

もう一度フレネル伝搬です。出力面の座標を $(x,y)$ とします。

$$ U_f(x,y) = \frac{e^{jkf}}{j\lambda f}\iint U_l^+(u,v)\exp\!\left\{\frac{jk}{2f}\left[(x-u)^2+(y-v)^2\right]\right\}du\,dv $$

3つをまとめる

以上を1本の4重積分に書き下します。

$$ U_f(x,y) = \frac{e^{j2kf}}{(j\lambda f)^2}\iint\!\!\iint U_o(\xi,\eta)\, \exp\!\left\{\frac{jk}{2f}\Big[\big\{(u-\xi)^2 – u^2 + (x-u)^2\big\} + \big\{(v-\eta)^2 – v^2 + (y-v)^2\big\}\Big]\right\}d\xi\,d\eta\,du\,dv $$

中括弧の第1項が $x$ 方向、第2項が $y$ 方向の位相です。それぞれ「入力面からの伝搬」「レンズの位相（マイナス）」「焦点面までの伝搬」の3つが並んでいることを確認してください。$x$ 方向と $y$ 方向は完全に分離しているので、以降は $x$ 方向だけを追います。カギは中括弧の中の整理です。まず展開します。

$$ (u-\xi)^2 – u^2 + (x-u)^2 = \left(u^2 – 2u\xi + \xi^2\right) – u^2 + \left(x^2 – 2xu + u^2\right) $$

$u^2$ が $+1$ 個、$-1$ 個、$+1$ 個で、合計 $u^2$ が1つだけ残ります。レンズが二次位相を1つ食べてくれたわけです。整理すると

$$ = u^2 – 2u(\xi + x) + \xi^2 + x^2 $$

次に $u$ について平方完成します。$u^2-2u(\xi+x) = \left[u-(\xi+x)\right]^2 – (\xi+x)^2$ ですから

$$ = \left[u-(\xi+x)\right]^2 – (\xi+x)^2 + \xi^2 + x^2 $$

ここで $(\xi+x)^2 = \xi^2 + 2x\xi + x^2$ を代入すると、$\xi^2$ と $x^2$ が両方きれいに消えます。

$$ \begin{equation} (u-\xi)^2 – u^2 + (x-u)^2 = \left[u-(\xi+x)\right]^2 – 2x\xi \end{equation} $$

左のグラフでは、伝搬・レンズ・伝搬の3つの二次位相（青・赤・橙）を足した緑の曲線が、恒等式の右辺である黒破線とぴったり重なっています。数値的な最大差は $1.1\times10^{-14}$ で、これは倍精度の丸め誤差そのものです。右のグラフは、そこから完全平方を引いた「残りかす」が $x$ について直線になることを示しており、傾きは $-2\xi$。$\xi$ を 0.4、0.8、1.6 と変えると傾きも比例して急になります。つまり残った位相は $x$ と $\xi$ の積だけであり、これがフーリエ核の指数部そのものです。

これが本記事の心臓部です。 $\xi^2$（入力面の二次位相）も $x^2$（出力面の残留二次位相）も、跡形もなく消えました。残ったのは $u$ に関する完全平方と、フーリエ核そのものである $-2x\xi$ だけです。前節のフラウンホーファー近似では $x^2$ の項が残ってしまいましたが、レンズを使うとそれが消える。しかも $\xi^2$ が消えたということは、フラウンホーファー近似で必要だった「遠く離れる」という条件が不要になったことを意味します。距離ではなく、レンズの位相が仕事をしているのです。

さて、$u$ に関する積分を実行します。$w = u-(\xi+x)$ と置くと

$$ \int_{-\infty}^{\infty}\exp\!\left[\frac{jk}{2f}w^2\right]dw $$

これはフレネル積分（複素ガウス積分）で、公式 $\int_{-\infty}^{\infty}e^{jat^2}dt = \sqrt{\pi/a}\;e^{j\pi/4}$（$a>0$）を $a=k/(2f)$ に適用します。$\pi/a = 2\pi f/k = \lambda f$ ですから

$$ \int_{-\infty}^{\infty}\exp\!\left[\frac{jk}{2f}w^2\right]dw = \sqrt{\lambda f}\;e^{j\pi/4} $$

$y$ 方向も同じなので、$u,v$ 両方の積分から出てくる因子は

$$ \left(\sqrt{\lambda f}\,e^{j\pi/4}\right)^2 = \lambda f\, e^{j\pi/2} = j\lambda f $$

（実際にはレンズ開口は有限なので積分範囲も有限ですが、開口が入力の遠方場を十分に受け止める大きさなら、この無限積分の値でよく近似できます。有限開口の効果は後で「伝達関数」として扱います。）

前係数をまとめます。$(j\lambda f)^2 = -\lambda^2f^2$ なので

$$ \frac{e^{j2kf}}{(j\lambda f)^2}\cdot j\lambda f = \frac{j\,e^{j2kf}}{-\lambda f} = \frac{-j\,e^{j2kf}}{\lambda f} = \frac{e^{j2kf}}{j\lambda f} $$

最後の等号は $1/j = -j$ を使いました。以上をすべて合わせると、目的の結果が得られます。

$$ \begin{equation} \boxed{\;U_f(x,y) = \frac{e^{j2kf}}{j\lambda f}\iint U_o(\xi,\eta)\,\exp\!\left[-j\frac{2\pi}{\lambda f}\left(x\xi + y\eta\right)\right]d\xi\,d\eta\;} \end{equation} $$

$\exp[jk(-2x\xi)/2f] = \exp[-j k x\xi/f] = \exp[-j2\pi x\xi/(\lambda f)]$ と書き換えたことに注意してください。これは紛れもなくフーリエ変換の定義式です。前係数 $e^{j2kf}/(j\lambda f)$ は $x,y$ に依存しない定数なので、位相まで含めて純粋なフーリエ変換だと言えます。

レンズを別の位置に置いたらどうなるか

$d=f$ という配置が特別だったことを確認するため、入力面をレンズから距離 $d$ 手前に置いた一般の場合を考えます。同じ計算を $d$ で行うと、$u^2$ の係数が $\frac{1}{d}-\frac{1}{f}+\frac{1}{f}=\frac{1}{d}$ ではなく、正確には

$$ \frac{(u-\xi)^2}{d} – \frac{u^2}{f} + \frac{(x-u)^2}{f} $$

という形になり、平方完成の結果に $x^2$ の項が残ります。結論だけ書くと

$$ \begin{equation} U_f(x,y) = \frac{e^{jk(d+f)}}{j\lambda f}\,\exp\!\left[\frac{jk}{2f}\left(1-\frac{d}{f}\right)(x^2+y^2)\right]\hat U_o\!\left(\frac{x}{\lambda f},\frac{y}{\lambda f}\right) \end{equation} $$

となります。読み取れることは3つです。

1. 振幅は $d$ によらない。$|U_f| \propto |\hat U_o|$ は、レンズをどこに置いても成り立ちます。強度だけ測る回折実験なら $d$ は自由。
2. 位相は $d$ に依存する。残留二次位相の係数は $\frac{k}{2f}\left(1-\frac{d}{f}\right)$ で、$d=0$（物体をレンズに密着）のとき最大 $\frac{k}{2f}$、$d=f$ でちょうどゼロになります。
3. したがって、複素振幅を扱う光学系（＝これから作る4f系）では $d=f$ が必須です。位相が歪んでいると、次のレンズで逆変換したときに像がぼける。

この $d$ 依存性は後でPythonで数値的に確かめます。ここまでで「レンズ＝フーリエ変換器」は証明できました。次は、そのフーリエ変換の目盛り、つまり焦点面のどこが何 lp/mm に対応するのかをはっきりさせましょう。

スケーリング則：焦点面の座標 = 空間周波数

導いた式のフーリエ核を、もう一度睨みます。

$$ \exp\!\left[-j2\pi\left(\frac{x}{\lambda f}\,\xi + \frac{y}{\lambda f}\,\eta\right)\right] $$

標準的なフーリエ変換の核は $\exp[-j2\pi(f_X\xi + f_Y\eta)]$ ですから、比較すれば対応関係は一目瞭然です。

$$ \begin{equation} f_X = \frac{x}{\lambda f},\qquad f_Y = \frac{y}{\lambda f} \end{equation} $$

つまり焦点面の物理座標 $x$ [m] と空間周波数 $f_X$ [1/m] は、$\lambda f$ という一つの定数で結ばれています。$\lambda f$ は長さの2乗の次元を持ち、フーリエ光学の「換算係数」として繰り返し登場します。

具体的な数値を入れてみます。He-Neレーザー $\lambda = 633$ nm、焦点距離 $f = 100$ mm のレンズなら

$$ \lambda f = 633\times10^{-9}\times 0.1 = 6.33\times10^{-8}\ \text{m}^2 $$

このとき、入力面にある 1 lp/mm（1 mmに1本の縞、$f_X = 1000$ 1/m）の格子は、焦点面の

$$ x = \lambda f\, f_X = 6.33\times10^{-8}\times1000 = 6.33\times10^{-5}\,\text{m} = 0.0633\ \text{mm} $$

に回折光として現れます。10 lp/mm なら 0.633 mm、100 lp/mm なら 6.33 mm。細かい模様ほど外側という直感どおりです。

左の3行は、同じ光学系（$\lambda=633$ nm、$f=100$ mm）に周期の違う格子を入れたときの焦点面の様子です。周期200 µm（5 lp/mm）の粗い格子では±1次が中心から0.32 mmしか離れませんが、周期20 µm（50 lp/mm）まで細かくすると3.17 mmまで飛び出します。右のグラフはこれを直線 $x = \lambda f\,f_X$ として一本にまとめたもので、3つの格子の実測位置がすべてこの直線上に乗ります。半径0.5 mmのピンホールを置くことは、この直線を $x=0.5$ mm で切ることに等しく、7.9 lp/mm という遮断周波数が幾何学的に読み取れるわけです。

逆に読むこともできます。焦点面に半径 $a$ のピンホールを置いたとすると、通過できる最大空間周波数は

$$ \begin{equation} f_c = \frac{a}{\lambda f} \end{equation} $$

です。$f=100$ mm、$a=0.5$ mm のピンホールなら $f_c = 0.5\times10^{-3}/6.33\times10^{-8} \approx 7900$ 1/m $= 7.9$ lp/mm。これは「1 mmあたり約8本より細かい模様は通さない」＝「$1/f_c\approx0.13$ mmより細かい構造は消える」という意味になります。

ここで実務的に大事なのは、$\lambda f$ が小さいほど焦点面が「縮む」という点です。短焦点レンズや短波長ではスペクトルが焦点近傍に密集するので、マスクの加工精度が厳しくなります。逆に長焦点レンズを使えばスペクトルが広がって扱いやすくなりますが、装置が長くなります。実験系の設計は、この $\lambda f$ の選び方がほぼすべてを決めます。

これで「フーリエ面」に目盛りが入りました。あとはそこにマスクを置き、もう一度フーリエ変換して像に戻せば、周波数フィルタリングが完成します。それが4f系です。

4f系 — 光でできるFFTフィルタ

構成

4f系は、次のように並べた4つの区間からなります。

$$ \text{入力面}\ \xrightarrow{\ f\ }\ \text{レンズ}L_1\ \xrightarrow{\ f\ }\ \text{フーリエ面}\ \xrightarrow{\ f\ }\ \text{レンズ}L_2\ \xrightarrow{\ f\ }\ \text{出力面} $$

全長が $4f$ なので4f系（4f correlator, 4f system）と呼ばれます。ポイントは、中央のフーリエ面が「$L_1$ の後側焦点面」であると同時に「$L_2$ の前側焦点面」でもあることです。前節の定理が2回続けて使える配置になっているわけです。

出力の導出

まず $L_1$ によって、フーリエ面直前の場は

$$ U_2^-(x,y) = \frac{e^{j2kf}}{j\lambda f}\,\hat U_1\!\left(\frac{x}{\lambda f},\frac{y}{\lambda f}\right) $$

ここにマスク（瞳関数と呼びます）$P(x,y)$ を置きます。$P$ は複素透過率で、$|P|$ が振幅の減衰、$\arg P$ が位相の遅れを表します。

$$ U_2^+(x,y) = P(x,y)\,U_2^-(x,y) $$

$L_2$ でもう一度フーリエ変換します。出力面座標を $(x_3,y_3)$ として

$$ U_3(x_3,y_3) = \frac{e^{j2kf}}{j\lambda f}\iint U_2^+(x,y)\,e^{-j\frac{2\pi}{\lambda f}(x_3x+y_3y)}\,dx\,dy $$

$U_2^+$ を代入し、変数を $f_X = x/(\lambda f)$、$f_Y = y/(\lambda f)$ に取り替えます。ヤコビアンは $dx\,dy = (\lambda f)^2 df_X\,df_Y$ です。

$$ U_3(x_3,y_3) = \frac{e^{j4kf}(\lambda f)^2}{(j\lambda f)^2}\iint \hat U_1(f_X,f_Y)\,P(\lambda f f_X,\lambda f f_Y)\,e^{-j2\pi(x_3f_X+y_3f_Y)}\,df_X\,df_Y $$

前係数は $(\lambda f)^2/(j\lambda f)^2 = 1/j^2 = -1$ と簡単になります。そこでコヒーレント伝達関数を

$$ \begin{equation} H(f_X,f_Y) \equiv P(\lambda f f_X,\; \lambda f f_Y) \end{equation} $$

と定義しましょう。これは「フーリエ面に置いた物理的なマスクを、空間周波数の目盛りで読み直したもの」です。すると

$$ U_3(x_3,y_3) = -e^{j4kf}\iint \hat U_1(f_X,f_Y)H(f_X,f_Y)\,e^{-j2\pi(x_3f_X+y_3f_Y)}df_X\,df_Y $$

積分の核が $e^{-j2\pi(\cdots)}$ であることに注意してください。逆変換の核は $e^{+j2\pi(\cdots)}$ なので、これは「逆変換の座標を反転したもの」です。$(x_3,y_3)\to(-x_3,-y_3)$ と読み替えれば逆変換になります。そこで出力面の座標軸を反転して取る（これは物理的には像が上下左右反転して結ばれることを意味します）と、畳み込み定理そのものが現れます。

$$ \begin{equation} \boxed{\;U_3 = -e^{j4kf}\,\left(U_1 * h\right),\qquad h = \mathcal{F}^{-1}\{H\}\;} \end{equation} $$

出力の複素振幅は、入力とマスクのインパルス応答の畳み込みです。$h$ を振幅点像分布関数（コヒーレントPSF）と呼びます。

これでデジタル画像処理との対応が完全に付きました。

違いは2つだけです。第一に、光の演算は複素振幅に対して行われます。カメラが記録するのは強度 $|U_3|^2$ なので、非負の実数に潰されてしまう。この「最後に $|\cdot|^2$ が付く」ことが、フーリエ光学の面白さと難しさの両方を生みます。第二に、レンズの開口が有限なので、$P$ には必ず開口による窓（$|P|=0$ の外側領域）が掛かります。つまり現実の結像系はどんなに頑張っても低域通過フィルタであり、これが解像限界の正体です。

ここまでで枠組みは完成しました。次は具体的にどんな $P$ を置くとどんな絵が出るのかを見ていきます。

代表的な空間フィルタ

(1) ローパス — ピンホールによるぼかしとビームクリーンアップ

$$ P(x,y) = \begin{cases}1 & \sqrt{x^2+y^2}\le a\\ 0 & \text{otherwise}\end{cases} $$

半径 $a$ の円形ピンホールです。伝達関数は半径 $f_c=a/(\lambda f)$ の円柱関数、そのインパルス応答は

$$ h(x,y) = \pi f_c^2\,\frac{2J_1(2\pi f_c r)}{2\pi f_c r},\qquad r=\sqrt{x^2+y^2} $$

というエアリーパターン（$J_1$ は1次ベッセル関数）になります。入力像をこれで畳み込むので、$1/f_c$ より細かい構造がぼけて消えます。同時に、鋭いエッジのところにはリンギング（明暗の波打ち）が出ます。矩形の窓関数で切ったせいで、時間信号を理想LPFに通したときと同じギブス現象が起きるからです。

このフィルタの最重要の応用がスペイシャルフィルタです。レーザービームは、ミラーの傷やレンズの気泡によって、本来の滑らかなガウス分布に細かい干渉縞（空間ノイズ）を背負っています。これを対物レンズで絞り、焦点にピンホールを置くと、ノイズ成分は高い空間周波数を持つので焦点面の外側に落ち、ピンホールに遮られます。中心を通り抜けるのは滑らかな低周波成分だけ。焦点に穴の空いた板を1枚置くだけでビームが綺麗になるというのは、フーリエ光学を知らないと魔法にしか見えません。

(2) ハイパス — 暗視野照明とエッジ強調

逆に中心を小さな円板で遮ります。

$$ P(x,y) = \begin{cases}0 & \sqrt{x^2+y^2}\le b\\ 1 & \text{otherwise}\end{cases} $$

$b$ を極小にして直流成分（＝一様な背景）だけを止めると、出力は「入力から平均値を引いたもの」に近くなります。べったり明るい領域は真っ暗になり、輪郭だけが光る。これが顕微鏡の暗視野照明です。透明な微粒子でも、散乱光（＝高周波成分）だけを拾えば漆黒の背景に星のように光って見えます。

数式でも確認できます。入力を $U_1 = \bar U + \Delta U$（平均＋変動）と分けると、直流遮断は $\bar U$ を消すので $U_3 \approx \Delta U$、強度は $|\Delta U|^2$。変動が小さいときは強度も小さいですが、背景がゼロなのでコントラストは無限大になります。感度と引き換えにダイナミックレンジを稼ぐ、という設計思想です。

(3) ナイフエッジ — シュリーレンとヒルベルト変換

フーリエ面に片刃のナイフを差し込み、半平面 $f_X<0$ を遮ります。

$$ H(f_X,f_Y) = \mathrm{step}(f_X) = \frac{1}{2}\left[1+\mathrm{sgn}(f_X)\right] $$

この効果を弱位相物体で調べましょう。位相のみを変える試料（気体の密度変化、細胞、光学ガラスの歪み）は

$$ U_1 = e^{j\phi(\xi,\eta)} \approx 1 + j\phi \qquad (|\phi|\ll 1) $$

と書けます。そのまま撮ると $|U_1|^2 = 1$ で完全に不可視です。位相物体が見えない、という冒頭の問題がここに数式で現れました。

ナイフエッジを通すと何が起きるか。$\hat U_1 = \delta + j\hat\phi$ にフィルタを掛け、逆変換します。ヒルベルト変換のフーリエ対 $\mathcal F^{-1}\{-j\,\mathrm{sgn}(f_X)\hat\phi\} = \mathcal H\{\phi\}$ を使うと $\mathcal F^{-1}\{\mathrm{sgn}(f_X)\hat\phi\} = j\mathcal H\{\phi\}$ ですから

$$ U_3 = 1 + \frac{j}{2}\phi + \frac{j}{2}\cdot j\,\mathcal H\{\phi\} = 1 – \frac{1}{2}\mathcal H\{\phi\} + \frac{j}{2}\phi $$

$\phi$ の1次までとって強度を取ると

$$ \begin{equation} I = |U_3|^2 \approx 1 – \mathcal H\{\phi\} \end{equation} $$

強度が位相のヒルベルト変換に比例して変調されるわけです。ヒルベルト変換は微分に似た働き（周波数領域で $-j\,\mathrm{sgn}$ を掛ける＝各周波数を90°位相シフト）をするので、$\phi$ が階段状に変化する場所で明暗の対が現れます。物体が浮き上がって見える「レリーフ調」の像になるのはこのためです。

風洞実験のシュリーレン写真はまさにこれです。空気の密度変化が屈折率変化を生み、光路長＝位相を変える。焦点にナイフエッジを置くだけで、目に見えない衝撃波が縞として写ります。

(4) ゼルニケ位相コントラスト — DC成分だけ90°回す

ナイフエッジは片側を捨てるので光量が半減し、しかも $x$ 方向に非対称な像になります。もっと賢いやり方をゼルニケが1934年に考案しました（1953年ノーベル物理学賞）。

発想の出発点は、$U_1 \approx 1+j\phi$ の $1$（背景光、直流成分）と $j\phi$（回折光）が複素平面上で直交していることです。だから足しても長さが変わらない：$|1+j\phi|^2 = 1+\phi^2\approx1$。強度が変化しないのは、この直交性のせいです。

ならば、片方を90°回して同じ向きに揃えればいい。フーリエ面の中心（直流成分が集中する1点）にだけ $\lambda/4$ 相当の位相段差を持つ小さな板を置き、直流成分に $e^{j\pi/2}=j$ を掛けます。

$$ U_3 = j\cdot 1 + j\phi = j(1+\phi) $$

強度は

$$ \begin{equation} I = |1+\phi|^2 \approx 1 + 2\phi \end{equation} $$

位相が強度に線形変換されました。$\phi$ の1次の項が出るので、微小な位相差でもしっかり見えます。位相が大きいところが明るくなるので「ポジティブ位相コントラスト」、$-j$ を掛ければ $I\approx1-2\phi$ で「ネガティブ位相コントラスト」です。

さらに一工夫あります。位相板に減衰も持たせ、直流成分を $\alpha\,e^{j\pi/2}$（$0<\alpha<1$）にすると

$$ I = |\alpha + \phi|^2 \approx \alpha^2 + 2\alpha\phi $$

コントラスト（変動／背景）は

$$ \frac{2\alpha\phi}{\alpha^2} = \frac{2\phi}{\alpha} $$

と、$\alpha$ を小さくするほど増幅されます。実際の位相差顕微鏡の位相板は、位相リングの部分に金属薄膜を蒸着して透過率を1〜2割に落としてあります。背景を暗くして相対的に信号を浮かせる、という発想です（暗視野照明の $\alpha\to0$ 極限が暗視野そのもの、と見ることもできます）。

なお、実際の位相板は点ではなく有限の大きさのリングです。そのため位相物体の低周波成分まで一緒に位相回転されてしまい、大きな構造の内部で明るさが落ちる「シェーディングオフ」や、輪郭の外側が光る「ハロー」というアーティファクトが出ます。これは4f系の枠組みで完全に説明でき、後のPythonでも再現できます。

4種類のフィルタを見てきました。どれも「フーリエ面に何を置くか」を変えただけで、光学系の構成は同一です。ここからはPythonで実際に動かして、これらの主張を数値で検証していきます。

Pythonでの実装

以降のコードは、離散フーリエ変換を使って前節までの理論をそのまま再現します。まず共通の道具を用意します。実装で気をつけるのは3点です。第一に、np.fft.fft2 は原点が配列の左上にある約束なので、fftshift/ifftshift で中心原点に揃えます。第二に、離散フーリエ変換は「和」なので、連続フーリエ変換に合わせるには面積要素 $dx\,dy$ を掛けます。第三に、周波数軸は np.fft.fftfreq(N, d=dx) で作り、フーリエ面の物理座標は $x=\lambda f f_X$ で変換します。

import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == fnt.name for fnt in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

def FT(u, dx):
    """連続フーリエ変換の離散近似（中心原点・面積要素つき）"""
    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(u))) * dx * dx

def IFT(A, dx):
    """逆変換"""
    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(A))) / (dx * dx)

def freq_axis(N, dx):
    """中心原点の空間周波数軸 [1/m]"""
    return np.fft.fftshift(np.fft.fftfreq(N, d=dx))

実験1：方形開口の遠方場を解析解と照合する

まず、レンズの焦点面に現れる分布が本当にフーリエ変換になっているかを、解析解が分かっている方形開口で確かめます。幅 $a\times b$ の方形開口のフーリエ変換は $ab\,\mathrm{sinc}(af_X)\mathrm{sinc}(bf_Y)$（$\mathrm{sinc}(t)=\sin\pi t/\pi t$）です。

import numpy as np

N, L = 1024, 4e-3          # 格子点数と計算窓 [m]
dx = L / N
x = (np.arange(N) - N // 2) * dx
X, Y = np.meshgrid(x, x)

a, b = 0.2e-3, 0.4e-3      # 方形開口の幅 [m]
u = ((np.abs(X) <= a / 2) & (np.abs(Y) <= b / 2)).astype(complex)

A = FT(u, dx)
fx = freq_axis(N, dx)
FX, FY = np.meshgrid(fx, fx)
A_analytic = a * b * np.sinc(a * FX) * np.sinc(b * FY)

err = np.max(np.abs(np.abs(A) - np.abs(A_analytic))) / np.max(np.abs(A_analytic))
print(f"FFT と解析解 sinc の最大相対誤差: {err:.3e}")   # -> 8.340e-03
print(f"第1ゼロの理論位置 fx = 1/a = {1/a:.0f} 1/m")     # -> 5000 1/m

中央の2次元パターンは、縦横それぞれの辺の長さに対応した2つのsincの積になっています。開口が縦長（0.4 mm）横短（0.2 mm）なので、焦点面ではその逆に横方向のほうがゆったり広がっているのが分かります。右の断面比較では、FFTの数値解（実線）と解析解（破線）が5桁のダイナミックレンジにわたって重なっており、5 lp/mm ごとに並ぶゼロ点の位置まで一致しています。この対数軸での一致は、振幅の大きい主ローブだけでなく、微弱なサイドローブまで正しく再現できていることの証拠です。

最大相対誤差は約 $8.3\times10^{-3}$、つまり0.8%程度に収まりました。この誤差は物理ではなく離散化に由来します。開口の縁が格子点上にちょうど乗らないため、実効的な開口幅が1画素分だけ揺らぐからです。格子を細かくすれば単調に減ります。第1ゼロの位置も理論値 $f_X = 1/a = 5000$ 1/m と数値解が一致しました。焦点面には確かに入力のフーリエ変換が並んでいることの、最初の確認です。

実験2：振幅格子の回折次数

もう少し構造のある例として、周期 $d$、デューティ50%の振幅格子（明暗の縞）を見ます。理論では、回折次数 $m$ の振幅は矩形波のフーリエ係数 $c_m = \frac{1}{2}\mathrm{sinc}(m/2)$ で与えられ、$c_0=1/2$、$c_{\pm1}=1/\pi$、$c_{\pm2}=0$、$c_{\pm3}=1/(3\pi)$ となります。したがって強度比の予測は $I_{\pm1}/I_0 = (2/\pi)^2 \approx 0.405$、偶数次は消滅、$I_{\pm3}/I_{\pm1} = 1/9 \approx 0.111$ です。

import numpy as np

d = 50e-6                       # 格子周期 [m]  -> 20 lp/mm
N, L = 2048, 4e-3
dx = L / N
x = (np.arange(N) - N // 2) * dx
X, Y = np.meshgrid(x, x)

W = 1.0e-3                                      # 照射窓 [m]
grating = (((X / d) % 1.0) < 0.5).astype(float)  # デューティ50%の振幅格子
u = grating * ((np.abs(X) <= W/2) & (np.abs(Y) <= W/2))

A = FT(u.astype(complex), dx)
fx = freq_axis(N, dx)
I = np.abs(A[N//2, :])**2

def peak_near(f0, w=8):
    j = np.argmin(np.abs(fx - f0))
    return I[j-w:j+w+1].max()

I0, I1, I2, I3 = [peak_near(m/d) for m in (0, 1, 2, 3)]
print(f"1次/0次 = {I1/I0:.4f}  (理論 (2/pi)^2 = {(2/np.pi)**2:.4f})")
print(f"2次/0次 = {I2/I0:.2e}  (理論 0 = 偶数次は消滅)")
print(f"3次/1次 = {I3/I1:.4f}  (理論 1/9 = {1/9:.4f})")

中央のスペクトルでは、0次・±1次・±3次に鋭いピークが立つ一方、±2次（±40 lp/mm）の位置だけが $10^{-4}$ 台の谷になっています。ピークが1本の線ではなく細かい構造を持つのは、格子を1 mm の窓で切っているためで、各次数が窓関数のsincで畳み込まれた結果です。右の棒グラフでは実測（青）と理論（橙）がほぼ同じ高さで並び、2次だけが両方ともゼロに張り付いています。フーリエ級数の教科書的な性質が、そのまま焦点面の明暗として目に見えるわけです。

出力は 1次/0次 = 0.4024（理論 0.4053）、2次/0次 = 1.09e-04（理論はゼロ）、3次/1次 = 0.1116（理論 0.1111）となりました。デューティ50%の矩形波に偶数次高調波が含まれないという、フーリエ級数の基本性質が、光の回折次数の消滅として物理的に現れているわけです。焦点面の1次スポットは $x = \lambda f/d$ に立ちますから、この位置を測れば格子周期が読めます。回折格子で分光ができる理由も同じ枠組みで説明できます。

実験3：残留二次位相の $d$ 依存性を確かめる

理論の山場だった「$d=f$ でだけ二次位相が消える」を、FFTを使わずフレネル積分を直接数値積分して確かめます。入力面からレンズ面へ、レンズ面から焦点面へ、と行列積で2回積分します。

import numpy as np

lam, f = 633e-9, 20e-3
k = 2*np.pi/lam
a = 0.10e-3                                   # スリット幅 [m]
xi = np.linspace(-a/2, a/2, 801); dxi = xi[1]-xi[0]
Uo = np.ones_like(xi, dtype=complex)

u = np.linspace(-2.0e-3, 2.0e-3, 26001); du = u[1]-u[0]   # レンズ面（細かく刻む）
xout = np.linspace(-0.12e-3, 0.12e-3, 241)                # 後側焦点面

def lens_to_focus(Ul):
    K = np.exp(1j*k/(2*f)*(xout[:, None] - u[None, :])**2)
    return (K @ Ul) * du / np.sqrt(1j*lam*f)

def simulate(d):
    """入力面 -> 距離 d -> レンズ -> 距離 f -> 後側焦点面（1次元フレネル）"""
    K = np.exp(1j*k/(2*d)*(u[:, None] - xi[None, :])**2)
    Ul = (K @ Uo) * dxi / np.sqrt(1j*lam*d)     # レンズ直前
    Ul = Ul * np.exp(-1j*k/(2*f)*u**2)          # レンズの位相因子
    return lens_to_focus(Ul)

レンズ面のサンプリングが粗いと二次位相のチャープをエイリアシングしてしまうので、$2.6\times10^4$ 点という細かい刻みを使っています。続いて $d=f$、$d=0$、$d=f/2$ の3通りを比較します。

import numpy as np

def object_on_lens():
    """d=0（物体をレンズに密着）"""
    Ul = np.where(np.abs(u) <= a/2, 1.0, 0.0).astype(complex)
    return lens_to_focus(Ul * np.exp(-1j*k/(2*f)*u**2))

FT_exact = a * np.sinc(a * xout / (lam*f))       # 解析的なフーリエ変換（実数）
m = np.abs(xout) < 0.11e-3                       # 主ローブ内だけで比較

for name, U, coef in [("d = f  ", simulate(f),   0.0),
                      ("d = 0  ", object_on_lens(), 1.0),
                      ("d = f/2", simulate(f/2), 0.5)]:
    amp = np.abs(U)
    s = amp.max() / np.abs(FT_exact).max()
    e = np.max(np.abs(amp[m]/s - np.abs(FT_exact[m]))) / np.abs(FT_exact).max()
    ph = np.unwrap(np.angle(U[m])); ph -= ph[len(ph)//2]     # 中心を基準にした位相
    pred = k/(2*f)*coef*xout[m]**2; pred -= pred[len(pred)//2]
    print(f"{name}: 振幅の相対誤差 {e:.2e} / 端の位相(実測) {ph[0]:+.3f} rad "
          f"/ 予測 {pred[0]:+.3f} rad")

結果は次のとおりです。

d = f  : 振幅の相対誤差 1.46e-03 / 端の位相(実測) -0.000 rad / 予測 +0.000 rad
d = 0  : 振幅の相対誤差 1.64e-03 / 端の位相(実測) +2.948 rad / 予測 +2.948 rad
d = f/2: 振幅の相対誤差 1.43e-03 / 端の位相(実測) +1.474 rad / 予測 +1.474 rad

左のグラフでは3本の曲線が完全に重なり、解析解の $|\mathrm{sinc}|$（黒点線）とも見分けがつきません。レンズをどこに置いても、強度だけ測るなら同じ回折像が得られるということです。ところが右のグラフを見ると事情は一変します。$d=f$（緑）は横一直線でゼロですが、$d=0$（赤）は端で約3 rad も反り返り、$d=f/2$（橙）はちょうどその半分。しかもいずれも理論曲線 $\frac{k}{2f}(1-d/f)x^2$（黒破線）に重なっています。振幅図だけ見ていては絶対に気づけない違いが、位相図では一目瞭然という点が、この2枚を並べる理由です。

3点が読み取れます。第一に、振幅はどの $d$ でも解析的なフーリエ変換（sinc）に0.15%程度で一致しました。予告どおり、振幅だけなら $d$ は自由です。第二に、$d=f$ のとき位相は完全に平坦（$|$誤差$|<10^{-3}$ rad）でした。二次位相の相殺が数値でも確認できます。第三に、$d=0$ と $d=f/2$ では残留位相が現れ、その大きさは理論式 $\frac{k}{2f}(1-d/f)x^2$ の予測と小数点以下3桁まで一致します（実測 2.948 rad に対し予測 2.948 rad）。$d=f/2$ の位相が $d=0$ のちょうど半分（1.474 rad）になっているのも予測どおりです。「なぜ入力を前側焦点面に置くのか」は、この平坦性のためだったわけです。

実験4：4f系のローパスフィルタ

ここからは4f系のシミュレータを作ります。$L_1$ の変換、瞳関数の掛け算、$L_2$ の変換をまとめると、結局 IFT(FT(u)*P) の1行になります（座標反転は本質ではないので無視します）。

import numpy as np

lam, f = 633e-9, 100e-3        # He-Ne + f=100mm
N, L = 512, 8e-3
dx = L / N
x = (np.arange(N) - N // 2) * dx
X, Y = np.meshgrid(x, x)

fa = freq_axis(N, dx)
FX, FY = np.meshgrid(fa, fa)
XF, YF = lam * f * FX, lam * f * FY     # フーリエ面の物理座標 [m]
R = np.sqrt(XF**2 + YF**2)
print(f"lambda*f = {lam*f:.3e} m^2, フーリエ面の1画素 = {(XF[0,1]-XF[0,0])*1e6:.2f} um")

def four_f(u_in, P):
    """4f系: L1でFT -> 瞳関数Pを掛ける -> L2でFT（=逆変換+座標反転）"""
    return IFT(FT(u_in, dx) * P, dx)

# 振幅物体（十字＋円）
obj = np.zeros((N, N))
obj[(np.abs(X) < 1.5e-3) & (np.abs(Y) < 0.18e-3)] = 1
obj[(np.abs(Y) < 1.5e-3) & (np.abs(X) < 0.18e-3)] = 1
obj[((X - 2.2e-3)**2 + (Y - 2.2e-3)**2) < (0.8e-3)**2] = 1
u_in = obj.astype(complex)

ピンホール半径を変えながら遮断周波数と像の変化を見ます。

import numpy as np
import matplotlib.pyplot as plt

radii = [1.5e-3, 0.5e-3, 0.2e-3, 0.1e-3]
fig, axes = plt.subplots(1, 5, figsize=(18, 4))
axes[0].imshow(np.abs(u_in)**2, cmap="gray", extent=[-4, 4, -4, 4])
axes[0].set_title("入力（フィルタなし）")
for ax, a_pin in zip(axes[1:], radii):
    P = (R <= a_pin).astype(complex)
    I = np.abs(four_f(u_in, P))**2
    fc = a_pin / (lam * f)
    ax.imshow(I, cmap="gray", extent=[-4, 4, -4, 4])
    ax.set_title(f"ピンホール半径 {a_pin*1e3:.2f} mm\n"
                 f"遮断 {fc/1000:.2f} lp/mm（分解能 {1/fc*1e3:.2f} mm）")
    print(f"半径 {a_pin*1e3:.2f} mm -> fc = {fc/1000:.2f} lp/mm, "
          f"分解能 1/fc = {1/fc*1e3:.3f} mm")
for ax in axes:
    ax.set_xlabel("x [mm]")
plt.tight_layout(); plt.show()

出力は次のようになります。

半径 1.50 mm -> fc = 23.70 lp/mm, 分解能 1/fc = 0.042 mm
半径 0.50 mm -> fc = 7.90 lp/mm, 分解能 1/fc = 0.127 mm
半径 0.20 mm -> fc = 3.16 lp/mm, 分解能 1/fc = 0.317 mm
半径 0.10 mm -> fc = 1.58 lp/mm, 分解能 1/fc = 0.633 mm

上段と下段を縦に見比べると、フィルタリングの正体がよく分かります。下段のスペクトルには十字の腕に対応する縦横の直線状の成分が伸びていますが、ピンホールを絞るとその外側が緑の円の外に出て捨てられます。すると上段の像では、まず円のエッジに同心円状のリンギング（半径0.5 mm）が現れ、次に十字の輪郭が丸まり（0.2 mm）、最後には腕が太い帯に溶けてしまいます（0.1 mm）。捨てた高周波成分の量と、像のぼけ具合が一対一で対応しているのが視覚的に確認できます。

画像を並べると、ピンホールを絞るほど像がぼけ、エッジにリンギングが出る様子がはっきり見えます。注目してほしいのは、ぼけ幅が $1/f_c$ とよく対応することです。半径0.1 mmのとき分解能は0.633 mm ですが、十字の腕の幅は0.36 mm しかありません。だからこの条件では腕が完全に潰れて、太さの情報が失われます。「マスクのサイズが像の分解能を決める」というスケーリング則が、絵として確認できたわけです。実験室でスペイシャルフィルタのピンホールを選ぶときの計算も、まさにこれと同じです。

実験5：ハイパスフィルタ（暗視野・エッジ強調）

今度は中心を小さな円板で遮ります。

import numpy as np
import matplotlib.pyplot as plt

b_stop = 0.15e-3                       # 遮蔽円板の半径 [m]
P_hp = (R > b_stop).astype(complex)
I_hp = np.abs(four_f(u_in, P_hp))**2

inside = I_hp[(np.abs(X) < 0.8e-3) & (np.abs(Y) < 0.10e-3)].mean()
print(f"遮蔽の遮断周波数: {b_stop/(lam*f)/1000:.2f} lp/mm")
print(f"物体内部の平均強度 {inside:.4f} / 像の最大強度 {I_hp.max():.4f} "
      f"→ 縁/内部 の比 {I_hp.max()/inside:.1f}")

fig, ax = plt.subplots(1, 2, figsize=(10, 4.5))
ax[0].imshow(np.abs(u_in)**2, cmap="gray", extent=[-4, 4, -4, 4])
ax[0].set_title("入力（べた塗りの十字と円）")
ax[1].imshow(I_hp, cmap="gray", extent=[-4, 4, -4, 4])
ax[1].set_title(f"ハイパス後（{b_stop/(lam*f)/1000:.2f} lp/mm 以下を遮断）\n"
                "内部は暗く、輪郭だけが光る")
plt.tight_layout(); plt.show()

中央の像は、もはや「べた塗り」の面影がなく、十字と円の輪郭だけをペンでなぞったような線画になっています。右の断面プロファイルがその理由を数字で語っていて、入力（灰色）では $|x|<1.5$ mm の全域が強度1だったのに、ハイパス後（赤）は内部が0.0078まで落ち、$x=\pm1.5$ mm のエッジ位置にだけ鋭いピークが立ちます。一様な領域は「変化がない」＝「低周波しか含まない」ので、直流を捨てると丸ごと消える。エッジ検出フィルタが何をしているかを、これ以上ないほど素直に見せてくれる図です。

遮蔽の遮断周波数: 2.37 lp/mm、物体内部の平均強度 0.0078 / 像の最大強度 0.5116 → 縁/内部 の比 65.7 と出ます。入力では十字も円も一様に明るい「べた塗り」でしたが、ハイパス後は内部の明るさが約1/128に落ち、輪郭部だけが65倍以上明るく残ります。低周波成分（＝広い面積を一様に埋める成分）を捨てたので、急激に値が変わるエッジしか生き残らない。デジタル画像処理のラプラシアンフィルタやアンシャープマスクと同じことが、円板1枚で起きています。

実験6：位相物体は本当に見えないのか、そしてゼルニケ位相板

いよいよ本命です。まず、位相しか変えない試料が本当に不可視であることを確認します。

import numpy as np

phi = np.zeros((N, N))
phi[((X + 1.2e-3)**2 + Y**2) < (1.0e-3)**2] = 0.30      # 大きい円: 0.30 rad
phi[((X - 1.4e-3)**2 + (Y - 1.2e-3)**2) < (0.6e-3)**2] = 0.15   # 小さい円: 0.15 rad
u_ph = np.exp(1j * phi)                                  # 純位相物体（吸収なし）

I_none = np.abs(u_ph)**2
print(f"フィルタなし: 強度の min={I_none.min():.6f}, max={I_none.max():.6f}")
# -> min=1.000000, max=1.000000  （完全に真っ平ら＝何も見えない）

強度は最小値も最大値も 1.000000。位相物体はどう頑張っても撮像できません。カメラは複素振幅の絶対値の2乗しか記録しないので、絶対値1の位相因子は完全に情報を失います。この一行の出力が、位相差顕微鏡が必要とされた理由そのものです。

そこで、フーリエ面の中心1画素（＝直流成分）だけを $\alpha\,e^{j\pi/2} = \alpha j$ にする位相板を入れます。ここで「中心1画素」でよい理由は、計算窓の全面にわたって背景 $=1$ が広がっているためで、その離散フーリエ変換はぴたり1画素のデルタ関数になるからです。

import numpy as np

def zernike_plate(alpha, sign=+1):
    """フーリエ面の直流成分だけを alpha * exp(±j pi/2) にする位相板"""
    P = np.ones((N, N), dtype=complex)
    P[N//2, N//2] = alpha * (1j * sign)
    return P

obj_mask = ((X + 1.2e-3)**2 + Y**2) < (0.7e-3)**2                 # 大円の内側
sml_mask = ((X - 1.4e-3)**2 + (Y - 1.2e-3)**2) < (0.35e-3)**2     # 小円の内側
bg_mask  = ((X - 2.8e-3)**2 + (Y + 2.6e-3)**2) < (0.7e-3)**2      # 背景

for alpha in [1.0, 0.5, 0.3]:
    I = np.abs(four_f(u_ph, zernike_plate(alpha)))**2
    Io, Is, Ib = I[obj_mask].mean(), I[sml_mask].mean(), I[bg_mask].mean()
    print(f"alpha={alpha}: 大円(phi=0.30) {Io:.4f} [理論 {alpha**2+2*alpha*0.30:.4f}] / "
          f"小円(phi=0.15) {Is:.4f} [理論 {alpha**2+2*alpha*0.15:.4f}] / "
          f"背景 {Ib:.4f} [理論 {alpha**2:.4f}] / コントラスト比 {Io/Ib:.3f}")

出力は次のとおりです。

alpha=1.0: 大円(phi=0.30) 1.6317 [理論 1.6000] / 小円(phi=0.15) 1.2774 [理論 1.3000] / 背景 0.9615 [理論 1.0000] / コントラスト比 1.697
alpha=0.5: 大円(phi=0.30) 0.6066 [理論 0.5500] / 小円(phi=0.15) 0.3986 [理論 0.4000] / 背景 0.2320 [理論 0.2500] / コントラスト比 2.614
alpha=0.3: 大円(phi=0.30) 0.3359 [理論 0.2700] / 小円(phi=0.15) 0.1865 [理論 0.1800] / 背景 0.0796 [理論 0.0900] / コントラスト比 4.220

上段左端は一様な灰色で、どこに物体があるかまったく分かりません（強度の表示レンジを $\pm1\%$ まで拡大しても真っ平らです）。ところが右の3枚では、$\phi=0.30$ の大円と $\phi=0.15$ の小円が明るさの違う2段階として現れています。下段中央の散布図は、その明るさが位相 $\phi$ に対して直線状に並ぶこと、すなわち $I\approx\alpha^2+2\alpha\phi$ という予測が成り立っていることを示しています。下段右の棒グラフでは、$\alpha$ を下げるほど実測のコントラスト比が 1.70 → 2.61 → 4.22 と上がり、理論の傾向（1.60 → 2.20 → 3.00）と同じ向きに動いています。実測が理論をやや上回るのは、$\phi=0.30$ が弱位相近似の想定より大きく、高次項が効いているためです。

見どころは3つあります。第一に、さっきまで完全に真っ平らだった像に、位相板1枚で明確なコントラストが生まれました。第二に、強度が $I \approx \alpha^2 + 2\alpha\phi$ という位相の1次式によく従っています。$\phi=0.30$ と $\phi=0.15$ の円がきちんと2段階の明るさに分かれ、しかも明るさの差が位相差に比例している。位相コントラスト法は定量的なのです。第三に、$\alpha$ を 1.0 → 0.5 → 0.3 と減らすとコントラスト比が 1.70 → 2.61 → 4.22 と単調に上がりました。理論の $1+2\phi/\alpha$（それぞれ 1.60、2.20、3.00）よりやや大きめですが、これは $\phi=0.3$ が「弱位相」としてはやや大きく、$\phi^2$ 以上の高次項が効いているためです。背景を暗くするほど感度が上がるという設計指針は、数値でもはっきり確認できました。

なお $\alpha=1$ での背景が理論値1.000ではなく0.9615になっているのは、物体の存在によって直流成分 $\langle e^{j\phi}\rangle$ の絶対値が1よりわずかに小さくなるためです。これも「位相物体が視野の一部を占めると背景の明るさが変わる」という実際の顕微鏡の挙動に対応します。

実験7：ネガティブ位相コントラスト・暗視野・ナイフエッジの比較

同じ位相物体に、残る3つのフィルタを掛けて比べます。

import numpy as np
import matplotlib.pyplot as plt

# ネガティブ位相コントラスト（DCに -j）
I_neg = np.abs(four_f(u_ph, zernike_plate(1.0, sign=-1)))**2
# 暗視野（DCを完全に遮断）
P_dark = np.ones((N, N), dtype=complex); P_dark[N//2, N//2] = 0
I_dark = np.abs(four_f(u_ph, P_dark))**2
# ナイフエッジ（半平面 fx<0 を遮断、DCは通す）
P_knife = (XF >= -1e-12).astype(complex)
I_knife = np.abs(four_f(u_ph, P_knife))**2

print(f"ネガティブ位相板: 物体部 {I_neg[obj_mask].mean():.4f} / 背景 {I_neg[bg_mask].mean():.4f}")
print(f"暗視野          : 物体部 {I_dark[obj_mask].mean():.5f} / 背景 {I_dark[bg_mask].mean():.5f}"
      f"  （理論 phi^2 = {0.30**2:.4f}）")
print(f"ナイフエッジ    : min {I_knife.min():.4f} / max {I_knife.max():.4f} / 背景 {I_knife[bg_mask].mean():.4f}")

fig, ax = plt.subplots(1, 4, figsize=(16, 4.2))
for a_, im, t in zip(ax, [np.abs(u_ph)**2, I_neg, I_dark, I_knife],
                     ["フィルタなし（完全に不可視）", "ネガティブ位相コントラスト",
                      "暗視野（直流遮断）", "ナイフエッジ（シュリーレン）"]):
    a_.imshow(im, cmap="gray", extent=[-4, 4, -4, 4]); a_.set_title(t)
plt.tight_layout(); plt.show()

結果は次のようになります。

ネガティブ位相板: 物体部 0.5179 / 背景 1.0301
暗視野          : 物体部 0.07928 / 背景 0.00030  （理論 phi^2 = 0.0900）
ナイフエッジ    : min 0.4950 / max 1.6382 / 背景 1.0000

上段の4枚は、入力もレンズ配置もまったく同じで、フーリエ面のマスクだけを取り替えたものです。ネガティブ位相板では円が背景より暗く沈み、暗視野では背景が真っ黒に落ちて円だけが浮かび、ナイフエッジでは円が左右で明暗の対を作る立体的な陰影になります。下段のプロファイルを見ると、ナイフエッジ像の実測（青）が理論の $1-\mathcal H\{\phi\}$（赤破線）にほぼ完全に重なっており、円の左端で1.64まで跳ね上がり右端で0.50まで沈む「レリーフ調」の正体がヒルベルト変換だと分かります。緑の参考線（$1+\phi$）が矩形のままなのと対比すると、ナイフエッジは位相の値ではなく位相の変わり目に応答していることがはっきりします。

ネガティブ位相板では明暗が反転し、物体部（0.518）が背景（1.030）のほぼ半分の暗さになりました。弱位相近似の予測は $I\approx1-2\phi = 0.40$ ですから、実測はやや明るめです。これは $\phi=0.30$ が「弱位相」としては大きく、$\phi^2/2$ 以上の項を無視できないためで、$\phi$ を 0.1 程度まで下げると予測にきちんと収束します。符号を変えるだけでコントラストが反転する、というゼルニケ法の性質そのものは明確に確認できます。暗視野では背景が 0.0003 とほぼ真っ黒になり、物体部だけが 0.079 で光ります。この値は理論値 $\phi^2 = 0.09$ に近く、「直流を捨てると残るのは回折光だけ、その強度は位相の2乗」という予測どおりです。絶対的な明るさは位相コントラストより2桁暗いのに、背景がゼロなのでコントラストは圧倒的、という暗視野の性格がよく出ています。

ナイフエッジの像は、位相が一定の内部では背景と同じ明るさなのに、円の縁の左右で明暗が対になって現れるという特徴的なパターンになります。まさにシュリーレン写真の見た目です。理論では $I\approx 1-\mathcal H\{\phi\}$ でした。これを直接確かめてみます。

import numpy as np

# ヒルベルト変換（x方向）: 周波数領域で -j*sgn(fx) を掛ける
H_phi = np.real(IFT(FT(phi.astype(complex), dx) * (-1j*np.sign(FX)), dx))
pred = 1 - H_phi
print(f"ナイフエッジ像 I : min {I_knife.min():.4f} / max {I_knife.max():.4f}")
print(f"予測 1 - H[phi]  : min {pred.min():.4f} / max {pred.max():.4f}")
print(f"両者の相関係数   : {np.corrcoef(I_knife.ravel(), pred.ravel())[0,1]:.4f}")
print(f"最大絶対差       : {np.abs(I_knife - pred).max():.4f}")

出力は ナイフエッジ像 I : min 0.4950 / max 1.6382、予測 1 - H[phi] : min 0.4254 / max 1.5746、相関係数 0.9959、最大絶対差 0.0807 でした。相関0.996 という一致は、ナイフエッジフィルタが位相のヒルベルト変換を強度として書き出しているという理論を、数値的に裏付けています。残差は $\phi$ の2次以上の項に由来し、位相を小さくすればさらに縮みます。

シュリーレン法が「密度勾配を見ている」と説明されるのも納得できます。ヒルベルト変換は微分そのものではありませんが、周波数領域で $-j\,\mathrm{sgn}(f_X)$ を掛ける（＝各成分を90°回す）という操作は、$-j2\pi f_X$ を掛ける微分と方向性を共有しています。だから急に位相が変わる場所ほど強く応答するのです。

実務上の注意点

シミュレーションと実験の橋渡しとして、押さえておきたい落とし穴を挙げます。

サンプリングとエイリアシング。 離散シミュレーションでは、フーリエ面で表現できる最大空間周波数が $f_{\max} = 1/(2dx)$、周波数の刻みが $\Delta f_X = 1/L$ に固定されます。物理座標に直せば、フーリエ面の視野は $\pm\lambda f/(2dx)$、画素サイズは $\lambda f/L$ です。上の実験4で組んだ4f系の設定（$\lambda f = 6.33\times10^{-8}$ m$^2$、$L=8$ mm）ではフーリエ面の1画素が 7.91 µm でした。位相板やピンホールをこれより小さく作ろうとしても表現できません。逆に実験では、$\lambda f/L$ より細かい構造をマスクに刻む必要があるかどうかで、レンズ焦点距離の選択が決まります。

レンズの有限開口は常にローパス。 理論の導出では $u$ の積分を $\pm\infty$ まで取りましたが、実際のレンズには縁があります。半径 $D/2$ の開口は $f_c = (D/2)/(\lambda f) = 1/(2\lambda F_\#)$（$F_\# = f/D$ はFナンバー）という遮断周波数を持ちます。$\lambda=633$ nm、$F_\#=4$ なら $f_c \approx 197$ 1/mm。どんなに設計を頑張っても、これより細かい構造は結像できません。この限界がアッベの回折限界で、露光装置が短波長化（i線 → KrF → ArF → EUV）と高NA化を進めてきた理由そのものです。

コヒーレントとインコヒーレント。 本記事はレーザーのようなコヒーレント照明を仮定しました。このとき系は複素振幅について線形で、伝達関数は $H(f_X,f_Y)=P(\lambda f f_X,\lambda f f_Y)$ でした。一方、ランプ照明のようなインコヒーレント照明では系は強度について線形になり、伝達関数は $H$ の自己相関を規格化したOTF（光学伝達関数）になります。遮断周波数はコヒーレント系の2倍に伸びますが、その代わり中間周波数のコントラストは低下します。「コヒーレントのほうが解像度が高い」とは単純に言えないのは、この違いのためです。

収差はフーリエ面の位相誤差。 理想レンズの位相は $-k\rho^2/2f$ でしたが、実際のレンズには4次以上の項が残ります。これは瞳関数を $P(x,y)=|P|e^{jW(x,y)}$ と書いたときの波面収差 $W$ にほかなりません。球面収差なら $W\propto\rho^4$、デフォーカスなら $W\propto\rho^2$。収差とは、フーリエ面に置いた「望まぬ位相マスク」である、というのがフーリエ光学の見方です。詳しくはレンズの収差の記事も参照してください。

位相板の有限サイズ。 実験6では位相板を1画素にしましたが、実際の位相差顕微鏡では位相リングは有限の幅を持ちます。すると位相物体の低周波成分まで位相回転を受け、大きな構造の内部が暗くなる「シェーディングオフ」や、輪郭の外側が過剰に光る「ハロー」が発生します。上のコードで P[N//2, N//2] を半径数画素の円マスクに変えれば、この劣化を再現できます。

以上の注意点を踏まえれば、シミュレーションと実験の食い違いのほとんどは説明が付きます。

まとめ

本記事では、フーリエ光学の中核である「レンズはフーリエ変換器である」という主張を、回折積分から証明し、4f系の空間フィルタリングまで一気通貫で辿りました。

• フレネル回折積分は、自由空間伝搬が二次位相チャープ $\exp[jk\rho^2/2z]$ との畳み込みであることを示す。フラウンホーファー近似で遠方場はフーリエ変換になるが、開口1 mmでも1 m以上離れる必要があり、しかも残留二次位相が残る。
• 薄肉レンズは二次位相板 $t_l = \exp[-jk(x^2+y^2)/2f]$ である。レンズメーカーの公式を近軸近似で導けば、この形が自然に出てくる。伝搬のチャープと符号が逆なのが本質。
• 入力を前側焦点面に置くと、位相因子が完全に相殺する。核心は恒等式 $(u-\xi)^2-u^2+(x-u)^2 = [u-(\xi+x)]^2-2x\xi$ であり、$\xi^2$ と $x^2$ が同時に消えることで純粋なフーリエ核 $-2x\xi$ だけが残る。数値実験でも $d=f$ のとき位相が $10^{-3}$ rad 以下で平坦になることを確認した。
• スケーリング則は $f_X = x/(\lambda f)$。焦点面の座標がそのまま空間周波数の目盛りになる。$\lambda f$ が光学系設計の換算係数。
• 4f系の出力は $U_3 = U_1 * h$（$h=\mathcal F^{-1}\{H\}$、$H(f_X,f_Y)=P(\lambda f f_X,\lambda f f_Y)$）。デジタルの「FFT → マスク → 逆FFT」を、レンズ2枚とマスク1枚が光速で実行する。
• フィルタの実例：ピンホール（ローパス、ビームクリーンアップ、遮断 $f_c=a/\lambda f$）、遮蔽円板（ハイパス、暗視野、縁/内部比65.7を確認）、ナイフエッジ（シュリーレン、$I\approx1-\mathcal H\{\phi\}$、相関0.996）、ゼルニケ位相板（位相コントラスト、$I\approx\alpha^2+2\alpha\phi$、背景を暗くするほど高感度）。
• 完全に不可視だった位相物体（強度 min=max=1.000000）が、フーリエ面に部品を1つ足すだけで定量的に可視化される。これがフーリエ光学の最も鮮やかな成果。

この視点を身につけると、光学系を「レンズと絞りの並び」ではなく「線形システムとその伝達関数」として設計できるようになります。収差は瞳の位相誤差、解像限界は伝達関数の遮断周波数、アポダイゼーションは窓関数の設計 — すべてが信号処理の語彙で語れるようになるのです。

次のステップとして、以下の記事も参考にしてください。

• 単スリットのフラウンホーファー回折 — 本記事が前提とした遠方場回折の基本形
• 回折格子の分光分解能 — 実験2で見た回折次数の応用
• ガウスビーム光学の基礎 — レーザービームの伝搬とビームウエストの扱い
• レンズの収差とその評価 — 瞳関数の位相誤差としての収差
• 畳み込み定理をわかりやすく解説 — 4f系の出力式を支える数学
• 高速フーリエ変換（FFT）のアルゴリズム — 本記事のシミュレーションを支える計算法

参考文献

• J. W. Goodman, Introduction to Fourier Optics, 4th ed., W. H. Freeman, 2017
• M. Born and E. Wolf, Principles of Optics, 7th ed., Cambridge University Press, 1999
• F. Zernike, “Phase contrast, a new method for the microscopic observation of transparent objects,” Physica, vol. 9, pp. 686–698, 1942
• G. S. Settles, Schlieren and Shadowgraph Techniques, Springer, 2001