# rekosntrukce fáze GS algoritmem

function sfa10_1

close all

graphics_toolkit('qt')

n_iteraci=100;
L=5; N=128;
del=L/N;
delf=1/N/del;
[x,y]=meshgrid((-N/2:N/2-1)*del);
[fx,fy]=meshgrid((-N/2:N/2-1)*delf);


# skutecna amplituda

faze=50*x.*(x.^2+y.^2);
uu=exp(I*faze).*circ(x,y,0.5);   # koma + pupila




figure(1)
subplot(1,2,1)
imshow(real(uu),[])
subplot(1,2,2)
imshow(imag(uu),[])   # re,im cast skutecne amplitudy

# skutecne spektrum

GG=ft2(uu,del);


# mereni

I_u=abs(uu).^2;
I_G=abs(GG).^2;

figure(2)

subplot(1,2,1)
imshow(sqrt(I_u),[])
subplot(1,2,2)
imshow(sqrt(I_G),[])



# GS algoritmus


dim=size(x);

fi1=zeros(dim);


res=zeros(n_iteraci,1);

for i=1:n_iteraci
i
u1=sqrt(I_u).*exp(I*fi1);
G1=ft2(u1,del);
fi2=arg(G1);
G2=sqrt(I_G).*exp(I*fi2);
u2=ift2(G2,delf);
fi1=arg(u2);

res(i)=norm(abs(GG)-abs(G1),"fro");

endfor

# uprava celkove faze


m1=arg(uu(N/2-1:N/2+1,N/2-1:N/2+1))
m2=arg(u1(N/2-1:N/2+1,N/2-1:N/2+1))
del=mean(mean(m1-m2))
u1=u1*exp(I*del);

figure(3)
h=plot(res);
set(h,'linewidth',2);
h1=xlabel("iterace");
h2=ylabel("chyba predikce");
set(h1,"fontsize", 16)
set(h2,"fontsize", 16)
set(gca, "linewidth", 1.5, "fontsize", 15)


figure(4)
subplot(1,2,1)
imshow(real(u1),[])
subplot(1,2,2)
imshow(imag(u1),[])




# predikce

G=ft2(u1,del);

figure(5)
imshow(abs(G),[])


endfunction



function G=ft(g,del)
 G=fftshift(fft(ifftshift(g)))*del;
end


function g=ift(G,delf)
 g=fftshift(ifft(ifftshift(G)))*delf*length(G);
end


function G=ft2(g,del)
 G=fftshift(fft2(ifftshift(g)))*del^2;
end


function g=ift2(G,delf)
 g=fftshift(ifft2(ifftshift(G)))*delf^2*length(G)^2;
end

function g=circ(x,y,radius)
r2=x.^2+y.^2;
g=double(r2<radius^2);
endfunction

