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