function [saesta, saestb, estscales,estA,estB] = MDF(X,F);

% Multi-Dimensional Folding algorithm for 2-D harmonic retrieval
% References: 
% X. Liu, and N.D. Sidiropoulos, Almost Sure Identifiability of Constant Modulus
% Multidimensional Harmonic Retrieval, IEEE Trans. on Signal Processing, 
% 50(9):2366-2368, Sep. 2002
% X. Liu, and N.D. Sidiropoulos, On Constant Modulus Multidimensional Harmonic
% Retrieval, in Proc. ICASSP2002, Orlando, Florida, May 13-17, 2002.
% Developed by X. Liu, N.D. Sidiropoulos, with contributions by Tao Jiang
% Last modified by X. Liu, ECE UofL, Nov. 8, 2004

% Input: 
%
% X: 2-D harmonic mixture
% F: number of components
%
% Output:
%
% saesta: sorted estimated frequencies along the row dimension in rads
% saestb: sorted estimated frequencies along the column dimension in rads
%         (sorting is according to est freqs along the row dimension)
% estscales: estimated complex amplitudes 
% estA, estB: reconstructed loading matrices (X ~ A C B.', where C is 
% diagonal, holding complex amplitudes)

[I,J]=size(X);

% select I1, I2, I3, J1, J2, according to Theorem:

if mod(I,2) == 0
 I1 = I/2;
 I2 = (I+2)/2;
else
 I1 = (I+1)/2;
 I2 = (I+1)/2;
end

if mod(J,2) == 0
 J1 = J/2;
 J2 = (J+2)/2;
else
 J1 = (J+1)/2;
 J2 = (J+1)/2;
end

X2 = zeros(I1*J1,I2*J2);
%%%%%%%%%%%%% An alternative way to form X2
for i1=1:I1
    for i2=1:I2
        for j1=1:J1
            for j2=1:J2
                   X2((i1-1)*J1+j1,(i2-1)*J2+j2)=X(i1+i2-1, ...
                    j1+j2-1);
            end
        end
    end
end

%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%


Y2=flipud(fliplr(conj(X2)));

[G,H,ratio]=LiuGE(X2,Y2,F);
ratios=diag(ratio);

% ratios contain e^{-j\omega_f(I-1)-j\nu_f(J-1)}

estaG = zeros(1,F);    
estbG = zeros(1,F);

estaG= mean(G(J1+[1:J1*(I1-1)],:)./G(1:J1*(I1-1),:));
%%%%%%%%%%%%%%%%%%%%%%
% From X. Liu, Nov. 8, 2004:
% When J1=1,  1:J1-1=1:1,  then G([1:J1-1]+k*J1+1,:) becomes a vector
% "sum" will cause problem, so use iteration as follows. 
%%%%%%%%%%%%%%%%%%%%%%
%for k=0:I1-1
%   estbG = estbG+sum(G([1:J1-1]+k*J1+1,:)./G([1:J1-1]+k*J1,:));
%end
for k=0:I1-1
    sum1=zeros(1,F);
    for jj=1:J1-1
    sum1=sum1+G(jj+k*J1+1,:)./G(jj+k*J1,:);
    end
    estbG = estbG+sum1;
end

estbG = estbG/(I1*J1-I1);

estaH = zeros(1,F);
estbH = zeros(1,F);

estaH= mean(H(J2+[1:J2*(I2-1)],:)./H(1:J2*(I2-1),:));
%for k=0:I2-1
%   estbH = estbH+sum(H([1:J2-1]+k*J2+1,:)./H([1:J2-1]+k*J2,:));
%end
for k=0:I2-1
    sum1=zeros(1,F);
    for jj=1:J2-1
    sum1=sum1+H(jj+k*J2+1,:)./H(jj+k*J2,:);
    end
    estbH = estbH+sum1;
end
estbH = estbH/(I2*J2-I2);


%%%% Average process

esta = (estaG+estaH)/2;
estb = (estbG+estbH)/2;

aesta = angle(esta);
aestb = angle(estb);


[saesta,estperm] = sort(aesta);
saestb = aestb(estperm);

esta = esta(estperm);
estb = estb(estperm);

for k=1:I,
	estA(k,:)=esta.^(k-1);
end;

for l=1:J,
	estB(l,:)=estb.^(l-1);
end;


BkrpA = [];
for f=1:F,
 BkrpA = [BkrpA kron(estB(:,f),estA(:,f))];
end

estscales = (pinv(BkrpA)*vec(X)).';

% Function LiuGE------------------------------------------

function [estB,estC,estD] = LiuGE(X1,X2,M)
[U,S,V] = svd([X1; X2]);
U = U(:,1:M);
S = S(1:M,1:M);
V = V(:,1:M);

%     [ U1 ]   [ B  ]
% U = [- - ] = [- - ] T.', where T:QF x QF  Eq. (2)
%     [ U2 ]   [ BD ]

[nrowsU,ncolsU] = size(U);
U1 = U(1:(nrowsU/2),:);
U2 = U(1+(nrowsU/2):nrowsU,:);


A=pinv(U1)*U2;
[invT,estD] = eig(A);
estB = U1*invT;
estC = (pinv(estB)*X1).'; 

% Function Vec-------------------------------------------
function vectorized = vec(A);

[R,C]=size(A);

vectorized = [];

for c=1:C,
 vectorized = [vectorized; A(:,c)];
end