人類能看懂的衍射光學公式（內含基爾霍夫公式，瑞利--索末菲公式，角譜衍射公式，菲涅爾積分公式，以及菲涅爾衍射的S-FFT,T-FFT,D-FFT演算法））

考慮一些有基礎的小伙伴的需求，開頭先給結論，感興趣的可以慢慢看后面的推導程序，

結論：

基爾霍夫，瑞利–索末菲公式： U ( x , y , d ) = 1 j λ ∫ ? ∞ ∞ ∫ ? ∞ ∞ U ( x 0 , y 0 , 0 ) e x p ( j k r ) r K ( θ ) d x 0 d y 0 U(x,y,d) = \frac{1}{j\lambda} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U(x_0,y_0,0) \frac{exp(jkr)}{r} K(\theta)dx_0dy_0 U(x,y,d)=jλ1?∫?∞∞?∫?∞∞?U(x0?,y0?,0)rexp(jkr)?K(θ)dx0?dy0?‘

1. K ( θ ) = c o s ( θ ) + 1 2 K(\theta) = \frac{cos(\theta)+1}{2} K(θ)=2cos(θ)+1? 稱為基爾霍夫公式;

2. K ( θ ) = c o s ( θ ) K(\theta) = cos(\theta) K(θ)=cos(θ) 稱為第一種瑞利–索末菲公式;

3. K ( θ ) = 1 K(\theta) = 1 K(θ)=1 稱為第二種瑞利–索末菲公式;

角譜： U ( x , y , d ) = ? ? [ ? [ U ( x , y , 0 ) ] H ( f x , f y ) ] U(x,y,d) = \digamma^- [\digamma[U(x,y,0)]H(f_x,f_y)] U(x,y,d)=??[?[U(x,y,0)]H(fx?,fy?)]

菲涅爾衍射積分： U ( x , y , d ) = e x p ( j k d ) j λ d ∫ ? ∞ ∞ ∫ ? ∞ ∞ U ( x 0 , y 0 , 0 ) e x p [ j k 2 d [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U(x_0,y_0,0) exp[\frac{jk}{2d} [(x-x_0)^2+(y-y_0)^2] ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?∫?∞∞?∫?∞∞?U(x0?,y0?,0)exp[2djk?[(x?x0?)2+(y?y0?)2]]dx0?dy0?

夫瑯禾費衍射： U ( x , y , d ) = e x p ( j k d ) j λ d e x p [ j k 2 d ( x 2 + y 2 ) ] ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x 0 , y 0 , 0 ) e x p [ ? j 2 π λ d ( x x 0 + y y 0 ) ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} exp[\frac{jk}{2d}(x^2+y^2)] \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x_0,y_0,0) exp[\frac{-j2\pi}{\lambda d} (xx_0+yy_0) ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?exp[2djk?(x2+y2)]∫?∞∞?∫?∞∞?U0?(x0?,y0?,0)exp[λd?j2π?(xx0?+yy0?)]dx0?dy0?

S-FFT： U ( x , y , d ) = e x p ( j k d ) j k d e x p [ j k 2 d ( x 2 + y 2 ) ] ? [ U ( x 0 , y 0 , 0 ) e x p [ j k 2 d ( x 0 2 + y 0 2 ) ] ] U(x,y,d) = \frac{exp(jkd)}{jkd} exp[\frac{jk}{2d}(x^2+y^2)] \digamma[U(x_0,y_0,0)exp[\frac{jk}{2d}(x_0^2+y_0^2)] ] U(x,y,d)=jkdexp(jkd)?exp[2djk?(x2+y2)]?[U(x0?,y0?,0)exp[2djk?(x02?+y02?)]]

T-FFT： U ( x , y , d ) = e x p ( j k d ) j λ d ? ? [ ? [ U ( x 0 , y 0 , 0 ) ] ? [ e x p [ j k 2 d ( x 2 + y 2 ) ] ] ] U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \digamma^-[\digamma[U(x_0,y_0,0)] \digamma [exp[\frac{jk}{2d}(x^2+y^2)]] ] U(x,y,d)=jλdexp(jkd)???[?[U(x0?,y0?,0)]?[exp[2djk?(x2+y2)]]]

D-FFT： U ( x , y , d ) = ? ? [ ? [ U ( x 0 , y 0 , 0 ) ] e x p [ j k d ( 1 ? λ 2 2 ( f x 2 + f y 2 ) ) ] ] U(x,y,d) = \digamma^- [\digamma [U(x_0,y_0,0)]exp[jkd(1-\frac{\lambda^2}{2}(f_x^2+f_y^2))] ] U(x,y,d)=??[?[U(x0?,y0?,0)]exp[jkd(1?2λ2?(fx2?+fy2?))]]








接下來是上述衍射公式的推導程序：

首先，對一些基本概念和衍射的意義進行說明，

基本概念和衍射的意義

衍射（diffraction）是指波遇到障礙物時偏離原來直線傳播的物理現象，

波面的任何變形（通過相位物體）或者說波面（波前）上光場的復振幅分布受到任何空間調制，都將導致衍射現象的發生，而使通過障礙物以后的光場的復振幅重新分布，

衍射屏：導致衍射發生的障礙物

如果不涉及光傳播與變換程序中障礙物或光學元件結構尺寸接近于光波長的情況，對衍射問題的研究不鄰接衍射平面，
(這句話我TM想了好久沒想懂鄰接什么意思，查閱的時候一個個參考這句話參考得飛起但是都沒解釋，以我淺薄的理解就是在納米級別的情況下衍射作用可能不會發生，也就是本來應該導致衍射發生的障礙物不再認為是障礙物，)

標量衍射理論：將電場強度E視為標量，描述光傳播的物理程序和求解波動方程的方法，

傍軸近似：當光束發散較小且觀測區域的寬度遠小于光傳播距離，

實際應用中需要計算的往往是垂直于光學系統光軸或光束傳播方向的空間平面上的光波場，并且通常將z軸視為光學系統的光軸，也就是說，需要表述的平面是與光軸垂直的平面，
所以我們先以如下三維圖建立一個坐標系方便后續公式推導，設衍射屏與觀察屏的距離為d，矢徑 r = ( x ? x 0 ) 2 + ( y ? y 0 ) 2 + d 2 r = \sqrt{(x-x_0)^2+(y-y_0)^2+d^2} r=(x?x0?)2+(y?y0?)2+d2 ?, θ \theta θ 代表點 ( x 0 , y 0 , 0 ) (x_0,y_0,0) (x0?,y0?,0) 到點 ( x , y , d ) (x,y,d) (x,y,d) 的矢徑 r 與 ( x 0 , y 0 , 0 ) (x_0,y_0,0) (x0?,y0?,0)點法線 n n n的夾角，波矢 k = 2 π / λ k = 2\pi / \lambda k=2π/λ ， U ( x , y , 0 ) U(x,y,0) U(x,y,0)以及 U ( x , y , d ) U(x,y,d) U(x,y,d)分別為衍射屏及觀察屏上的光波的復振幅， G 0 ( f x , f y ) G_0(f_x,f_y) G0?(fx?,fy?), G d ( f x , f y ) G_d(f_x,f_y) Gd?(fx?,fy?)分別為其對應的頻譜函式，

基爾霍夫公式，瑞利–索末菲公式，以及角譜傳播公式分別是亥姆霍茲方程的三個解（光是一種電磁波，故其傳播特性可由經典電磁場理論的麥克斯韋方程推出，這里不具體展開，其實是筆者對于波動方程真的深惡痛絕），在標量衍射理論框架下，光傳播的物理程序可以嚴格地由基爾霍夫公式和瑞利–索末菲公式以及衍射的角譜的傳播公式表示，

如果知道空間中垂直于光傳播方向的一個平面上的光波場，便能計算該平面前后的空間光波場分布尤其是在觀察屏上的分布，
（xd后面沒有圖片了，全是密密麻麻的公式和文字）

基爾霍夫，瑞利–索末菲公式：

如果只是告訴大家基爾霍夫，瑞利–索末菲公式是亥姆霍茲方程的解而來可能難以理解其物理含義，但如果大家通過惠更斯–菲涅爾原理，以及場論中的格林定理電磁場的邊值條件等一步步推導就能逐漸窺見那些大佬的邏輯所在，這里就不展開來了，移步工程光學一書吧，

基爾霍夫，瑞利–索末菲公式兩者在數學上有統一的公式： U ( x , y , d ) = 1 j λ ∫ ? ∞ ∞ ∫ ? ∞ ∞ U ( x 0 , y 0 , 0 ) e x p ( j k r ) r K ( θ ) d x 0 d y 0 U(x,y,d) = \frac{1}{j\lambda} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U(x_0,y_0,0) \frac{exp(jkr)}{r} K(\theta)dx_0dy_0 U(x,y,d)=jλ1?∫?∞∞?∫?∞∞?U(x0?,y0?,0)rexp(jkr)?K(θ)dx0?dy0?

（其中 r = ( x ? x 0 ) 2 + ( y ? y 0 ) 2 + d 2 r = \sqrt{(x-x_0)^2+(y-y_0)^2+d^2} r=(x?x0?)2+(y?y0?)2+d2 ?, θ \theta θ 代表點 ( x 0 , y 0 , 0 ) (x_0,y_0,0) (x0?,y0?,0) 到點 ( x , y , d ) (x,y,d) (x,y,d) 的矢徑 r 與 ( x 0 , y 0 , 0 ) (x_0,y_0,0) (x0?,y0?,0)點法線 n n n的夾角，, K ( θ ) K(\theta) K(θ)為傾斜因子（由惠更斯–菲涅爾原理可知， K ( θ ) K(\theta) K(θ)表示子波的振幅隨面元法線與子波所在點到傳播點的夾角 θ \theta θ的變化），而不同的傾斜因子對應不同的公式：

1. K ( θ ) = c o s ( θ ) + 1 2 K(\theta) = \frac{cos(\theta)+1}{2} K(θ)=2cos(θ)+1? 稱為基爾霍夫公式;

2. K ( θ ) = c o s ( θ ) K(\theta) = cos(\theta) K(θ)=cos(θ) 稱為第一種瑞利–索末菲公式;

3. K ( θ ) = 1 K(\theta) = 1 K(θ)=1 稱為第二種瑞利–索末菲公式;

在實際問題中， θ \theta θ通常較小，三個公式的傾斜因子均接近1，

觀察屏上的光波的復振幅 U ( x , y , d ) = 1 j λ ∫ ? ∞ ∞ ∫ ? ∞ ∞ U ( x 0 , y 0 , 0 ) e x p ( j k r ) r d x 0 d y 0 U(x,y,d) = \frac{1}{j\lambda} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U(x_0,y_0,0) \frac{exp(jkr)}{r}dx_0dy_0 U(x,y,d)=jλ1?∫?∞∞?∫?∞∞?U(x0?,y0?,0)rexp(jkr)?dx0?dy0?

在傍軸情況（當光束發散較小且觀測區域的寬度遠小于光傳播距離）下，積分函式 e x p ( j k r ) r \frac{exp(jkr)}{r} rexp(jkr)?中的分母r可以由d來代替，因為這對于光波振幅無大的影響， r = ( x ? x 0 ) 2 + ( y ? y 0 ) 2 + d 2 r = \sqrt{(x-x_0)^2+(y-y_0)^2+d^2} r=(x?x0?)2+(y?y0?)2+d2 ?，這是由于傍軸情況下光束與z軸平行傳播時不會發散，也就是在x軸y軸上的分量不會有移動故 ( x ? x 0 ) 2 + ( y ? y 0 ) 2 = 0 (x-x_0)^2+(y-y_0)^2 = 0 (x?x0?)2+(y?y0?)2=0， r ≈ d r \approx d r≈d，
但是對于指數部分的r不能進行這樣簡單的處理，其原因是指數部分r的輕微變化會引起甚至大于 2 π 2\pi 2π的相位變化，

由 r = ( x ? x 0 ) 2 + ( y ? y 0 ) 2 + d 2 r = \sqrt{(x-x_0)^2+(y-y_0)^2+d^2} r=(x?x0?)2+(y?y0?)2+d2 ?(起初建立坐標系時就有的)

對 e x p ( j k r ) exp(jkr) exp(jkr)中的 r 進行泰勒展開， r ≈ d [ 1 + ( x ? x 0 ) 2 + ( y ? y 0 ) 2 2 d 2 ? [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] 2 8 d 4 ] r \approx d[1+ \frac{(x-x_0)^2+(y-y_0)^2}{2d^2} - \frac{[(x-x_0)^2+(y-y_0)^2]^2}{8d^4}] r≈d[1+2d2(x?x0?)2+(y?y0?)2??8d4[(x?x0?)2+(y?y0?)2]2?]……,第兩項之后的值過小，只保留泰勒展開的前兩項，故

e x p ( j k r ) exp(jkr) exp(jkr)中的 r ≈ d + ( x ? x 0 ) 2 + ( y ? y 0 ) 2 2 d r \approx d+ \frac{(x-x_0)^2+(y-y_0)^2}{2d} r≈d+2d(x?x0?)2+(y?y0?)2?
綜上， e x p ( j k r ) r ≈ e x p [ j k [ d + ( x ? x 0 ) 2 + ( y ? y 0 ) 2 2 d ] ] d \frac{exp(jkr)}{r} \approx \frac{exp[jk[d+ \frac{(x-x_0)^2+(y-y_0)^2}{2d}]]}{d} rexp(jkr)?≈dexp[jk[d+2d(x?x0?)2+(y?y0?)2?]]?

化簡后觀察屏上的光波的復振幅 U ( x , y , d ) = e x p ( j k d ) j λ d ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x 0 , y 0 , 0 ) e x p [ j k 2 d [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x_0,y_0,0) exp[\frac{jk}{2d} [(x-x_0)^2+(y-y_0)^2] ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?∫?∞∞?∫?∞∞?U0?(x0?,y0?,0)exp[2djk?[(x?x0?)2+(y?y0?)2]]dx0?dy0?，

上述那條公式也就是菲涅爾衍射積分，

而指數部分 r ≈ d + ( x ? x 0 ) 2 + ( y ? y 0 ) 2 2 d r \approx d+ \frac{(x-x_0)^2+(y-y_0)^2}{2d} r≈d+2d(x?x0?)2+(y?y0?)2?一般被稱為菲涅爾近似，滿足菲涅爾近似條件的衍射稱為菲涅爾衍射，

角譜衍射：

亥姆霍茲方程的第三個解便是角譜公式，在角譜理論中，光波沿z方向傳播的結果，在頻域內表現為將衍射屏上的光波場的頻譜 G 0 （ f x , f y ） G_0（f_x,f_y） G0?（fx?,fy?）乘以一個與z有關的相位延遲因子 e x p [ j 2 π λ z ( 1 ? ( λ f x ) 2 ? ( λ f y ) 2 ] exp[j\frac{2\pi}{\lambda}z\sqrt{(1-(\lambda f_x)^2-(\lambda f_y)^2}] exp[jλ2π?z(1?(λfx?)2?(λfy?)2 ?]，等效于通過一個半徑為1/ λ \lambda λ的理想低通濾波器，在線性系統理論中，該相位延遲因子即為衍射在頻域的傳遞函式 H ( f x , f y ) H(f_x,f_y) H(fx?,fy?)，表明衍射問題可以視為是光波場通過一個線性空間不變系統的變換程序，（ f x f_x fx? ， f y f_y fy?為頻域譜上對應的坐標系，也就是空間頻率，）

叨叨叨一堆話是不是不理解為什么呢？其實角譜理論就是將平面上的復振幅分布分解為多個不同傳播方向的平面波的線性疊加，類似于傅里葉變換把一個時域分量分解成多個頻率分量的總和那樣，

總結成公式就更為簡潔易懂了，就是

G d ( f x , f y ) = G 0 ( f x , f y ) e x p [ j 2 π λ z ( 1 ? ( λ f x ) 2 ? ( λ f y ) 2 ] G_d(f_x,f_y)=G_0(f_x,f_y)exp[j\frac{2\pi}{\lambda}z\sqrt{(1-(\lambda f_x)^2-(\lambda f_y)^2}] Gd?(fx?,fy?)=G0?(fx?,fy?)exp[jλ2π?z(1?(λfx?)2?(λfy?)2 ?]

（ G 0 ( f x , f y ) G_0(f_x,f_y) G0?(fx?,fy?), G d ( f x , f y ) G_d(f_x,f_y) Gd?(fx?,fy?)分別為 U ( x , y , 0 ) U(x,y,0) U(x,y,0)和 U ( x , y , d ) U(x,y,d) U(x,y,d)的頻譜函式，）

由上式可知若 1 ? ( λ f x ) 2 ? ( λ f y ) 2 < 0 1-(\lambda f_x)^2-(\lambda f_y)^2<0 1?(λfx?)2?(λfy?)2<0，則角譜分量將隨z的增大按指數規律急劇衰減，光波只存在于鄰近衍射屏的一個非常薄的區域，故只有 1 ? ( λ f x ) 2 ? ( λ f y ) 2 > 0 1-(\lambda f_x)^2-(\lambda f_y)^2>0 1?(λfx?)2?(λfy?)2>0，即 f x 2 + f y 2 < 1 λ 2 f_x^2+f_y^2<\frac{1}{\lambda^2} fx2?+fy2?<λ21?，角譜分量才會到達觀察屏，

（以 f x = f y f_x=f_y fx?=fy?， λ \lambda λ=632nm為例， f x = f y < 1.1188 e + 06 f_x=f_y<1.1188e+06 fx?=fy?<1.1188e+06 m）

簡寫下上面的公式，設衍射距離為d，定義角譜衍射的傳遞函式為

H ( f x , f y ) = e x p [ j 2 π λ d ( 1 ? ( λ f x ) 2 ? ( λ f y ) 2 H(f_x,f_y) = exp[j\frac{2\pi}{\lambda}d\sqrt{(1-(\lambda f_x)^2-(\lambda f_y)^2} H(fx?,fy?)=exp[jλ2π?d(1?(λfx?)2?(λfy?)2 ?

則角譜衍射公式可簡寫為

$U(x,y,d) = \digamma^- [\digamma[U(x,y,0)]H(f_x,f_y)] $

當傍軸近似（當光束發散較小且觀測區域的寬度遠小于光傳播距離）時，即 f x , f y f_x,f_y fx?,fy?較小時， ( 1 ? ( λ f x ) 2 ? ( λ f y ) 2 \sqrt{(1-(\lambda f_x)^2-(\lambda f_y)^2} (1?(λfx?)2?(λfy?)2 ?泰勒系數展開，同樣只保留前兩項

（ ( 1 ? ( λ f x ) 2 ? ( λ f y ) 2 = 1 ? 1 2 λ 2 ( f x 2 + f y 2 ) + 1 8 λ 4 ( f x 2 + f y 2 ) \sqrt{(1-(\lambda f_x)^2-(\lambda f_y)^2} = 1-\frac{1}{2}\lambda^2(f_x^2+f_y^2)+\frac{1}{8}\lambda^4(f_x^2+f_y^2) (1?(λfx?)2?(λfy?)2 ?=1?21?λ2(fx2?+fy2?)+81?λ4(fx2?+fy2?)……）,故

角譜衍射的傳遞函式 H ( f x , f y ) = e x p [ j k d ( 1 ? λ 2 2 ( f x 2 + f y 2 ) ) ] H(f_x,f_y) = exp[jkd(1-\frac{\lambda^2}{2}(f_x^2+f_y^2))] H(fx?,fy?)=exp[jkd(1?2λ2?(fx2?+fy2?))]

這只是對 H ( f x , f y ) H(f_x,f_y) H(fx?,fy?)的一個小近似，然后根據傅里葉變換性質，頻域分量的相乘等于空域分量的卷積，

G d ( f x , f y ) = G 0 ( f x , f y ) H ( f x , f y ) G_d(f_x,f_y)=G_0(f_x,f_y)H(f_x,f_y) Gd?(fx?,fy?)=G0?(fx?,fy?)H(fx?,fy?)便變換為

U ( x , y , d ) = U ( x 0 , y 0 , 0 ) ? ? ? [ H ( f x , f y ) ] U(x,y,d) = U(x_0,y_0,0) * \digamma^-[H(f_x,f_y)] U(x,y,d)=U(x0?,y0?,0)???[H(fx?,fy?)]

而 ? ? [ H ( f x , f y ) ] \digamma^-[H(f_x,f_y)] ??[H(fx?,fy?)]有決議解，即 ? ? [ H ( f x , f y ) ] = e x p ( j k d ) j λ d e x p [ j k 2 d ( x 2 + y 2 ) ] \digamma^-[H(f_x,f_y)] = \frac{exp(jkd)}{j\lambda d} exp[\frac{jk}{2d}(x^2+y^2)] ??[H(fx?,fy?)]=jλdexp(jkd)?exp[2djk?(x2+y2)]

角譜衍射卷積形式的公式即可改寫為菲涅爾衍射積分

U ( x , y , d ) = e x p ( j k d ) j λ d ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x , y , 0 ) e x p [ j k 2 d [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x,y,0) exp[\frac{jk}{2d} [(x-x_0)^2+(y-y_0)^2] ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?∫?∞∞?∫?∞∞?U0?(x,y,0)exp[2djk?[(x?x0?)2+(y?y0?)2]]dx0?dy0?

兜兜轉轉，最后還是你，盡管角譜衍射公式，基爾霍夫公式以及瑞利–索末菲公式有不同的形式，但它們的傍軸近似具有相同的形式，也就是菲涅爾積分形式，

夫瑯禾費衍射：

若菲涅爾衍射積分的傅里葉變換中，如果 d > > k ( x 0 2 + y 0 2 ) m a x 2 d>>\frac{k(x_0^2+y_0^2)_{max}}{2} d>>2k(x02?+y02?)max??，即衍射屏距離觀察屏無限遠，則 j k ( x 0 2 + y 0 2 ) 2 d ≈ 0 \frac{jk(x_0^2+y_0^2)}{2d} \approx 0 2djk(x02?+y02?)?≈0， e x p [ j k ( x 0 2 + y 0 2 ) 2 d ] ≈ 1 exp[\frac{jk(x_0^2+y_0^2)}{2d}] \approx 1 exp[2djk(x02?+y02?)?]≈1，故夫瑯禾費衍射場可簡單地變為 U 0 ( x 0 , y 0 , 0 ) U_0(x_0,y_0,0) U0?(x0?,y0?,0)的傅里葉變換，

U ( x , y , d ) = e x p ( j k d ) j λ d e x p [ j k 2 d ( x 2 + y 2 ) ] ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x , y , 0 ) e x p [ ? j 2 π d λ ( x x 0 + y y 0 ) ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} exp[\frac{jk}{2d}(x^2+y^2)] \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x,y,0) exp[\frac{-j2\pi}{d\lambda} (xx_0+yy_0) ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?exp[2djk?(x2+y2)]∫?∞∞?∫?∞∞?U0?(x,y,0)exp[dλ?j2π?(xx0?+yy0?)]dx0?dy0?

將上式對比于基爾霍夫–瑞利索末菲公式 U ( x , y , d ) = 1 j λ ∫ ? ∞ ∞ ∫ ? ∞ ∞ U ( x 0 , y 0 , 0 ) e x p ( j k r ) r d x 0 d y 0 U(x,y,d) = \frac{1}{j\lambda} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U(x_0,y_0,0) \frac{exp(jkr)}{r} dx_0dy_0 U(x,y,d)=jλ1?∫?∞∞?∫?∞∞?U(x0?,y0?,0)rexp(jkr)?dx0?dy0?， e x p [ j k ( x 0 2 + y 0 2 ) 2 d ] ≈ 1 exp[\frac{jk(x_0^2+y_0^2)}{2d}]\approx1 exp[2djk(x02?+y02?)?]≈1也就是對指數部分r進行近似，

r ≈ d + x 2 + y 2 2 d ? x x 0 + y y 0 d r \approx d+\frac{x^2+y^2}{2d}-\frac{xx_0+yy_0}{d} r≈d+2dx2+y2??dxx0?+yy0??

也就是在菲涅爾近似 r ≈ d + ( x ? x 0 ) 2 + ( y ? y 0 ) 2 2 d r \approx d+ \frac{(x-x_0)^2+(y-y_0)^2}{2d} r≈d+2d(x?x0?)2+(y?y0?)2? 的基礎上 x 0 2 + y 0 2 = 0 x_0^2 + y_0^2 = 0 x02?+y02?=0 ,可以理解為對比于在無限遠的觀察屏，衍射屏上的坐標分量已經小到可以忽略不計了，

上述近似即被稱為夫瑯禾費近似，滿足夫瑯禾費近似的衍射即被稱為夫瑯禾費衍射，

夫瑯禾費近似成立所要求的條件較為苛刻，但如果物平面光波場為一個會聚的球面波，則夫瑯禾費近似很容易滿足，而會聚球面波應用上可以很容易通過透鏡實作，這也是觀察夫瑯禾費衍射情況下在觀察屏前加透鏡的原因，








接下來是對菲涅爾衍射積分的三種常用演算法的推導：

S-FFT（single fast Fourier transform algorithm）:單次傅里葉變換演算法

菲涅爾衍射積分： U ( x , y , d ) = e x p ( j k d ) j λ d ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x , y , 0 ) e x p [ j k 2 d [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x,y,0) exp[\frac{jk}{2d} [(x-x_0)^2+(y-y_0)^2] ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?∫?∞∞?∫?∞∞?U0?(x,y,0)exp[2djk?[(x?x0?)2+(y?y0?)2]]dx0?dy0?

對上式二次相位因子 e x p [ j k 2 d [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] ] exp[\frac{jk}{2d} [(x-x_0)^2+(y-y_0)^2] ] exp[2djk?[(x?x0?)2+(y?y0?)2]]進行展開，波矢 k = 2 π / λ k = 2\pi / \lambda k=2π/λ ，并將與積分變數無關的項提到積分號前則得到

U ( x , y , d ) = e x p ( j k d ) j λ d e x p [ j k 2 d ( x 2 + y 2 ) ] ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x , y , 0 ) e x p [ j k 2 d [ ( x 0 2 + y 0 2 ] ] e x p [ ? j 2 π ( x λ d x 0 + y λ d y 0 ) ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} exp[\frac{jk}{2d}(x^2+y^2)] \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x,y,0) exp[\frac{jk}{2d} [(x_0^2+y_0^2] ] exp[-j2\pi(\frac{x}{\lambda d}x_0 + \frac{y}{\lambda d}y_0)] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?exp[2djk?(x2+y2)]∫?∞∞?∫?∞∞?U0?(x,y,0)exp[2djk?[(x02?+y02?]]exp[?j2π(λdx?x0?+λdy?y0?)]dx0?dy0?

對比于傅里葉變換的公式

F ( f x , f y ) = ∫ ? ∞ ∞ ∫ ? ∞ ∞ f ( x , y ) e x p [ ? j 2 π ( f x x + f y y ) ] d x d y F(f_x,f_y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f(x,y)exp[-j2\pi(f_x x+f_y y)]dxdy F(fx?,fy?)=∫?∞∞?∫?∞∞?f(x,y)exp[?j2π(fx?x+fy?y)]dxdy

便可以將二重積分項等同為 U 0 ( x , y , 0 ) e x p [ j k 2 d ( x 0 2 + y 0 2 ) ] U_0(x,y,0)exp[\frac{jk}{2d}(x_0^2+y_0^2)] U0?(x,y,0)exp[2djk?(x02?+y02?)]的傅里葉變換，設 G 1 ( f x , f y ) G_1(f_x,f_y) G1?(fx?,fy?)是 U 0 ( x , y , 0 ) e x p [ j k 2 d ( x 0 2 + y 0 2 ) ] U_0(x,y,0)exp[\frac{jk}{2d}(x_0^2+y_0^2)] U0?(x,y,0)exp[2djk?(x02?+y02?)]的傅里葉變換

且 f x = x λ d f_x = \frac{x}{\lambda d} fx?=λdx? ， f y = y λ d f_y = \frac{y}{\lambda d} fy?=λdy?

故菲涅爾衍射積分的S-FFT演算法為： U ( x , y , d ) = e x p ( j k d ) j λ d e x p [ j k 2 d ( x 2 + y 2 ) ] G 1 ( f x , f y ) U(x,y,d) = \frac{exp(jkd)}{j \lambda d} exp[\frac{jk}{2d}(x^2+y^2)] G_1(f_x,f_y) U(x,y,d)=jλdexp(jkd)?exp[2djk?(x2+y2)]G1?(fx?,fy?)

T-FFT(triple fast Fourier transform algorithm):三次傅里葉變換演算法

（兩次傅里葉變換一次逆傅里葉變換），這種演算法常被稱為菲涅爾衍射的卷積演算法，

菲涅爾衍射積分： U ( x , y , d ) = e x p ( j k d ) j λ d ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x , y , 0 ) e x p [ j k 2 d [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x,y,0) exp[\frac{jk}{2d} [(x-x_0)^2+(y-y_0)^2] ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?∫?∞∞?∫?∞∞?U0?(x,y,0)exp[2djk?[(x?x0?)2+(y?y0?)2]]dx0?dy0?

類比于卷積公式 y ( t ) = ∫ ? ∞ ∞ x ( p ) h ( t ? p ) d p = x ( t ) ? y(t) = \int_{-\infty}^{\infty} x(p)h(t-p)dp = x(t)\bigotimes y(t)=∫?∞∞?x(p)h(t?p)dp=x(t)?h(t)

定義脈沖回應

h ( x , y ) = e x p ( j k d ) j λ d e x p [ j k 2 d ( x 2 + y 2 ) ] h(x,y) = \frac{exp(jkd)}{j\lambda d} exp[\frac{jk}{2d}(x^2+y^2)] h(x,y)=jλdexp(jkd)?exp[2djk?(x2+y2)]

便可將菲涅爾衍射積分可轉化為 U ( x 0 , y 0 , 0 ) U(x_0,y_0,0) U(x0?,y0?,0)和 h ( x , y ) h(x,y) h(x,y)的卷積，

U ( x , y , d ) = U ( x 0 , y 0 , 0 ) ? h ( x , y ) U(x,y,d) = U(x_0,y_0,0) \bigotimes h(x,y) U(x,y,d)=U(x0?,y0?,0)?h(x,y)

把脈沖項的無關項提前，則式子變為

U ( x , y , d ) = e x p ( j k d ) j λ d U ( x 0 , y 0 , 0 ) ? e x p [ j k 2 d ( x 2 + y 2 ) ] U(x,y,d) = \frac{exp(jkd)}{j\lambda d} U(x_0,y_0,0) \bigotimes exp[\frac{jk}{2d}(x^2+y^2)] U(x,y,d)=jλdexp(jkd)?U(x0?,y0?,0)?exp[2djk?(x2+y2)]

但由于在空域完成卷積運算比較花機器時間，所以通常要轉換到頻域進行計算，根據頻譜特性，空域的卷積等于頻域的乘積，故需要對 U ( x 0 , y 0 , 0 ) U(x_0,y_0,0) U(x0?,y0?,0)和 e x p [ j k 2 d ( x 2 + y 2 ) ] exp[\frac{jk}{2d}(x^2+y^2)] exp[2djk?(x2+y2)]各做一次傅里葉變換，再對它們的乘積做一次逆傅里葉變換，

故菲涅爾衍射積分的T-FFT演算法為： U ( x , y , d ) = e x p ( j k d ) j λ d ? ? [ ? [ U ( x 0 , y 0 , 0 ) ] ? [ e x p [ j k 2 d ( x 2 + y 2 ) ] ] ] U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \digamma^-[\digamma[U(x_0,y_0,0)] \digamma [exp[\frac{jk}{2d}(x^2+y^2)]] ] U(x,y,d)=jλdexp(jkd)???[?[U(x0?,y0?,0)]?[exp[2djk?(x2+y2)]]]

與S-FFT演算法不同的是，觀察面的尺寸與衍射面的尺寸是一樣的，是因為在計算中在頻域進行相乘時，要求頻譜中相同頻率分量對應相乘，故最高頻率也就必須相同，抽樣數一樣對應的幾何尺寸也就相等了，

D-FFT（double fast Fourier transform algorithm）:兩次傅里葉變換演算法

與T-FFT演算法一樣，該演算法也常被稱為卷積演算法，

菲涅爾衍射積分： U ( x , y , d ) = e x p ( j k d ) j λ d ∫ ? ∞ ∞ ∫ ? ∞ ∞ U 0 ( x , y , 0 ) e x p [ j k 2 d [ ( x ? x 0 ) 2 + ( y ? y 0 ) 2 ] ] d x 0 d y 0 U(x,y,d) = \frac{exp(jkd)}{j\lambda d} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} U_0(x,y,0) exp[\frac{jk}{2d} [(x-x_0)^2+(y-y_0)^2] ] dx_0dy_0 U(x,y,d)=jλdexp(jkd)?∫?∞∞?∫?∞∞?U0?(x,y,0)exp[2djk?[(x?x0?)2+(y?y0?)2]]dx0?dy0?

把D-FFT放在T-FFT后面講是因為這兩者有相似之處，

同樣設脈沖回應 h ( x , y ) = e x p ( j k d ) j λ d e x p [ j k 2 d ( x 2 + y 2 ) ] h(x,y) = \frac{exp(jkd)}{j\lambda d} exp[\frac{jk}{2d}(x^2+y^2)] h(x,y)=jλdexp(jkd)?exp[2djk?(x2+y2)]

與T-FFT演算法類似，菲涅爾積分變換為 U ( x , y , d ) = U ( x 0 , y 0 , 0 ) ? h ( x , y ) U(x,y,d) = U(x_0,y_0,0) \bigotimes h(x,y) U(x,y,d)=U(x0?,y0?,0)?h(x,y)

從這里開始，D-FFT與T-FFT的演算法的推導就不一樣了，

設 H ( f x , f y ) H(f_x,f_y) H(fx?,fy?)為 h ( x , y ) h(x,y) h(x,y)對應的頻譜函式（ H ( f x , f y ) H(f_x,f_y) H(fx?,fy?)也被稱為菲涅爾傳遞函式，），由于 h ( x , y ) h(x,y) h(x,y) 的傅里葉變換有決議解，即 H ( f x , f y ) = e x p [ j k d ( 1 ? λ 2 2 ( f x 2 + f y 2 ) ) ] H(f_x,f_y) = exp[jkd(1-\frac{\lambda^2}{2}(f_x^2+f_y^2))] H(fx?,fy?)=exp[jkd(1?2λ2?(fx2?+fy2?))]

則菲涅爾衍射積分的D-FFT演算法： U ( x , y , d ) = ? ? [ ? [ U ( x 0 , y 0 , 0 ) ] e x p [ j k d ( 1 ? λ 2 2 ( f x 2 + f y 2 ) ) ] ] U(x,y,d) = \digamma^- [\digamma [U(x_0,y_0,0)] exp[jkd(1-\frac{\lambda^2}{2}(f_x^2+f_y^2))]] U(x,y,d)=??[?[U(x0?,y0?,0)]exp[jkd(1?2λ2?(fx2?+fy2?))]]

與T-FFT演算法不同，D-FFT，已經知道傳遞函式 H ( f x , f y ) H(f_x,f_y) H(fx?,fy?)的決議運算式，相當于把T-FFT對 h ( x , y ) h(x,y) h(x,y) 的FFT變換先行計算了出來，只需要將空間頻率 f x f_x fx? ， f y f_y fy?代入決議式即可得到傳遞函式的值而不用做FFT變換，但同時也要求必須通過實驗引數給出正確的 f x , f_x, fx?, f y f_y fy?

對比于角譜函式，可以發現角譜衍射的傳遞函式在傍軸近似條件下等于菲涅爾衍射的傳遞函式，也就是角譜衍射的公式在傍軸條件下等于菲涅爾積分的D-FFT演算法，

總結：

到此為止，已對常見的衍射公式進行了大致的介紹和簡單的推導，建議你再到開頭去看下結論，希望你能對這些經典的衍射公式有更深的理解，（然后還是要吐槽下國內的很多光學教材都編得不咋地，很多陳述句都是抄來抄去，參考起來也沒有相關解釋，很多淺顯的公式道理也講得繁瑣復雜(╯‵□′)╯︵┻━┻）

（附注：因為是手敲公式所以可能會有些許疏漏還望見諒QAQ，此外，也歡迎感興趣有疑問的小伙伴在評論區留言討論）

參考文獻：
《資訊光學數字實驗室（Matlab版）》——錢曉凡
《衍射計算及數字全息》—— 李俊昌
《計算機制全息圖》——金國藩，張浩，蘇萍，蔣強，曹良才
《工程光學》——郁道銀，談恒英
《資訊光學》——蘇顯渝