Joint SVD

%2017.6.24%JSVD%Joint SVD of 2 cross_correlation matrices to achieve automatic paring in%2-D angle estimation problemsclcclose all;clear allMx=10;                                                                     %x轴阵元数Mz=10;                                                                     %z轴阵元数d_lamda=1/2;                                                               %阵元间距K=2;                                                                       %信号数phix=[45,52];                                                              %x轴方位角thetaz=[15,25];                                                            %z轴俯仰角fc=[200*10^6,300*10^6];                                                    %信号频率fs=1000*10^6;                                                              %采样频率L=1000;                                                                    %快拍数snr=10;                                                                    %信噪比for k=1:L    S(:,k)=sqrt(10.^(snr/10))*randn(1)*exp(1j*2*pi*fc'/fs*(k-1));endfor kk=1:K    A1x(:,kk)=exp(-1j*2*pi*[0:Mx-1]'*d_lamda*cos(phix(kk)/180*pi));    A1z(:,kk)=exp(-1j*2*pi*[0:Mz-1]'*d_lamda*cos(thetaz(kk)/180*pi));      %1维导向矢量endZZ=A1z*S+wgn(Mz,L,0,'dBm','complex');                                      %z轴接收到的信号XX=A1x*S+wgn(Mz,L,0,'dBm','complex');                                      %x轴接收到的信号Rxz=XX*ZZ'/L;                                                              %(4)获取协方差矩阵Rxz1=Rxz(:,1:Mx-1);Rxz2=Rxz(:,2:Mx);R=[Rxz1;Rxz2];                                                             %(10)[U,D,V]=svd(R);                                                            %(11)U1=U(:,1:K);U11=U1(1:Mx,:);U12=U1(Mx+1:2*Mx,:);W=pinv(U11)*U12;[T,w]=eig(W);                                                              %(18)gama=(diag(w)).';                                                          %(20)doa_theta=acos(angle(gama)/pi)*180/pi;doa_theta=sort(doa_theta,'ascend');                                        %估计出来方位角Aphi=U11*(T);phii=0:0.0001:90;for k=1:Kfor  kk=1:length(phii)      a=exp(-1j*2*pi*d_lamda*cos(phii(kk)/180*pi)*[0:Mx-1]');   p(k,kk)=(abs(a'*Aphi(:,k)))/(norm(a,2).^2);  end   [Pv(k,:),Pd(k,:)]=max(p(k,:));   doa_phi(k)=(Pd(k,:)-1)/10000;enddisp('（俯仰角，方位角）');for k=1:K    disp([num2str(doa_phi(k)) ',' num2str(doa_theta(k))]);end