clear all
close all
ph=peaks(256).*3; %预设原始相位
mu=0.3; %噪声系数
noise=rands(256,256).*mu.*pi; %噪声
ph=ph+noise; %叠加有噪声的预设原始相位
figure,imshow(ph,[]); %显示原始含噪声连续相位
colormap(jet);
U=exp(j*ph);
wrapped_phs =angle(U); %对应的包裹相位
figure,imshow(wrapped_phs,[]); %显示包裹相位
colormap(jet);
phi=wrapped_phs; %将包裹相位wrapped_phs赋值给a
[m,n]=size(phi); %计算二维包裹相位的大小（行、列数）
M=m-1;
N=n-1;
coosai=angle(U); %cosi为缠绕相位矩阵
X=zeros(M,N+1);
for P=1:N+1
for i=1:M
X(i,P)=mod(coosai(i+1,P)-coosai(i,P)+pi,2*pi)-pi;
end
end
XX=zeros(M,N+2);
XX(:,1)=X(:,2);
XX(:,2:N+2)=X;
derta_X=zeros(M+2,N+2);%derta_X 为缠绕相位X方向上的差分矩阵
derta_X(1,:)=-XX(1,:);
derta_X(2:M+1,:)=XX;
derta_X(M+2,:)=-XX(M,:);
for jj=1:N
for ii=1:M+1
Y(ii,jj)=mod(coosai(ii,jj+1)-coosai(ii,jj)+pi,2*pi)-pi;
end
end
YY=zeros(M+1,N+2);
YY(:,1)=-Y(:,1);
YY(:,2:N+1)=Y;
YY(:,N+2)=-Y(:,N);
derta_Y=zeros(M+2,N+2);%derta_Y 为缠绕相位Y方向上的差分矩阵
derta_Y(1,:)=YY(2,:);
derta_Y(2:M+2,:)=YY;
%%%%%%%%%%%%%%以上是求出x和y方向的差分矩阵%%%%%%%%%%%%%%%%%%%%%%
Lou=zeros(M+1,N+1);
for iii=1:M+1
for jjj=1:N+1
Lou(iii,jjj)=derta_X(iii+1,jjj+1)-derta_X(iii,jjj+1)+derta_Y(iii+1,jjj+1)-derta_Y(iii+1,jjj);
end
end
f=size(2*M,2*N);
f(1:(M+1),1:(N+1))=Lou;
f(1:(M+1),(N+2):2*N)=Lou(1:(M+1),N:-1:2);
f((M+2):2*M,1:(N+1))=Lou(M:-1:2,1:(N+1));
f((M+2):2*M,(N+2):2*N)=Lou(M:-1:2,N:-1:2);
P=fft2(f);
%%%%%%%%%以上是求出P(m,n)，即lou的二维离散傅里叶变换%%%%%%%%%%
for p=1:2*M
for q=1:2*N
k1=2*cos((p-1)*pi/(M));
k2=2*cos((q-1)*pi/(N));
PH(p,q)=P(p,q)/(k1+k2-4);
end
end
PH(1,1)=0;
phi1=ifft2(PH);
phi=phi1(1:M+1,1:N+1);
figure,imshow(real(phi),[]) %显示解包裹相位
colormap(jet);
%为检验解包裹效果，对解包裹相位作再包裹运算
% U=exp(j.*real(phi3)); %用解包裹相位构建光场
% rewrap_ph=atan(imag(U)./real(U)); %作再包裹运算
% figure,imshow(rewrap_ph,[]) %显示再包裹结果
% colormap(jet);
RMSE2=sqrt((sum(sum((real(phi)-ph).^2)))/(m*n));
RMSE2
