music.m

多波束系统中高分辨算法－－MUSIC算法的MATLAB实现的源程序

%二：空间有三个独立源，用Music,R-Music,进行DOA估计
%　　能分辨的角度和不能分辨的各两个结果进行比较。
clc
clear
%--------------------基本参数----------------------------------
M=8;                                                %阵元数
N=256;                                              %快拍数
SNR=25;                                             %信噪比
f=1500;                                             %信号频率
fs=5*f;                                             %采样频率                             
rsj=[-6 1 3];                                     %信号源入射角
Q1=-25;                                             %扫描扇区
Q2=15;
xys=3;                                              %信源数
z=0:M-1;
n=1:N;
%-------------------产生三个独立信号----------------------------
s1=exp(-j*2*pi*f/fs*n);
s2=2*exp(-j*2*pi*2*f/fs*n);
s3=5*exp(-j*2*pi*3*f/fs*n);
s=[s1;s2;s3];                                  %构造三个独立信号数据s(3x256)
A=exp(-j*pi*z'*sin(rsj*pi/180));               %导向矢量(8x3)
S=(10^(SNR/20))*A*s;                           %信号空间(8x256)
noise=hilbert(randn(M,N))/sqrt(2);             %噪声空间(8x256)
X=S+noise;                                     %接受数据空间
R=1/N*X*X';                                    %数据协方差矩阵
%
[C,D]=eig(R);                 %特征分解，C为特征向量，D为特征值
[D1,pos1]=max(D);             %使D变为单行矩阵
for l=1:3                     %去信号,1:3--指去三个信号
    [D2,pos2]=max(D1);
     D1(:,pos2)=[];
     C(:,pos2)=[];
end
C;                            %C为噪声空间Un（4x1)
%－－－－－－－－－－－－谱峰搜索－－－－－－－－－－－－－－
k=length(Q1:0.1:Q2);
for n=1:k
    i=1:M;                
    a=exp(-j*pi*(i-1)'*sin((Q1+(n-1)*0.1)*pi/180));
    Pm(n)=1/(a'*C*C'*a);
end
Pu=10*log10(abs(Pm));
%－－－－－－－－－－－－绘制图形－－－－－－－－－－－－－
figure
plot(Q1:0.1:Q2,abs(Pu));xlabel('入射角(度)');ylabel('空间方位谱(dB)')
title('MUSIC')
grid on
%－－－－－－－－－－－－求根MUSIC－－－－－－－－－－－－－－
syms z;
p=[1;z;z^2;z^3;z^4;z^5;z^6;z^7];
pt=[1,z^(-1),z^(-2),z^(-3),z^(-4),z^(-5),z^(-6),z^(-7)];
f=z^7*pt*C*C'*p;                   %定义多项式
xs=sym2poly(f);                    %将符号多项式转化为系数多项式
gen=roots(xs);                     %求根
arg=angle(gen);                    %求复角主值
b=[9,11,13];
Q=-arg(b);
DOA=asin(1/pi*Q')*180/pi            %求信号源角度,公式φ=(2π|Δ|sina）/λ