Gerchberg–Saxton algorithm

23 次查看（过去 30 天）
maria
maria 2015-9-27
评论： Isaac Oguntoye 2018-5-31
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
  2 个评论
Image Analyst
Image Analyst 2015-9-27
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
maria 2015-9-28
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

PNZ BDCB
PNZ BDCB 2017-10-25
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))
  1 个评论
Isaac Oguntoye
Isaac Oguntoye 2018-5-31
Thanks for the response. Did you consider using your code for a simple image like a dot? Thanks.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Images 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by