DOA——MUSIC算法

一、均匀圆阵(UCA, Uniform Circular Array)的MUSIC算法

假设一个半径为R的M元均匀圆阵的所有阵元均位于坐标系X-Y平面内，第k-1个阵元坐标为，第i个窄带信号波长为，来波方向为，如图1，则第k-1个阵元到圆心(即原点)的波程差为：

均匀圆阵

存在P个入射信号均匀圆阵的接收模型可表示为：

其他步骤与基于ULA的MUSIC算法一致。

令任意两阵元间的波程差为：

当时，即产生相位模糊。将均匀圆阵各阵元投影到入射方向，得到一个随入射方向变动的非均匀线阵。需要保证在任意入射方向上投影出的非均匀线阵，其最小间隔总是小于信号波长，模糊谱峰对测向结果影响较小，即：

在方向角相同时，水平入射()信号的波程差最长，且投影出的非均匀线阵随方向角不同周期变化，因此只需要讨论水平入射信号对应投影线阵的不同情况。

在MUSIC算法中，阵元的最小间隔越大模糊谱峰峰值就越大。但在均匀圆阵中，阵元间隔随着入射波方向变化，因此算法性能受到最小间隔最大值的影响。根据来波方向不同，入射方向上的第k个阵元投影间隔分别为：

当M为奇数时，对着阵元入射，投影点重合为${\rm{(M + 1)/2}}$个；当M为偶数时，对着相邻阵元连线中点入射，投影点重合为M/2个，此时投影线阵的非零最小间隔的值最大，且取得该最大值k=1时。进一步可求得半径的选取关系：

选取半径时，按上式等比例缩放，即能使对应的奇数阵与偶数阵有近似的抗相位模糊的性能。

每一路接收的结构图：

二、MUSIC

1、

clear all;
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsym=500;%符号个数
Fsym=1;%符号速率
M=3;%一个符号对应的比特数
Fbit=M*Fsym;%比特速率
Nsour=3;%信源数
Angle=[5,15,35];%信源的来波方向
Fc=10;%载波频率
Fs=100;%抽样频率
R=0.5;%滚降因子
Del=5;%群延迟因子
% Nsamp=50;%采样点数或者快拍数
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%构造相干信源--信源1、信源2与信源3
Rcos2=Rcos11+Rcos21+Rcos31;%构造相干信源--信源1、信源2与信源3
Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相干信源--信源1、信源2与信源3
save xyc3 Rcos1 Rcos2 Rcos3
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsamp=1024;%采样点数或者快拍数 
i=sqrt(-1);
j=i;
Ntx=8;%阵列数
SNR=[2,2,2];%三信源的信噪比
% sn=10;                        %----单信号源
Lamda=2;%信号波长
D=Lamda/2;%线性阵列的距离
p=3;%子阵个数
L=Ntx-p+1;%子阵阵元数
nr=randn(Ntx,Nsamp);
ni=randn(Ntx,Nsamp);
n=nr+j*ni;%产生背景噪声
load xyc3;
t=1:Nsamp;
% s1=[Rcos1(t).'];%接收信号的采样点数%----单信号源
s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩阵维数=信源数*抽样点数
ps=diag((s1*s1')/Nsamp);%无噪声信号功率--%矩阵维数=信源数*1
delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1
% delta1=ps./(2*10.^(sn/10));        %----单信号源
delta2=diag(delta1);%矩阵维数=1*1
delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1
Rev_s1=(1./delta')*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数
%计算各信源SNR比条件下，阵列接收到的信号幅度%
Pn=zeros(Nsamp,1);
pn=zeros(Ntx,Nsamp);
Pn=diag(n'*n);
for h=1:Nsamp
    pn(:,h)=n(:,h)./sqrt(Pn(h,:));
end
Rev_n=pn;
%计算各阵列接收到的背景噪声下的信号幅度%
tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数
% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda;   %----单信号源
tmp1=[0:Ntx-1]';%矩阵维数=阵元数*1
tmp4=[0:L-1]';%子矩阵维数=子矩阵阵元数*1
a1=tmp1*tmp;%矩阵维数=阵元数*信源数
A=exp(a1);%方向矩阵--%矩阵维数=阵元数*信源数
X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数
Rxx=(X*X')/Nsamp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sub_FRxx=zeros(L,L);
sub_BRxx=zeros(L,L);
for i=1:p
    sub_FR=zeros(L,Nsamp);
    sub_BR=zeros(L,Nsamp);
    sub_FR=X(i:1:i+L-1,:);
    last=Ntx+1-i;
    first=last-L+1;
    sub_BR=conj(X(last:-1:first,:));
    sub_FRxx=sub_FRxx+((sub_FR*sub_FR')./Nsamp);
    sub_BRxx=sub_BRxx+((sub_BR*sub_BR')./Nsamp);
end
sub_FRxx=sub_FRxx./p;
sub_BRxx=sub_BRxx./p;
sub_Rxx=(sub_FRxx+sub_BRxx)./2;
[VFB,HFB]=eig(sub_Rxx);
[HFB,IFB]=sort(diag(HFB),1);
VFB=VFB(:,IFB);
VnFB=VFB(:,1:L-Nsour);
VsFB=VFB(:,L-Nsour+1:L);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ScanAng=[-90:1:90];
for i=1:length(ScanAng)
    tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;
    tmp3=tmp2*tmp1;
    tmp5=tmp2*tmp4;
    A_Sita=exp(tmp3);
    Sub_Sita=exp(tmp5);
    Sub_FBsita(i)=(Sub_Sita'*Sub_Sita)/(Sub_Sita'*VnFB*VnFB'*Sub_Sita);
end 
figure(1);
semilogy(ScanAng,real(Sub_FBsita),'bo-');
axis([-60 60 0.1 1e7]);
xlabel('M_Angle(deg)');
ylabel('M_Spectrum');
grid on

　　

2、

clear all;
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsym=500;%符号个数
Fsym=1;%符号速率
M=3;%一个符号对应的比特数
Fbit=M*Fsym;%比特速率
Nsour=3;%信源数
Angle=[10,40,80];%信源的来波方向
Fc=10;%载波频率
Fs=100;%抽样频率
R=0.5;%滚降因子
Del=5;%群延迟因子
% Nsamp=50;%采样点数或者快拍数
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=Rcos11;%构造相干信源--信源1、信源2与信源3
Rcos2=Rcos21;%构造相干信源--信源1、信源2与信源3
Rcos3=Rcos31;%构造相干信源--信源1、信源2与信源3
save xyc3 Rcos1 Rcos2 Rcos3
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsamp=1024;%采样点数或者快拍数 
i=sqrt(-1);
j=i;
Ntx=8;%阵