数字全息学习记录

文件：my_DDBFT.m

1、S-IFFT运算

在程序中有这么一段

%% 建立无干扰全息图
Uf = ifft2(Uf,N,N);	% 这里为什么要做一次傅里叶逆变换，可以参考S-FFT傅里叶正变换的时候，是怎么来的，具体推导参考本文件夹下面的README.md文件
Uf = fftshift(Uf);
n = 1:N;
x = -L/2 + L/N*(n-1);
y = x;
[xx,yy] = meshgrid(x,y);
Frensel = exp(1j*k/2/zi *(xx.^2+yy.^2));
U0 = Uf.*conj(Uf);	%  U0,S-IFFT

刚开始1时候，对于这一段非常的不理解，直到后来，仔细观察S-FFT衍射运算的公式才明白过来（原来严格来说，代码错了！），具体推导如下

1、S-FFT正向衍射公式

$$ \begin{align} U_Z(x,y) &= \frac{exp(jkd)} { (j\lambda d)} \\ \\ exp\big[ \frac{jk}{2d}(x^2+y^2)\big]\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\Big{ U_0(x_0,y_0)\\ exp\Big[\frac{jk}{2d}(x_0^2+y_0^2)\Big] \Big}exp\Big[ -j\frac{2\pi}{\lambda d} \big( x_0x + y_0y \big) \Big]\\ &=\frac{exp(jkd)} { (j\lambda d)} \\ \\ exp\big[ \frac{jk}{2d}(x^2+y^2)\big]\times FFT\Big{ U_0(x_0,y_0) exp\Big[\frac{jk}{2d}(x_0^2+y_0^2)\Big] \Big} \end{align} \tag{1.1.1} $$

变量代换，令 $$ phase = \frac{exp(jkd)} { (j\lambda d)} \\ \\ exp\big[ \frac{jk}{2d}(x^2+y^2)\big] \tag{1.1.2} $$

将公式(1.1.2)回带到公式(1.1.1)中去，可得 $$ U_Z(x,y) = phase \times FFT\Big{ U_0(x_0,y_0) exp\Big[\frac{jk}{2d}(x_0^2+y_0^2) \Big}\tag{1.1.3} $$ 下面开始，逆运算的推导

2、S-IFFT逆向衍射公式

首先，通过公式(1.1.3)可得 $$ U_Z(x,y)\\ \\ phase^* =phase\\ \\ phase^* \\ \\ FFT\Big{ U_0(x_0,y_0) exp\Big[\frac{jk}{2d}(x_0^2+y_0^2) \Big}\tag{1.2.1} $$ 那么将公式(1.2.1)两边做一次傅里叶逆变换可得 $$ \begin{align} FFT^{-1}\Big{ U_Z(x,y)\\ \\ phase^* \Big} &=|phase|^2\\ \\ FFT^{-1}\Big{ FFT\Big{ U_0(x_0,y_0) exp\Big[\frac{jk}{2d}(x_0^2+y_0^2) \Big} \Big}\\ FFT^{-1}\Big{ U_Z(x,y) phase^\Big} &=|phase|^2\\ \\ U_0(x_0,y_0) exp\Big[\frac{jk}{2d}(x_0^2+y_0^2)\Big] \end{align} \tag{1.2.2} $$ 变量带换，令 $$ Frensel = exp\Big[\frac{jk}{2d}(x_0^2+y_0^2)\Big]\tag{1.2.3} $$ 将公式(1.2.3)回带到公式(1.2.2)中，得 $$ FFT^{-1}\Big{ U_Z(x,y)\\ phase^ \Big} = U_0(x_0,y_0)\\ Frensel\\ \\ |phase|^2 \tag{1.2.4} $$ 将公式（1.2.4）两边公式乘以Frensel的共轭 $$ Frensel^* \times FFT^{-1}\Big{ U_Z(x,y)\\ phase^* \Big} = U_0(x_0,y_0)\\ Frensel\times Frensel^* |phase|^2 = U_0(x_0,y_0) \times |Frensel|^2 |phase|^2 \tag{1.2.5} $$ 将公式（1.2.5）整理的好看一点，就是下面公式（1.2.6）的样子 $$ U_0(x_0,y_0)\\ \\ |Frensel|^2\\ \\ |phase|^2 = Frensel^* \times FFT^{-1}\Big{ U_Z(x,y)\\ phase^* \Big} \tag{1.2.6} $$ 所以，最终的matlab代码，绝对不应该是他给的那样的，他没有将前面的相位因子减掉，正确的代码应该是这样的

%% 建立无干扰全息图
Uf = Uf.*conj(phase);	% phase 在上面有
Uf = ifft2(Uf,N,N);	% 这里为什么要做一次傅里叶逆变换，可以参考S-FFT傅里叶正变换的时候，是怎么来的，具体推导参考本文件夹下面的README.md文件
Uf = fftshift(Uf);
n = 1:N;
x = -L/2 + L/N*(n-1);
y = x;
[xx,yy] = meshgrid(x,y);
Frensel = exp(1j*k/2/zi *(xx.^2+yy.^2));
U0 = Uf.*conj(Uf);	%  U0,S-IFFT

这里注意一下，$Frensel$这个变量涉及到的坐标系，都是$x_0,y_0$，也就是原始物平面

3、抽样定理

满足振幅抽样定理的条件 $$ \Delta x_0^2 \ge \frac{\lambda d}{N} \tag{1.3.1} $$ 满足相位抽样定理的条件 $$ \Delta x_0^2 \le\frac{\lambda d}{N}\tag{1.3.2} $$ 同时满足振幅抽样定理和相位抽样定理的条件 $$ \Delta x_0 = \sqrt{\frac{\lambda d}{N}} \or\L_0 = L =\sqrt{\lambda d N} \tag{1.3.3} $$

4、改变光的传播方向

[xc,yc] = ginput(1);	% 括号中的1代表选择1个点
Lx = Li/N * (N-xc);	% 计算中心点到最右边的物理宽度
Ly = Li/N * (N-yc);	% 计算中心点到最下面的物理宽度

r = sqrt(Lx*Lx + z0*z0);	% 这一步是在计算位矢吗？课本上没有提到这一步啊
Qx = Lx / r;
Qy = Ly / r;
%	……
%	……
%	……
%	……
intvar = intvar.*exp(1j*k*(Qx.*x + Qy.*y));		% 变换传播方向让物光沿着光轴传播

在光的复振幅表达式中，波矢$\overrightarrow{k}$代表的是光的传播方向

这里有几点疑问：

1. 物光以前的传播反向难道不是衍射光轴的传播方向吗？物光以前的传播方向是沿着什么传播？
2. 乘上去之后，会怎样更改光的传播方向

最后，在程序里面测试了一下，发现，如果，在中间不成上那个改变光的传播方向的相位因子的话，最终再现的图像，会差一个fftshift

4、最终成图

最终呈现的图，发现会出现重影的现象，初步判断这是全息图生成的问题