function testJStB

% Corroborate AS Uniqueness of 2-D HR claim in Jiang-Sidiropoulos-ten Berge
% 
% Reference: 
% 
% 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 (see also Proc. ICASSP2001)
% Developed by N.D. Sidiropoulos
% NOTE: This version does not incorporate averaging of estimates (for this, 
% use MDE.m) - present program is simply to illustrate the ID theorem. 
% Last modified by N. Sidiropoulos/ECE/UMN, Jan 25, 2001.

clear all;

% generate 2-D harmonic mixture:

%F=9;
%I=6;
%J=6;

% try something a bit more challenging:

F=16;
I=8;
J=8;

% test w/ unit-circle generators for simplicity, but works as is for 
% damped harmonics as well:

a = exp(sqrt(-1)*2*pi*rand(F,1));
b = exp(sqrt(-1)*2*pi*rand(F,1));

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

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

scales = randn(F,1)+sqrt(-1)*randn(F,1);

X2 = A*diag(scales)*B.';

% 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];

esta = D(J1+1,:)./D(1,:);

estb = D(2,:)./D(1,:);

aa = angle(a);
ab = angle(b);
aesta = angle(esta);
aestb = angle(estb);

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

[saa,perm] = sort(aa);
sab = ab(perm);
[saesta,estperm] = sort(aesta);
saestb = aestb(estperm);

disp('Column frequencies, true');
saa.'
disp('Column frequencies, recovered');
saesta
disp('Row frequencies, true');
sab.'
disp('Row frequencies, recovered');
saestb

a = a(perm);
b = b(perm);
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);

disp('Scale factors, true');
scales = scales(perm).'
scales
disp('Scale factors, recovered');
estscales = estscales.'
estscales 

figure(1);
plot(a,'o');
hold on
plot(esta,'*r');
axis([-1 1 -1 1]);
xlabel('real');
ylabel('imag');
legend('True column frequencies','Recovered column frequencies');
hold off;
figure(2);
plot(b,'o');
hold on
plot(estb,'*r');
axis([-1 1 -1 1]);
xlabel('real');
ylabel('imag');
legend('True row frequencies','Recovered row frequencies');
hold off;
figure(3);
plot(scales,'o');
hold on
plot(estscales,'*r');
xlabel('real');
ylabel('imag');
legend('True scales','Recovered scales');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;