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;