clear;load lena.mat; 
XR=X(121:141,121:141);XI=X(121:141,151:171);
X=XR+j*XI;M=size(X,1);N=ceil(2.6*M);MN=(M-1)*N;
Y=conv2(X,fliplr(flipud(conj(X))));%Same as this:
%Y=fftshift(ifft2((abs(fft2(X,2*M-1,2*M-1))).^2))
FY1=fft(Y,N);FY2=fft(Y.',N);for I=1:N;
RY1(:,I)=sort(roots(FY1(I,:)));
RY2(:,I)=sort(roots(FY2(I,:)));
end;%RY?(:,I) in conjugate reciprocal pairs
FX1=fft(X,N);FX2=fft(X.',N);for I=1:N;
RX1(:,I)=sort(roots(FX1(I,:)));
RX2(:,I)=sort(roots(FX2(I,:)));
end;for I=1:N;for J=1:M-1;
FR1=unwrap(angle((ones(1,N)-RY1(J,I)*exp(-j*2*pi*[0:N-1]/N)))+pi*[0:N-1]/N)';
A((I-1)*N+1:I*N,(J-1)*N+I)=FR1(1:N);
FR2=unwrap(angle((ones(1,N)-RY2(J,I)*exp(-j*2*pi*[0:N-1]/N)))+pi*[0:N-1]/N)';
A(I:N:N*N,(J-1)*N+I+MN)=-FR2(1:N);end;end
for I=0:N-1;B([1:N-1]+(N-1)*I,:)=A([2:N]+N*I,:)-A([1:N-1]+N*I,:);end;
for I=0:N-2;C([1:N-1]+(N-1)*I,:)=B([1:N-1]+(N-1)*(I+1),:)-B([1:N-1]+(N-1)*I,:);end;
CC=C';DD=CC(:,sum(abs(CC))<6.28);D=DD';
X1=D(:,1:2*MN-1)\D(:,2*MN);
%[U S V]=svd(D);X1=V(:,2*MN)/V(1,2*MN);
P1=A(:,1:MN)*X1(1:MN);
Q1=reshape(P1,N,N)-pi*(M-1)*[0:N-1]'*ones(1,N)/N;
P2=A(:,MN+1:2*MN)*[X1(MN+1:2*MN-1);-1];
Q2=-reshape(P2,N,N)-pi*(M-1)*ones(N,1)*[0:N-1]/N;
Q3=exp(j*Q1)./exp(j*Q2);[U1 S1 V1]=svd(Q3);%rank=1
P=reshape(P1,N,N)-pi*(M-1)*[0:N-1]'*ones(1,N)/N+ones(N,1)*(angle(V1(:,1)))';
EST=ifft2(abs(fft2(X,N,N)).*exp(j*P'));EST=EST(1:M,1:M)/EST(1,1)*X(1,1);
subplot(661),imagesc(real(X)),colormap(gray),axis off,title('ORIG')
subplot(662),imagesc(real(EST)),colormap(gray),axis off,title('RECON')
subplot(334),imagesc(abs(Y)),colormap(gray),axis off,title('AUTOCORRE')
subplot(667),imagesc(imag(X)),colormap(gray),axis off
subplot(668),imagesc(imag(EST)),colormap(gray),axis off
