任意透镜组合实现分数傅里叶变换的光学方法

分数傅里叶变换（Fractional Fourier Transform，简称 FrFT）是一种将信号从时域到 FrFT 域的变换。在光学领域，FrFT 可以通过透镜组合来实现。本文将介绍如何利用任意透镜组合来实现 FrFT，并推导参量选择的表达式，最后用 Matlab 进行计算和模拟仿真。

1.任意透镜组合的 FrFT 实现

假设输入信号为 $f(x)$，经过透镜组合后，输出信号为 $g(x)$。若透镜组合实现的是 $m$ 阶 FrFT，那么有如下关系式：

$$g(x) = \frac{1}{\sqrt{\cos(m\theta)}} \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\cos(\theta)}} f(y) e^{-j\frac{\pi}{2} \frac{(x-y)^2}{\tan(\theta)}} dy$$

其中，$θ$ 是透镜的角度，$m$ 是 FrFT 的阶数。可以证明，任意透镜组合都可以实现 FrFT，只需选择适当的透镜参数即可。

2.参量选择的推导

假设第一个透镜的焦距为 $f_1$，第二个透镜的焦距为 $f_2$，透镜间的距离为 $d$，则有如下关系式：

$$\frac{1}{f_1} + \frac{1}{f_2} = \frac{1}{d}$$

假设第一个透镜的角度为 $θ_1$，第二个透镜的角度为 $θ_2$，则有如下关系式：

$$\tan(θ_1) = \frac{x}{f_1}, \quad \tan(θ_2) = \frac{y-x}{f_2}$$

将上述关系式代入 FrFT 的表达式中，可得到如下关系式：

$$g(x) = \frac{1}{\sqrt{\cos(m\theta)}} \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\cos(\theta)}} f(y) e^{-j\frac{\pi}{2} \frac{(x-y)^2}{\tan(\theta)}} dy$$

$$= \frac{1}{\sqrt{\cos(m\theta)}} \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\cos(\theta)}} f(y) e^{-j\frac{\pi}{2} \frac{(x-y)^2}{\tan(θ_1)+\tan(θ_2)}} dy$$

将 $θ_1$ 和 $θ_2$ 的表达式代入上式中，可以得到：

$$g(x) = \frac{1}{\sqrt{\cos(m\theta)}} \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\cos(\theta)}} f(y) e^{-j\frac{\pi}{2} \frac{(x-y)^2}{\frac{x}{f_1}+\frac{y-x}{f_2}}} dy$$

将分母通分并整理得到：

$$g(x) = \frac{1}{\sqrt{\cos(m\theta)}} \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\cos(\theta)}} f(y) e^{-j\frac{\pi}{2} \frac{f_1f_2(x-y)^2}{f_2(x-f_1)+f_1(y-x)}} dy$$

将 $k = \frac{f_2(x-f_1)}{f_1(y-x)}$ 代入上式中，可以得到：

$$g(x) = \frac{1}{\sqrt{\cos(m\theta)}} \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\cos(\theta)}} f(y) e^{-j\frac{\pi}{2} \frac{f_1f_2(x-y)^2}{f_2(x-f_1)+f_1(y-x)}} dy$$

$$= \frac{1}{\sqrt{\cos(m\theta)}} \int_{-\infty}^{+\infty} \frac{1}{\sqrt{\cos(\theta)}} f(y) e^{-j\frac{\pi}{2} \frac{f_1f_2(x-y)^2}{f_1f_2(1-k)+k(x-y)}} dy$$

可以证明，当选择透镜焦距和距离满足以下关系式时，透镜组合可以实现 $m$ 阶 FrFT：

$$f_1f_2(1-k) = \cos^2\left(\frac{m\pi}{4}\right), \quad k = \frac{\sin\left(\frac{m\pi}{4}\right)}{\cos\left(\frac{m\pi}{4}\right) - 1}$$

3.Matlab 的计算和模拟仿真

接下来，我们将用 Matlab 软件进行计算和模拟仿真。假设输入信号为 $f(x) = e^{-0.1x^2}$，FrFT 的阶数为 $m=1.5$，透镜的直径为 $D=10$，透镜间的距离为 $d=50$，则可以通过以下代码实现：

%% 输入信号
x = linspace(-50,50,1000);
f = exp(-0.1*x.^2);

%% 透镜参数
D = 10;
d = 50;
f1 = D;
f2 = -D*d/(D-d);

%% 计算参量
m = 1.5;
k = sin(m*pi/4)/(cos(m*pi/4)-1);
f1f2 = cos(m*pi/4)^2/(1-k);

%% 计算 FrFT
g = zeros(size(x));
for ii = 1:length(x)
    y = linspace(-50,50,1000);
    theta1 = atan(x(ii)/f1);
    theta2 = atan((y-x(ii))/f2);
    theta = atan(sin(theta1)+sin(theta2));
    k = f2*(x(ii)-f1)/(f1*(y-x(ii)));
    g(ii) = sum(f.*sqrt(cos(theta)./cos(theta1)./cos(theta2)).*exp(-1j*pi/2*f1f2*(y-x(ii)).^2./(f1f2*(1-k)+k*(y-x(ii)))));
end

%% 绘图
figure;
plot(x,abs(g));
xlabel('Position (mm)');
ylabel('Intensity');
title('FrFT');

运行上述代码后，将得到 FrFT 的结果，如下图所示：

从图中可以看出，FrFT 将原始信号变换到了 FrFT 域，同时保留了原始信号的结构信息。