Gerchberg–Saxton algorithm

Hello! I have the following code:

for k=1:1:20;
G_pr=absP.*exp(i.*theta); 
g_pr=ifft2(ifftshift(G_pr)); 
absPhase=abs(angle(g_pr)); 
maxPh=max(max(absPhase)); 
minPh=min(min(absPhase)); 
g_pr(absPhase>=(minPh+0.2*(maxPh-minPh)))=0; 
g_pr=real(g_pr); 
gg=255*g_pr/(max(max(g_pr))); 
figure(1),imshow(uint8(gg)); title(num2str(k)) 
G=fftshift(fft2(g_pr)); 
G=G./abs(G); 
theta=angle(G); 
end

The first theta is the phase of my model image (angle(model)) However, this code diverges instead of converge, Does someone knows why?

Thank you

[Image Analyst]

You haven't given us enough code to even run your snippet that you posted here. Can't you step through it with the debugger to find out why?

[maria]

Gerchberg–Saxton algorithm is made in four steps:

• step 1: Fourier transform of the model
• step 2: Replace the modulus of the model-FT by the signal acquired
• step 3: Fourier-inverse step 2
• step 4: Create a new model with modulus of the FT-signal and the phase is step 3. Am I right?Since yesterday I have simplified a lot my code to find the mistake. So now it looks like:

model=zeros(150);
model(5:20, 20:121)=1;
model(24:140, 40:50)=1;
model(24:140, 85:97)=1;
figure(1);imagesc(model); colormap(gray)
DS=model; 
model=2*pi.*rand(150)-pi; 

pi_signal=double(imread('pi.bmp'));
pi_signal=rgb2gray(pi_signal);
pi_signal=imresize(pi_signal, [150 150]);

j=fftshift(fft2(pi_signal));
j=abs(j)/(max(max(abs(j))));
j=imadjust(abs(j), [0; 0.01], [0; 1]);

% 
for k=1:1:150; 
    step1=fftshift(fft2(model));
    phase=angle(step1); 
    Gu=medfilt2(j).*exp(i*phase);
    gx=fft2(fftshift(Gu));
    model=pi_signal.*exp(i*angle(gx));
    figure(7),imagesc(abs(gx)); colormap(gray), title(num2str(k));
end

where "model" is the model in step 1, "pi_signal" is the Fourier Transform of my signal and "j" is my signal.

Answer 1

[PNZ BDCB]

I'm not sure where in your code is the error. The following code works perfectly for me (adopted from wikipedia):

A = fftshift(ifft2(fftshift(Target)));
for i=1:25
  B = abs(Source) .* exp(1i*angle(A));
  C = fftshift(fft2(fftshift(B)));
  D = abs(Target) .* exp(1i*angle(C));
  A = fftshift(ifft2(fftshift(D)));

    imagesc(abs(C)) %Present current pattern
    title(sprintf('%d',i));
    pause(0.5)
end

Before running the code, make sure 'Source' contains your input beam, for example:

Source = exp(-1/2*(xx0.^2+yy0.^2)/sigma^2);

And 'Target' contains your requested pattern.

The phase mask can be presented at the end of the for loop:

imagesc(angle(A))

[Isaac Oguntoye]

Thanks for the response. Did you consider using your code for a simple image like a dot? Thanks.