beam_music.asv

阵元间隔为半波长的均匀分布16元线阵

% %阵元间隔为半波长的均匀分布16元线阵，预成指向士arcsin（(2i-1)/16）(i=l，2，…，16)的16个均匀加权常规波束,远场有5个互
% 不相关的目标源发射信号，每个源相对于基阵的方位为-40，-42，-5，30和33，且每个信号源到达基阵的信噪比相同，均为5DB，快拍数
%为1000.
clear all
close all
C=340;                                      %%声速
d=2;                                        %%声源数
num=16;                                     %% 阵元数
k=1000;                                     %% 快拍数
Bearings=[-40 -42];                         %入射信号方位角
D=0.1;                                      %阵元间距
fc=1700;                                    %信号频率
fs=2*fc;                                    %采样频率
snr=5;                                      %信噪比
e_position=[0:num-1]';
s_position=[0:d-1]';

%%波束形成
as=-13/16:1/8:-9/16;                        %波束指向角
as=asin(as);
vs=exp(j*2*pi*fc*e_position*D*sin(as)/C);
w=1/num*vs';                                %波束形成矩阵

%%入射信号
aa=90*[-1:0.002:1];                         %从-90到90采样1001个点
s=exp(j*pi*sin(aa/180*pi));                 
S=repmat(s,[d 1]);                          %声源信号

%%方向向量
Bearings=Bearings*pi/180;
TimeDelay=D*[0:num-1].'*sin(Bearings)/C;    %线阵的延迟
A=exp(sqrt(-1)*2*pi*fc*TimeDelay);          %方向向量

%%基阵接收信号
x=A*S;                                      %阵列输出信号
X=awgn(x,snr);                              %加噪声

%波束输出
Y=w*X;                                      %方位估计

%%波束域MUSIC算法
B=3;                                        %波束数
Ryy=(Y*Y')/k;                             %协方差矩阵
[EigenVectors,EigenValues]=eig(Ryy);
Lemda=diag(EigenValues);                    %计算矩阵特征值
[SortedLemda,Index]=sort(Lemda);
Index=flipud(Index);                        %将特征值降序排列
NoiseSubspace(1:B,1:B-d)=EigenVectors(1:B,Index(d+1:B));%把噪声空间分成两部分，分别为加入信号空间的一部分和噪声

%%进行搜索
for i=-90:1:90
    
a=exp(sqrt(-1)*2*pi*fc*D*[0:num-1]'*sin(i*pi/180)/C);

    MUSIC_Spec(i+91)=a'*w'*w*a/(a'*w'*NoiseSubspace*(w'*NoiseSubspace)'*a);

end

figure(1);
Delta1=[-90:1:90];
plot(Delta1,10*log10(MUSIC_Spec)); %NonNormative  Spatial Spectrum




%画图
% length=55*[-1:2/15:1];
% figure(1)
% plot(20*log10(abs(Y)/max(abs(Y))));
% axis([-90 90 -20 0]);
% hold on
% % 


   
% e_position1=[0:num-1]';
% aa1=90*[-1:0.002:1];                                            %% scan vector
% vs=exp(j*pi*e_position1*sin(as));
% s=exp(j*pi*e_position1*sin(aa1/180*pi));  %%  入射信号
% 
% w1=1/num*vs';     %%加权矩阵
% 
% figure(1),plot(aa,20*log10(abs(w1*s))),xlabel('角度（度）'),ylabel('指向性(DB)');
% axis([-90 90 -40 0]);
%