Code for Email 2016 IETSP_克拉-美罗界_MUSIC算法_克拉美罗_克拉美罗界_

function [Phi_Est,Theta_Est]=SSJD_NLSA(X,Z,K,X_Location,L,flag) % flag=1 2 3, they have the same performance. [M,N]=size(X); C=kr_sum(-X_Location,X_Location); if mod(M,2)~=0 M1=(M-1)/2;M2=(M+1)/2; DOF=(M^2-1)/2+M; else M1=M/2;M2=M/2; DOF=(M^2-2)/2+M; end M0=(DOF+1)/2; J=zeros(DOF,M^2); for m=1:DOF J(m,C==(m-M0))=1; end Js=sum(J,2); for m=1:DOF J(m,:)=J(m,:)/Js(m); end Rxx=zeros(M^2,2*L-1); Rzz=zeros(M^2,2*L-1); Rxz=zeros(M^2,2*L-1); for jj=1:L for m2=1:M for m1=1:M Rxx(m1+(m2-1)*M,L-jj+1)=X(m1,1:1+N-L)*X(m2,jj:jj+N-L)'/(N-L+1); Rxx(m1+(m2-1)*M,jj+L-1)=X(m1,jj:jj+N-L)*X(m2,1:1+N-L)'/(N-L+1); Rzz(m1+(m2-1)*M,L-jj+1)=Z(m1,1:1+N-L)*Z(m2,jj:jj+N-L)'/(N-L+1); Rzz(m1+(m2-1)*M,jj+L-1)=Z(m1,jj:jj+N-L)*Z(m2,1:1+N-L)'/(N-L+1); Rxz(m1+(m2-1)*M,L-jj+1)=X(m1,1:1+N-L)*Z(m2,jj:jj+N-L)'/(N-L+1); Rxz(m1+(m2-1)*M,jj+L-1)=X(m1,jj:jj+N-L)*Z(m2,1:1+N-L)'/(N-L+1); end end end Rxx_hat=J*Rxx; Rzz_hat=J*Rzz; J2=fliplr(eye(2*L-1)); Rxz_hat=conj(Rxz)*J2; Y=[Rxx_hat;Rzz_hat;Rxz;Rxz_hat]; if flag==1 R=Y*Y'/(2*L-1); [U,D]=eig(R); [d,ind]=sort(diag(D),'descend'); D=diag(d); U=U(:,ind); Us=U(:,1:K); Usx1=zeros(DOF-1+2*(M-1)*M1,K); Usx2=zeros(DOF-1+2*(M-1)*M1,K); Usz1=zeros(DOF-1+2*(M-1)*M1,K); Usz2=zeros(DOF-1+2*(M-1)*M1,K); Usx1(1:DOF-1,:)=Us(1:DOF-1,:); Usx2(1:DOF-1,:)=Us(2:DOF,:); Usz1(1:DOF-1,:)=Us(DOF+2:2*DOF,:); Usz2(1:DOF-1,:)=Us(DOF+1:2*DOF-1,:); for m=1:M-1 Usx1(DOF-1+(m-1)*M1+1:DOF-1+m*M1,:)=... Us(2*DOF+m*M+1:2*DOF+m*M+M1,:); Usx1(DOF-1+(M-1)*M1+(m-1)*M1+1:DOF-1+(M-1)*M1+m*M1,:)=... Us(2*DOF+M^2+m*M+2:2*DOF+M^2+m*M+M1+1,:); Usx2(DOF-1+(m-1)*M1+1:DOF-1+m*M1,:)=... Us(2*DOF+m*M+2:2*DOF+m*M+M1+1,:); Usx2(DOF-1+(M-1)*M1+(m-1)*M1+1:DOF-1+(M-1)*M1+m*M1,:)=... Us(2*DOF+M^2+m*M+1:2*DOF+M^2+m*M+M1,:); Usz1(DOF-1+(m-1)*M1+1:DOF-1+m*M1,:)=... Us(2*DOF+m+1:M:2*DOF+(M1-1)*M+m+1,:); Usz1(DOF-1+(M-1)*M1+(m-1)*M1+1:DOF-1+(M-1)*M1+m*M1,:)=... Us(2*DOF+M^2+M+m+1:M:2*DOF+M^2+M1*M+m+1,:); Usz2(DOF-1+(m-1)*M1+1:DOF-1+m*M1,:)=... Us(2*DOF+M+m+1:M:2*DOF+M1*M+m+1,:); Usz2(DOF-1+(M-1)*M1+(m-1)*M1+1:DOF-1+(M-1)*M1+m*M1,:)=... Us(2*DOF+M^2+m+1:M:2*DOF+M^2+(M1-1)*M+m+1,:); end G1=pinv(Usx1*D(1:K,1:K)^0.5)*Usx2*D(1:K,1:K)^0.5; H1=pinv(Usz1*D(1:K,1:K)^0.5)*Usz2*D(1:K,1:K)^0.5; [V,D]=joint_diag([G1,H1],1e-10); Psi=diag(D(1:K,1:K));Ksi=diag(D(1:K,K+1:2*K)); Phi_Est=acosd(angle(Psi)/pi); Theta_Est=acosd(-angle(Ksi)/pi); [Theta_Est,ind]=sort(Theta_Est); Phi_Est=Phi_Est(ind); elseif flag==2 Y_hat=zeros(2*DOF+2*(M^2-2*M+1),2*L-1); Y_hat(1:2*DOF,:)=Y(1:2*DOF,:); for m=1:M-1 Y_hat(2*DOF+(m-1)*(M-1)+1:2*DOF+(m-1)*(M-1)+M-1,:)... =Y(2*DOF+M+(m-1)*M+2:2*DOF+M+(m-1)*M+M,:); Y_hat(2*DOF+(M^2-2*M+1)+(m-1)*(M-1)+1:2*DOF+(M^2-2*M+1)+(m-1)*(M-1)+M-1,:)=... Y(2*DOF+M^2+M+(m-1)*M+2:2*DOF+M^2+M+(m-1)*M+M,:); end R=Y_hat*Y_hat'/(2*L-1); [U,D]=eig(R); [d,ind]=sort(diag(D),'descend'); D=diag(d); U=U(:,ind); Us=U(:,1:K); Usx1=zeros(DOF-1+2*(M-1)*(M1-1),K); Usx2=zeros(DOF-1+2*(M-1)*(M1-1),K); Usz1=zeros(DOF-1+2*(M-1)*(M1-1),K); Usz2=zeros(DOF-1+2*(M-1)*(M1-1),K); Usx1(1:DOF-1,:)=Us(1:DOF-1,:); Usx2(1:DOF-1,:)=Us(2:DOF,:); Usz1(1:DOF-1,:)=Us(DOF+2:2*DOF,:); Usz2(1:DOF-1,:)=Us(DOF+1:2*DOF-1,:); for m=1:M-1 Usx1(DOF-1+(m-1)*(M1-1)+1:DOF-1+m*(M1-1),:)=... Us(2*DOF+(m-1)*(M-1)+1:2*DOF+(m-1)*(M-1)+M1-1,:); Usx1(DOF-1+(M-1)*(M1-1)+(m-1)*(M1-1)+1:DOF-1+(M-1)*(M1-1)+m*(M1-1),:)=... Us(2*DOF+(M^2-2*M+1)+(m-1)*(M-1)+2:2*DOF+(M^2-2*M+1)+(m-1)*(M-1)+M1,:); Usx2(DOF-1+(m-1)*(M1-1)+1:DOF-1+m*(M1-1),:)=... Us(2*DOF+(m-1)*(M-1)+2:2*DOF+(m-1)*(M-1)+M1,:); Usx2(DOF-1+(M-1)*(M1-1)+(m-1)*(M1-1)+1:DOF-1+(M-1)*(M1-1)+m*(M1-1),:)=... Us(2*DOF+(M^2-2*M+1)+(m-1)*(M-1)+1:2*DOF+(M^2-2*M+1)+(m-1)*(M-1)+M1-1,:); Usz1(DOF-1+(m-1)*(M1-1)+1:DOF-1+m*(M1-1),:)=... Us(2*DOF+m:M-1:2*DOF+(M1-2)*(M-1)+m,:); Usz1(DOF-1+(M-1)*(M1-1)+(m-1)*(M1-1)+1:DOF-1+(M-1)*(M1-1)+m*(M1-1),:)=... Us(2*DOF+(M^2-2*M+1)+M-1+m:M-1:2*DOF+(M^2-2*M+1)+(M1-1)*(M-1)+m,:); Usz2(DOF-1+(m-1)*(M1-1)+1:DOF-1+m*(M1-1),:)=... Us(2*DOF+M-1+m:M-1:2*DOF+(M1-1)*(M-1)+m,:); Usz2(DOF-1+(M-1)*(M1-1)+(m-1)*(M1-1)+1:DOF-1+(M-1)*(M1-1)+m*(M1-1),:)=... Us(2*DOF+(M^2-2*M+1)+m:M-1:2*DOF+(M^2-2*M+1)+(M1-2)*(M-1)+m,:); end G1=pinv(Usx1*D(1:K,1:K)^0.5)*Usx2*D(1:K,1:K)^0.5; H1=pinv(Usz1*D(1:K,1:K)^0.5)*Usz2*D(1:K,1:K)^0.5; [V,D]=joint_diag([G1,H1],1e-10); Psi=diag(D(1:K,1:K));Ksi=diag(D(1:K,K+1:2*K)); Phi_Est=acosd(angle(Psi)/pi); Theta_Est=acosd(-angle(Ksi)/pi); [Theta_Est,ind]=sort(Theta_Est); Phi_Est=Phi_Est(ind); elseif flag==3 Ryy=Y*Y'/(2*L-1); [U,D]=eig(Ryy); [d,ind]=sort(diag(D),'descend'); D=diag(d); U=U(:,ind); Us=U(:,1:K); H01=[eye(DOF-1),zeros(DOF-1,DOF+1)]; H02=[zeros(DOF-1,1),eye(DOF-1),zeros(DOF-1,DOF)]; H03=[zeros(DOF-1,DOF+1),eye(DOF-1)]; H04=[zeros(DOF-1,DOF),eye(DOF-1),zeros(DOF-1,1)]; H1=kron(eye(M),[eye(M1),zeros(M1,M-M1)]); H2=kron(eye(M),[zeros(M1,1),eye(M1),zeros(M1,M-M1-1)]); H3=kron([eye(M1),zeros(M1,M-M1)],eye(M)); H4=kron([zeros(M1,1),eye(M1),zeros(M1,M-M1-1)],eye(M)); G1=blkdiag(H01,H1,H2); G2=blkdiag(H02,H2,H1); G3=blkdiag(H04,H4,H3); G4=blkdiag(H03,H3,H4); Q_phi=pinv(G1*Us*D(1:K,1:K)^0.5)*G2*Us*D(1:K,1:K)^0.5; Q_theta=pinv(G3*Us*D(1:K,1:K)^0.5)*G4*Us*D(1:K,1:K)^0.5; [V,D]=joint_diag([Q_phi,Q_theta],1e-10); Psi_phi=diag(D(1:K,1:K));Psi_theta=diag(D(1:K,K+1:2*K)); Phi_Est=acosd(angle(Psi_phi)/pi); Theta_Est=acosd(angle(Psi_theta)/pi); [Theta_Est,ind]=sort(Theta_Est); Phi_Est=Phi_Est(ind); end end