使用matlab,numpy和tensorflow2实现快速傅里叶变换FFT 交叉验证

以光学中圆孔的弗朗禾费衍射为例

matlab

N = 256;  % grid number
D = 1;    % diameter of aperture pupil [m]
L = 5.0*D/1.22;  % Airy disk occupies 5 pixels/grids [m]
x = linspace(-L/2.0, L/2.0, N);
y = linspace(-L/2.0, L/2.0, N);
[xm, ym] = meshgrid(x, y);
r = sqrt(xm.^2 + ym.^2);  % 注意加.
AP = double((r <= D/2.0));  
psf = abs(fftshift(fft2(AP))).^2;
subplot(121);imagesc(AP);title('Pupil');
axis square; colorbar;
subplot(122);imagesc(psf);title('Point Spread Function');
axis square; colorbar;

python

先生成孔径函数

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.io import savemat
print(tf.__version__)  # 2.1.0

# 生成Pupil函数
N = 256
D = 1
L = 10.0*D/1.22
x = np.linspace(-L/2.0, L/2.0, N)  
y = np.linspace(-L/2.0, L/2.0, N)
[xm, ym] = np.meshgrid(x, y)
r = (xm**2.0 + ym**2.0)**0.5
AP = np.double((r <= D/2.0))  
print(AP.dtype)  # float64

hf = plt.figure()
ax = Axes3D(hf)
ax.plot_surface(xm, ym, AP, cmap='jet')
ax.view_init(90, 0)

savemat('pupil.mat', mdict={'AP': AP})  # 保存数据到matlab中验证

输出
2.1.0
float64

numpy ===> np.fft.fft2

psf = np.fft.fftshift(np.fft.fft2(AP))

psf = (np.abs(psf))**2
print(psf.dtype)
print(np.max(psf))

hf2 = plt.figure()
ax2 = Axes3D(hf2)
ax2.plot_surface(xm, ym, psf, cmap='jet')
ax2.view_init(90, 0)
savemat('psf.mat', mdict={'psf': psf})

输出
float64
571536.0 与matlab的输出结果是对的上的

tensorflow ====>tf.fft2d(v1.) 或tf.signal.fft2d(v2.)

AP_t = tf.complex(AP,tf.zeros_like(AP))  # 必须转化为复数
print(AP_t.dtype)
psf_t = tf.signal.fftshift(tf.signal.fft2d(AP_t))  
psf_t = tf.square(tf.abs(psf_t))
print(psf_t.dtype)
psf_t = psf_t.numpy()
print(psf_t.dtype)
print(np.max(psf_t))

hf3 = plt.figure()
ax3 = Axes3D(hf3)
ax3.plot_surface(xm, ym, psf_t, cmap='jet')
ax3.view_init(90, 0)
savemat('psf_t.mat', mdict={'psf_t': psf_t})

输出
<dtype: ‘complex128’>
<dtype: ‘float64’> tensorflow输出的数据类型
float64 numpy输出的数据类型
571536.0 与matlab的输出结果是对的上的

python的绘图效果不太好，但经与matlab的数值结果对比验证，两者是一致。matlab的绘图比较好看是进行了插值运算