example-of-gerchbergsaxton-algorithm.py

#We focused on implementing a frame-theoretic Gerchberg-Saxton algorith, where the input of the algorithm is the magnitude of frame coefficients. In this example, we use a full spark frame and a randomly generated vector x to construct frame coefficients. 
#This is a simple example in R^3. 

x = np.array([0.2994, 0.3455, 1.1124, 0.3081, 1.159])
n=len(x)
phi = np.random.random(n)*2*np.pi
MX=np.abs(x)
MF=np.abs(fft(x))
y_o= MX*np.exp(1j*phi)
y_n=y_o
for i in range(20): 
    y_next=find_next(y_n, MX, MF)
    print(np.linalg.norm(y_next-y_n),np.linalg.norm(np.abs(fft(y_n))-MF))
    y_n = y_next
print(y_next,x)

#With output:

1.1709102853402724 2.8513418802733734
0.12922534683061873 0.5605597368432665
0.09736084446208405 0.4135659439365416
0.07890094837842093 0.285629698321044
0.045874928781473555 0.17284790855657625
0.024755710674660564 0.11405220909893807
0.01627721680169051 0.08770176448852708
0.012302498109894823 0.072370009453755
0.009807994665868373 0.06171591427122309
0.0080605787572477 0.05379479984558585
0.006776266818753968 0.04767060220838524
0.005799722588593616 0.04279481591453913
0.005036815800519569 0.03882151031396196
0.004427433481753667 0.03552165483083709
0.0039315801213237545 0.03273761285643242
0.0035217220355200905 0.03035734869503288
0.0031783505137311873 0.028299058654739646
0.002887297888205702 0.0265016231093695
0.002638053401924713 0.02491846084389613
0.002422672500878094 0.023513450934189856
[-0.29155495-0.06808869j -0.335974  -0.0805712j  -1.09621826-0.18904839j
 -0.30670602+0.02927504j -1.14239313-0.19549665j] [0.2994 0.3455 1.1124 0.3081 1.159 ] 
 
 #Note: in real phase retrieval, we recover vector up to a plus or minus sign.