%% *******************苏德志论文仿真********************
clear;clc;close all;
M = 4;
for k = 1:M
delta(k) = (k-1)*pi/2 + (2*rand-1)*pi/4;
end
%delta = [pi/10 pi/2-pi/7 pi+pi/7 pi*3/2-pi/8];
N = 200;
xmax = 1;
ymax = 1;
x = linspace(-xmax,xmax,N);
y = linspace(-ymax,ymax,N);
[X,Y] = meshgrid(x,y);
A = 145;%*exp(-1.8*(X.^2+Y.^2));
B = 100;%*exp(-0.2*(X.^2+Y.^2));
phi = 2*pi*(X.^2+Y.^2);
phi_wrapped = angle(exp(j*phi));
figure(1)
surf(phi,'FaceColor','interp', 'EdgeColor','none','FaceLighting','phong');
camlight left, axis tight
xlabel('X/Pixels','FontSize',14);ylabel('Y/Pixels','FontSize',14);zlabel('Phase/Radians','FontSize',14);title('初始相位图','FontSize',14)
set(figure(1),'name','Initial Phase 3D','Numbertitle','off');
I0 = A + B.*cos(phi+delta(1));
I1 = A + B.*cos(phi+delta(2));
I2 = A + B.*cos(phi+delta(3));
I3 = A + B.*cos(phi+delta(4));
figure(2);subplot(221);imshow(I0,[]);title('四步非等步长相位图delta=0+noise');
subplot(222);imshow(I1,[]);title('四步非等步长相位图delta=pi/2+noise');
subplot(223);imshow(I2,[]);title('四步非等步长相位图delta=pi+noise');
subplot(224);imshow(I3,[]);title('四步非等步长相位图delta=3pi/2+noise');
%% *********************GPSA算法相位提取***********************
% P = K-1.*L
M = 4;
cos_delta = 0;
sin_delta = 0;
sin_cos_delta = 0;
cos_sin_delta = 0;
cos_2_delta = 0;
sin_2_delta = 0;
L1 = 0; L2 = 0; L3 = 0;
delta_G = [0 pi/2 pi 3*pi/2];
tic
for k = 1:M
cos_delta = cos_delta + cos(delta_G(k));
sin_delta = sin_delta + sin(delta_G(k));
cos_sin_delta = cos_sin_delta + cos(delta_G(k))*sin(delta_G(k));
sin_cos_delta = sin_cos_delta + sin(delta_G(k))*cos(delta_G(k));
sin_2_delta = sin_2_delta + sin(delta_G(k))^2;
cos_2_delta = cos_2_delta + cos(delta_G(k))^2;
end
K = [M cos_delta sin_delta;cos_delta cos_2_delta cos_sin_delta;
sin_delta sin_cos_
% AIA算法和GPSA算法仿真对比(随机步长,移相未知)
% 作者:James_Ray_Murphy
% 参考文献:高精度干涉测量随机移相技术研究_苏志德
%****************************%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;clc;close all;
M = 4;
for k = 1:M
delta(k) = (k-1)*pi/2 + (2*rand-1)*pi/4;
end
%delta = [pi/10 pi/2-pi/7 pi+pi/7 pi*3/2-pi/8];
N = 200;
xmax = 1;
ymax = 1;
x = linspace(-xmax,xmax,N);
y = linspace(-ymax,ymax,N);
[X,Y] = meshgrid(x,y);
A = 145;%*exp(-1.8*(X.^2+Y.^2));
B = 100;%*exp(-0.2*(X.^2+Y.^2));
phi = 2*pi*(X.^2+Y.^2);
phi_wrapped = angle(exp(j*phi));
figure(1)
surf(phi,'FaceColor','interp', 'EdgeColor','none','FaceLighting','phong');
camlight left, axis tight
xlabel('X/Pixels','FontSize',14);ylabel('Y/Pixels','FontSize',14);zlabel('Phase/Radians','FontSize',14);title('初始相位图','FontSize',14)
set(figure(1),'name','Initial Phase 3D','Numbertitle','off');
I0 = A + B.*cos(phi+delta(1));
I1 = A + B.*cos(phi+delta(2));
I2 = A + B.*cos(phi+delta(3));
I3 = A + B.*cos(phi+delta(4));
figure(2);subplot(221);imshow(I0,[]);title('四步非等步长相位图delta=0+noise');
subplot(222);imshow(I1,[]);title('四步非等步长相位图delta=pi/2+noise');
subplot(223);imshow(I2,[]);title('四步非等步长相位图delta=pi+noise');
subplot(224);imshow(I3,[]);title('四步非等步长相位图delta=3pi/2+noise');
%% *********************GPSA算法相位提取***********************
% P = K-1.*L
M = 4;
cos_delta = 0;
sin_delta = 0;
sin_cos_delta = 0;
cos_sin_delta = 0;
cos_2_delta = 0;
sin_2_delta = 0;
L1 = 0; L2 = 0; L3 = 0;
delta_G = [0 pi/2 pi 3*pi/2];
tic
for k = 1:M
cos_delta = cos_delta + cos(delta_G(k));
sin_delta = sin_delta + sin(delta_G(k));
cos_sin_delta = cos_sin_delta + cos(delta_G(k))*sin(delta_G(k));
sin_cos_delta = sin_cos_delta + sin(delta_G(k))*cos(delta_G(k));
sin_2_delta = sin_2_delta + sin(delta_G(k))^2;
cos_2_delta = cos_2_delta + cos(delta_G(k))^2;
end
K = [M cos_delta sin_delta;cos_delta cos_2_delta cos_sin_delta;
sin_delta sin_cos_