基于圆阵的music算法，想做一些改进,希望能识别出比较近的角度或是去相干，试了几种方法，但效果都不对，希望大神帮忙看一下。
clc;
clear all;
M=7; %相位模式数
N=16; %阵元数
r=1/(4*sin(pi/N)); %圆阵半径
snapshots=1024; %快拍数
SNR=20; %信噪比
targets=3; %信源数
angle=[30 50
70 20
150 25]*pi/180; %方位角
w=[pi/3 pi/4 pi/6].'; %信号频率（频率相同则相干信号）
W=zeros(N,2*M+1);
for p=1:2*M+1
for i=1:N
W(i,p)=1/N.*exp(j*2*pi*(p-8)*(i-1)/N);
end
end
%生成矩阵V
V=zeros(2*M+1,2*M+1);
for p=1:2*M+1
for q=1:2*M+1
V(p,q)=1/sqrt(2*M+1).*exp(-j*(q-8)*2*pi*(p-8)/(2*M+1));
end
end
%生成矩阵C
C=eye(2*M+1,2*M+1);
for n=1:M+1
C(n,n)=j^(-M+n-1);
end
for n=(M+2):(2*M+1)
C(n,n)=j^(M-n+1);
end
%生成矩阵F
F=V'*C*W';
%生成矩阵A
for k=1:snapshots
fai0=2*pi*rand(1,targets);%设置到达信号的随机相位
for p=1:targets
for q=1:N
beta(q)=2*pi*(q-1)/N;
A(q,p)=exp(j*(2*pi*r.*sin(angle(p,2)).*cos(angle(p,1)-beta(q))+fai0(1,p))); %方向矩阵
end
end
end
S=10.^(SNR/20)*exp(j*w*[1:snapshots]); %信号
X=A*S+(randn(N,snapshots)+i*randn(N,snapshots))/sqrt(2); %加入噪声后的信号
X=F*X;
Rx=X*X'./snapshots; %计算矩阵Y的自相关矩阵Rx
%% 原来的程序
Rx=real(Rx);
[e_vector,e_value]=eig(Rx); %对Rx进行特征分解
[EVA,I]=sort(e_value); % 调整Rx的特征值排列顺序，并调整相应的特征向量，按从大到小排列
EVA=fliplr(EVA);
e_vector=fliplr(e_vector);
G=e_vector(:,targets+1:2*M+1); %生成噪声的特征向量
%%
%%自己改的程序 根据书上的方法
J=fliplr(eye(2*M+1));
Rx=Rx+J*conj(Rx)*J;
[U,S,V1]=svd(Rx);
Vu=V1(:,targets+1:2*M+1);
S((2*M+1-targets):2*M+1,(2*M+1-targets):2*M+1)=0;
Ss=S;
R1=U*Ss*conj(V1);
[Uu,Sss,V2]=svd(R1);
Vuu=V2(:,targets+1:2*M+1);
G=(Vu+Vuu)/2;
G=real(G); %这句根据原来程序中G是实数矩阵写的
%%
%% 另一种修改方法 将原来程序的自相关矩阵重构成托普利兹矩阵 这里的15根据Rx的维度写的，也可以是2*M+1
dd=zeros(2*15-1);
for i=-(15-1):(15-1)
d=sum(diag(Rx,i))/(15-abs(i));%每一对角线取平均
dd(i+15)=d;
end
for k=1:15
Rx(k,k)=dd(15);
end
for i=1:14
for k=1:(15-i)
Rx(k+i,k)=dd(15-i);
Rx(k,k+i)=dd(15+i);
end
end
Rx=real(Rx);
[e_vector,e_value]=eig(Rx);
[EVA,I]=sort(e_value);
EVA=fliplr(EVA);
e_vector=fliplr(e_vector);
G=e_vector(:,targets+1:2*M+1);
%%
%进行二维谱峰搜索
a1=zeros(N,1);
for theta=0:90 %俯仰角
for fai=0:360 %方位角
for q=1:N
beta1(q)=2*pi*(q-1)/N;
a1(q,:)=exp(j*2*pi*r.*sin(theta/180*pi).*cos(fai/180*pi-beta1(q)));
end
b=F*a1;
P(theta+1,fai+1)=1/norm(b.'*G*G.'*b); %norm((G'*b),2);
end
end
%绘制二维功率谱图
figure(1);
xx=0:360;%方位角
yy=0:90;%俯仰角
[X1,Y1]=meshgrid(xx,yy);
mesh(X1,Y1,20*log10(P));
xlabel('方位角（度）');
ylabel('俯仰角（度）');
zlabel('空间谱（dB）');
title('均匀圆阵的二维DOA估计');
比如将第三个角度改成（70,25），原来的程序就不能很好的识别出来。大神帮忙看一下是程序哪里错了，还是这个方向错了，谢谢解答。
