我也在做扩展互质阵列的空间平滑，对着论文方法检查修改多次，角度估计有问题，希望能帮指正，互相学习。谢！

Question

clc;
clear all;
f=110^8;
c=310^8;
lamda=c/f;
d=lamda/2;
n=3;%实际
M=4; %虚拟阵元
theta=[50,70];
thita=0:1:90;
N=500;
S1=randn(1:N);%来向信号
S2=randn(1:N);
% S3= 6exp(1ipi/4*[1:N]);
% S4=6exp(1ipi/8*[1:N]);
S=[S1;S2;];
Rs = SS'/N; %协方差估计值 44
%扩展互质阵列流型矩阵
d1=0:nd:(2M-1)nd;
d2=0:Md:(n-1)Md;
A1=exp(-1i2pid1'sin(thetapi/180));
A2=exp(-1i2pi*d2'sin(thetapi/180));
A=[A1;A2];

%实际协方差矩阵
Rxx=ARsA';
Rxx=awgn(Rxx,10);
%互质阵列对应的虚拟阵列的均匀线性子阵
Aa=-Mnd:d:Mnd;
AA=[];
%%得到均匀线性阵列所对应的流型矩阵
%%差集数组
[mat1,mat2]=ndgrid([0 3d 6d 9d 12d 15d 18d 21d 0 4d 8d]);
R_b=mat2-mat1;
%%2Mn+1个行向量重新排序
for i=1:1:(2Mn+1)
aA=find(R_b==Aa(i));
AA=[AA;aA(1)];
end
Z1=Rxx(AA);
R=zeros(13);
%%空间平滑协方差矩阵
for i=1:1:Mn+1
Zi=Z1(Mn+2-i:2Mn+2-i,:);
Ri= Zi Zi';
R=R+Ri;
end
R1=R/(Mn+1);
R=zeros(13);
for i=1:1:Mn+1
Zi=Z1(i:Mn+i,:);
Ri= Zi Zi';
R=R+Ri;
end
R2=R/(Mn+1);
R=(R1+R2)/2;
%%特征值，特征向量，噪声子空间
[V,D]=eig(R);
[Y,I]=sort(diag(D));
Un=V(:,I(1:3));
%%
dd=0:d:Mnd;
for i=1:length(thita)
Al=[exp(-1i2pidd*sin(thita(i)*pi/180)/lamda)]';
Pmusics(i)=abs(1/(Al'UnUn'*Al));
end

hold on;
plot(thita,10*log(Pmusics),'b--');
xlabel('入射方向角（度）');
ylabel('输出功率（dB）');