musicsf.m

近场MUSIC算法

clear all;clc;
M=5;%阵列个数
P=1;%信号个数
N=1024;%快拍数
C=3.0E+8;
f0=1.0E+9;%载波频率
d=(C/f0)/2;%两个相邻阵元距离，不能大于载波的半波长
Radius=2*d;%均匀圆阵半径
Bearings=[40];%入射角方位角
Pitchings=[60];%入射信号俯仰角


AmpCoeff=[100];%阵元放大系数
f=[500];%基带信号频率
fs=4*max(f)+1;%采样频率


AmpCoeff=repmat(AmpCoeff,[1 N]);
S=AmpCoeff.*sin(2*pi*f*[1:N]/fs);
Noise=randn(M,N);
Bearings=Bearings*pi/180;
Pitchings=Pitchings*pi/180;


TimeDelay=d*[0:M-1].'*sin(Bearings)/C;%线阵的延迟

A=exp(-1*sqrt(-1)*2*pi*f0*TimeDelay);

X=A*S+50*Noise;

Rxx=(X*X')/N;  %correlation matrix

[EigenVectors,EigenValues]=eig(Rxx);

Lemda=diag(EigenValues);  %计算矩阵特征值
[SortedLemda,Index]=sort(Lemda);
Index=flipud(Index);%将特征值降序排列


NoiseSubspace(1:M,1:M-P)=EigenVectors(1:M,Index(P+1:M));


Delta=[-90:1:90];
for  k=1:length(Delta)
    sum_spec(k)=0;
end

for i=1:1:180
    a(1,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*0/8-Delta*(pi/180))*cos(i*(pi/180))/C);
    a(2,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*1/8-Delta*(pi/180))*cos(i*(pi/180))/C);
    a(3,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*2/8-Delta*(pi/180))*cos(i*(pi/180))/C);
    a(4,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*3/8-Delta*(pi/180))*cos(i*(pi/180))/C);
    a(5,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*4/8-Delta*(pi/180))*cos(i*(pi/180))/C);
    % a(6,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*5/8-Delta*(pi/180))*cos(i*(pi/180))/C);
%     a(7,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*6/8-Delta*(pi/180))*cos(i*(pi/180))/C);
%     a(8,:)=exp(-1*sqrt(-1)*2*pi*f0*Radius*cos(2*pi*7/8-Delta*(pi/180))*cos(i*(pi/180))/C);


for k=1:length(Delta)
    MUSIC_Spec(k)=1/(a(:,k)'*NoiseSubspace*NoiseSubspace'*a(:,k));
end

for m=1:length(Delta)
    sum_spec(m)=sum_spec(m)+MUSIC_Spec(m);
end

for k=1:length(Delta)
    single_spec(i,k)=MUSIC_Spec(k);
end


end

figure(1);
Delta1=[-45:1:135];
plot(Delta1,10*log10(sum_spec)); %NonNormative  Spatial Spectrum
figure(2);
mesh(single_spec);

r=1;
for k=2:length(Delta)-1
    if(sum_spec(k)>=sum_spec(k-1)&&sum_spec(k)>=sum_spec(k+1))
        LocalPeak(r)=sum_spec(k);
        LocalPeak_Angle(r)=Delta(k);
        r=r+1;
    end
end

[LocalPeak,Index_Angle]=sort(LocalPeak);
p=1:P
EstimatedBearing(p)=LocalPeak_Angle(Index_Angle(r-p))+45