基于MATLAB的MUSIC算法程序

admin

%MUSIC spectrum
clear all;
clc;
%参数赋值
c=3.0e+8;                                             %光速
lmda=0.24;                                            %L波段
fc=c/lmda;                                            %载频
N=16;                                                 %阵元数
DOA_sigA=0;                                         %信号源1来波方向与阵法线方向夹角
DOA_sigB=20;                                          %信号源2来波方向与阵法线方向夹角
d=lmda/2;                                             %阵元间距
fs=1.0e+3;                                            %采样率
chirp_rateA=20;
chirp_rateB=10;                                       %调频率
P=2;                                                  %信源数
%计算信源复包络
t=[0:1/fs:1]';                                        %采样时间
s1=exp(j*pi*chirp_rateA*t.^2);                                    
s2=exp(j*pi*chirp_rateB*t.^2);
% carw=exp(j*2*pi*fc*t);                              %载波
% sc1=s1.*carw;                                       %信源复包络                                       
% sc2=s2.*carw;

%产生方向矢量1(steering vector)
sv_sigA=[];
for i=1:N
    sv_sigA=[sv_sigA;exp(j*2*pi*d*(i-1)*sin(DOA_sigA*pi/180)/lmda)];
end
T1=[];
for i=1:N
    Ei=i;
    sitai=i+1;
    T1=[T1;Ei*exp(j*sitai)];
end
sv_sigA1=diag(T1)*sv_sigA;
%产生方向矢量2(steering vector)
sv_sigB=[];
for i=1:N
    sv_sigB=[sv_sigB;exp(j*2*pi*d*(i-1)*sin(DOA_sigB*pi/180)/lmda)];
end
T2=[];
for i=1:N
    Ei=i-0.5;
    sitai=i;
    T2=[T2;Ei*exp(j*sitai)];
end
sv_sigB1=diag(T2)*sv_sigB;
%N元等距线阵接收信号矢量
x=[sv_sigA1,sv_sigB1]*[transpose(s1);transpose(s2)]+randn(N,length(t));
%阵列输出协方差矩阵
Rx=zeros(N,N);
for i=1:length(t)
    Rx=Rx+x(:,i)*x(:,i)';
end
Rx=Rx/length(t);
%对R做特征值分解
[V,D]=eig(Rx);%按照特征值从小到大排列在D的对角线上。
%get the noise subspace.
noise_subspace=V(:,1:N-P);
%搜寻谱峰
MUSICSpec=[];
for sita=-pi/2:1/100:pi/2;
    sv_sita=[];
    for i=1:N
        sv_sita=[sv_sita;exp(j*2*pi*d*(i-1)*sin(sita)/lmda)];
    end
   
    s=0;
    for k=1:N-P;
        s=abs((noise_subspace(:,k))'*sv_sita)^2+s;
    end
    MUSICSpec=[MUSICSpec;1/s];
end
%搜寻谱峰
MUSICSpec_en=[];
for sita=-pi/2:1/100:pi/2;
    sv_sita=[];
    for i=1:N
        sv_sita=[sv_sita;exp(j*2*pi*d*(i-1)*sin(sita)/lmda)];
    end
   
    s=0;
    for k=1:N-P;
        s=abs((noise_subspace(:,k))'*sv_sita)^2/D(k,k)+s;%按照特征值的大小进行加权处理
    end
    MUSICSpec_en=[MUSICSpec_en;1/s];
end
%搜寻谱峰(maximal likelihood method)
%估计噪声功率.
deltn=0;
for m=1:N-P
    delt=abs(D(m+1,m+1)-D(m,m))/D(m,m);
    deltn=deltn+D(m,m);
end
deltn=deltn/(N-P);
clear delt;
%maximal likelihood signal subspace
ML_sigsubspace=zeros(N,N);
for m=N-P+1:N
    ML_sigsubspace=ML_sigsubspace+deltn*(D(m,m)/(D(m,m)-deltn)^2)*V(:,m)*V(:,m)';
end
MUSICSpec_ML=[];
for sita=-pi/2:1/100:pi/2;
    sv_sita=[];
    for i=1:N
        sv_sita=[sv_sita;exp(j*2*pi*d*(i-1)*sin(sita)/lmda)];%sv_sita 以列向量形式表示。
    end
    s=(sv_sita'*noise_subspace*noise_subspace'*sv_sita)/(sv_sita'*ML_sigsubspace*sv_sita);
    MUSICSpec_ML=[MUSICSpec_ML;1/s];
end
%绘图
sita=-pi/2:1/100:pi/2;
sita_deg=sita*180/pi;
MUSICSpec_log=10*log10(abs(MUSICSpec/max(MUSICSpec)));
MUSICSpec_en_log=10*log10(abs(MUSICSpec_en/max(MUSICSpec_en)));
MUSICSpec_ML_log=10*log10(abs(MUSICSpec_ML/max(MUSICSpec_ML)));

figure(1)
plot(sita_deg,MUSICSpec_log)
grid
axis([-90,90,-40,0])
title('DOA估计图-MUSICSpec')
xlabel('方位角（度）')
ylabel('输出功率（dB）')
hold on
stem(DOA_sigA,-10,'*')
stem(DOA_sigB,-10,':')

figure(2)
plot(sita_deg,MUSICSpec_en_log)
grid
axis([-90,90,-40,0])
title('DOA估计图-MUSICSpeEn')
xlabel('方位角（度）')
ylabel('输出功率（dB）')
hold on
stem(DOA_sigA,-10,'*')
stem(DOA_sigB,-10,':')

figure(3)
plot(sita_deg,MUSICSpec_ML_log)
grid
axis([-90,90,-60,0])
title('DOA估计图-MUSICSpeML')
xlabel('方位角（度）')
ylabel('输出功率（dB）')
hold on
stem(DOA_sigA,-10,'*')
stem(DOA_sigB,-10,':')