close all;clc;clear
A=zeros(64);
A(15:20,20:40)=1;A(15:50,20:25)=1;
A(45:50,20:40)=1;A(30:34,20:35)=1;
% ppp=exp(rand(64)*pi*2*i);A=A.*ppp;
% Author's email: zjliu2001@163.com
figure;imshow(abs(A),[]);
Fa=fft2(fftshift(A));Fs=fftshift(Fa);
Am=abs(Fs);   % amplitude
Ph=angle(Fs); % phase
s=11; % 这表示边长吗？
cgh=zeros(64*s);
th=max(max(abs(Fs)));
qq=th/1.2; % 为什么除1.2，1.2表示什么意思？
Am(Am>qq)=qq; %表示什么意思？
q=1:s;
w=(s+1)/2; %这表示的是宽度吧？
for m=1:64;
    for n=1:64;
        h=round(Am(m,n)/qq*(w-1)-0.5); %表示什么意思？
        md=zeros(s);
        if h>0;
          td=ones(h*2+1,3); %表示什么意思？h*2+1,3是什么意思？
          Pm=round(Ph(m,n)/pi*3); %表示什么意思？为什么要除pi再乘3，乘3是啥意思
          kz=Pm+w; %表示什么意思？
          md(w-h:w+h,kz-1:kz+1)=td; %表示什么意思？
        end
        cgh((m-1)*s+q,(n-1)*s+q)=md; %表示什么意思？
    end
end
figure;imshow(cgh,[]);
Re=ifft2(cgh);  Re=fftshift(Re);
figure;imshow(abs(Re),[]);