iterative=300;????????????%設迭代次數為300次吧??
imagename=‘cat.jpg‘;????%你想要提取相位的圖像名稱??
phaseimage=‘phase.png‘;??%要保存的相位圖像名稱??
??
??
%空域輸入圖像的幅度（是已知的，也就是清晰的圖像，它的灰度就是幅值）和相位圖像（待恢復）??
known_abs_spatial=imread（imagename）;????????????%作為輸入圖像的幅度，是已知的??
known_abs_spatial?=rgb2gray（known_abs_spatial）;?%注意要用單通道圖像做實驗，如果你讀取的是彩色圖像，那就吧這行取消注釋變成灰度圖像吧??
known_abs_spatial=im2double（known_abs_spatial）;?%將圖像灰度映射到0～1??
??
unknown_phase=known_abs_spatial;????????????????%Peppers圖像作為輸入圖像的相位，也即為待恢復的數據，??
????????????????????????????????????????????????%要求它和known_abs_spatial大小一致，所以這里直接賦值就好了??
unknown_phase=im2double（unknown_phase）;?????????%將圖像灰度映射到0～1??
unknown_phase2=unknown_phase*2*pi;??????????????%相位范圍映射到0-2*pi??
unknown_phase2（unknown_phase2>pi）=unknown_phase2（unknown_phase2>pi）-2*pi;%進一步映射至[-pi+pi]??
??
[widthlength]=size（known_abs_spatial）;?????????%獲取Lena圖像的大小??
input=known_abs_spatial.*exp（i*unknown_phase2）;?%最終輸入圖像:幅度*e^（i*相位角度），它是復數圖像??
known_abs_fourier=abs（fft2（input））;?????????????%先將input圖像進行傅立葉變換，然后取模，就是傅氏變換后的幅度??
%以下開始迭代求相位??
phase_estimate=pi*rand（widthlength）;???????????%這是生成了一副大小為（width*length）的圖像??
????????????????????????????????????????????????%它的像素值是在[0pi]范圍內隨機生成的。??
imshow（phase_estimate）??
%以下開始迭代??
for?p=1:iterative??
????signal_estimate_spatial=known_abs_spatial.*exp（i*phase_estimate）;???%Step?1??構造estimated?signal：還是幅度*e^（i*相位角度）變成復數形式??
????temp1=fft2（signal_estimate_spatial）;????????????????????????????????%傅立葉變換到頻域??
????temp_ang=angle（temp1）;??????????????????????????????????????????????%求相位弧度，它的范圍是[-pipi]??
????signal_estimate_fourier=known_abs_fourier.*exp（i*temp_ang）;?????????%Step?2??替換傅氏變換后的幅度，產生estimate?Fourier?transform??
????temp2=ifft2（signal_estimate_fourier）;???????????????????????????????%Step?3??對Step?2產生的estimate?Fourier?transform進行傅立葉反變換，又變換到空域了??
????phase_estimate=angle（temp2）;????????????????????????????????????????%Step?4:estimated?phase??
end??
%以上循環就是通過隨便預設一個相位圖像，在循環中不斷調整逼近真實的相位，直到滿足條件（也就是我們求的相位和真實相位非常接近的時候）??
%不過這里我們只需要設定一個比較大的循環就可以了，基本上都可以滿足條件了，這個激光原理就講過了。??
??
phase_estimate（phase_estimate<0）=phase_estimate（phase_estimate<0）+2*pi;?%把estimate_phase從[-pi+pi]，映射到[02pi]??
retrieved=phase_estimate/（2*pi）;%再映射到[01]??
??
figure?（1）??
imshow（retrieved）;title（‘相位圖像‘）%顯示我們提取到的相位圖像??
imwrite（retrievedphaseimage）?

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????目錄???????????0??2018-05-08?06:42??Gerchberg-Saxton-GS-algorithm-master\
?????文件????????2820??2018-05-08?06:42??Gerchberg-Saxton-GS-algorithm-master\GS_phase_retrieval.m
?????文件?????????228??2018-05-08?06:42??Gerchberg-Saxton-GS-algorithm-master\README.md
?????文件???????83724??2018-05-08?06:42??Gerchberg-Saxton-GS-algorithm-master\cat.jpg
?????文件??????174076??2018-05-08?06:42??Gerchberg-Saxton-GS-algorithm-master\phase.png