伪代码:
Let:
FT – forward Fourier transform
IFT – inverse Fourier transform
i – the imaginary unit, √−1 (square root of −1)
exp – exponential function (exp(x)=ex)
Target and Source be the Target and Source Amplitude planes respectively
A, B, C & D be complex planes with the same dimension as Target and Source
Amplitude – Amplitude-extracting function:
e.g. for complex z = x + iy, amplitude(z) = sqrt(x·x + y·y)
for real x, amplitude(x) = |x|
Phase – Phase extracting function:
e.g. Phase(z) = arctan(y/x)
end Let
Gerchberg–Saxton Algorithm(Source, Target, Retrieved_Phase)
A = IFT(Target)
while error criterion is not satisfied
B = Amplitude(Source) * exp(i*Phase(A))
C = FT(B)
D = Amplitude(Target) * exp(i*Phase(C))
A = IFT(D)
end while
Retrieved_Phase = Phase(A)
end Gerchberg–Saxton Algorithm
python代码
import math
import numpy as np
def Ger_Sax_algo(img, max_iter):
h, w = img.shape
pm_s = np.random.rand(h, w)
pm_f = np.ones((h, w))
am_s = np.sqrt(img)
am_f = np.ones((h, w))
signal_s = am_s*np.exp(pm_s * 1j)
for iter in range(max_iter):
signal_f = np.fft.fft2(signal_s)
pm_f = np.angle(signal_f)
signal_f = am_f*np.exp(pm_f * 1j)
signal_s = np.fft.ifft2(signal_f)
pm_s = np.angle(signal_s)
signal_s = am_s*np.exp(pm_s * 1j)
pm =pm_f
return pm