function [esta,estb,estscales,estA,estB] = MDE(X2,F);

% Multi-Dimensional Embedding algorithm for 2-D damped harmonic retrieval
% References: 
% N.D. Sidiropoulos, X. Liu, and A. Swami, A New 2-D Harmonic Retrieval Algorithm,
% in Proc. 39th Allerton Conference on Communication, Control, and Computing,
% Oct. 3-5, 2001 (MDE implements the algebraic part only, not ALS refinement)
% T. Jiang, N.D. Sidiropoulos, and J. M.F. ten Berge, Almost Sure Identifiability
% of Multidimensional Harmonic Retrieval, IEEE Trans. on Signal Processing, 49(9):1849-
% 1859, Sep. 2001.
% Developed by N.D. Sidiropoulos and Tao Jiang
% Last modified by N. Sidiropoulos/ECE/UMN, Jan 25, 2001.
%
% Input: 
%
% X2: 2-D mixture of damped harmonics
% F: number of components
%
% Output:
%
% esta: sorted estimated cplx generators along the row dimension
% estb: sorted estimated cplx generators along the column dimension
%         (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(X2);

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

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

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

% construct 5-way array:

X5 = zeros(I1,I2,I3,J1,J2);

for i1=1:I1,
 for i2=1:I2,
  for i3=1:I3,
   for j1=1:J1,
    for j2=1:J2,
     X5(i1,i2,i3,j1,j2) = X2(i1+i2+i3-2,j1+j2-1);
    end
   end
  end
 end
end

% nest into 3-way array:

X3 = zeros(I3,I1*J1,I2*J2);

for i3=1:I3,
 for k=1:I1*J1,
  for l=1:I2*J2,
   X3(i3,k,l) = X5(ceil(k/J1),ceil(l/J2),i3,k-(ceil(k/J1)-1)*J1,l-(ceil(l/J2)-1)*J2);
  end
 end
end

[D,E,ratio]=betGE(squeeze(X3(1,:,:)),squeeze(X3(2,:,:)),F);
ratios = zeros(1,F);
for f=1:F,
 ratios(f) = ratio(f,f);
end
A3 = [zeros(1,F); ratios];

%%%% Estimates can be significantly improved,
%%%% by averaging out all available ratios:

%%%% Recovered from D

estaD = 0;
estbD = 0;

for k=1:J1*(I1-1)
   estaD = estaD+(D(J1+k,:)./D(k,:))';
end

for k=0:I1-1
   for l=1:J1-1
      estbD = estbD+(D(l+k*J1+1,:)./D(l+k*J1,:))';
   end
end

estaD = estaD/(I1*J1-J1);
estbD = estbD/(I1*J1-I1);

%%%% Recovered from E

E=E';

estaE = 0;
estbE = 0;

for k=1:J2*(I2-1)
   estaE = estaE+conj((E(J2+k,:)./E(k,:))');
end

for k=0:I2-1
   for l=1:J2-1
      estbE = estbE+conj((E(l+k*J2+1,:)./E(l+k*J2,:))');
   end
end

estaE = estaE/(I2*J2-J2);
estbE = estbE/(I2*J2-I2);

E=E';

%%%% Recovered from A3

estaA3=A3(2,:)';

%%%% Average process

esta = ((estaD+estaE+estaA3)/3)';
estb = ((estbD+estbE)/2)';

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

for f=1:F,
 if aesta(f) < 0
  aesta(f) = aesta(f) + 2*pi;
 end 
 if aestb(f) < 0
  aestb(f) = aestb(f) + 2*pi;
 end 
end

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

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

for f=1:F,
 for i=1:I,
  estA(i,f) = esta(f).^(i-1);
 end
end

for f=1:F,
 for j=1:J,
  estB(j,f) = estb(f).^(j-1);
 end
end

estscales = pinv(krp(estB,estA))*vec(X2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of main function body; auxiliary functions follow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function vectorized = vec(A);

[R,C]=size(A);

vectorized = [];

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

function AkrpB = krp(A,B);

[I F] = size(A);
[J F1] = size(B);

if (F1 ~= F)
 disp('krp.m: column dimensions do not match!!! - exiting matlab');
 exit;
end

AkrpB = [];
for f=1:F,
 AkrpB = [AkrpB kron(A(:,f),B(:,f))];
end

function [estB,estC,estD] = betGE(X1,X2,M)

% M = QF the rank of estB and estC.
% X1 = B*C.'   X2 = B*D*C.'                 Eq. (1)
% estB: signal matrix
% estC: effective code matrix
% estD: D = D2*inv(D1)

[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,:);

R1 = U1.'*U1; % = T*B.'*B*T.' = G*T.'       Eq. (3)
R2 = U1.'*U2; % = T*B.'*B*D*T.' = G*D*T.'   Eq. (4)

% inv(G)*R1 = inv(D)*inv(G)*R2 <=> D*inv(G) = inv(G)*R2*inv(R1)
% Let Gamma = inv(G) and A = R2*inv(R1). => D*Gamma = Gamma*A
%
% In MATLAB [V,D] = eig(X) means X*V = V*D. So, transpose the above
% equation. A.'*Gamma.' = Gamma.'*D (since D is a diagonal matrix).

A = R2*inv(R1);
[Gamma,estD] = eig(A.');
Gamma = Gamma.';

% Eq. (3) => T.' = inv(G)*R1 = Gamma*R1.
% Eq. (2) => B = U1*inv(T.') = U1*inv(Gamma*R1)
estB = U1*inv(Gamma*R1);

estC = pinv(estB)*X1; % from Eq. (1), but here estC is C.'.

% sort estD, estB and estC

   clear d de p;
   d = diag(estD);
   [de,p] = sort(d);
   clear pB pC;
   t = 1:1:M; % M = Q*F
   t = t.';
   estD = diag(de);
   pB(:,t) = estB(:,p);
   estB = pB;
   pC(t,:) = estC(p,:);
   estC = pC;