dcxhflmusic.txt

music算法进行雷达信号方向角的谱估计

%MUSIC  方法
clc
m=8;                                                                                                                      %等距线阵阵元的个数
p=3;                                                                                                                       %远场信号的个数
M=5;                                                                                                                      %快拍数据的个数
p_x=1024;                                                                                                            %每个远场信号的长度
l=1.8;                                                                                                                     %每个远场信号的波长
d=0.5*l;                                                                                                                 %相邻阵元之间的距离
w=[pi/6,pi/10,pi/8]';                                                                                              %三个远场信号的频率
theta1=pi/4;
theta2=0;
theta3=-pi/4;                                                                                                         %三个远场信号的波达方向角
n=[0:m-1]';                                                                                                             %生成导向矢量矩阵的行数
thetap=[theta1,theta2,theta3];                                                                             %生成三个远场信号波达方向角矩阵，以便生成导向矢量矩阵
Btheta=exp(-j*2*pi*n*d*sin(thetap)/l);                                                                 %生成m x p的导向矢量矩阵，阵元数决定改矩阵的行数，远场信号数决定列数
s=2*exp(j*w*[0:p_x-1]);                                                                                        %生成相互独立的仿真信号，每行代表一个远场信号


x1=Btheta*s+randn(8,p_x)+j*randn(8,p_x);                                                       %生成加入随机噪声的远场信号1
x2=Btheta*s+randn(8,p_x)+j*randn(8,p_x);                                                       %生成加入随机噪声的远场信号2
x3=Btheta*s+randn(8,p_x)+j*randn(8,p_x);                                                       %生成加入随机噪声的远场信号3
x4=Btheta*s+randn(8,p_x)+j*randn(8,p_x);                                                       %生成加入随机噪声的远场信号4
x5=Btheta*s+randn(8,p_x)+j*randn(8,p_x);                                                       %生成加入随机噪声的远场信号5


%计算输入信号的自相关矩阵
 R=(x1*x1'+x2*x2'+x3*x3'+x4*x4'+x5*x5')/M;


%对R进行特征值分解
[U,D,V]=svd(R);
G=U(:,p+1:m);
C=G*G';                                                                                                                


%进行估计远场信号的波达方向角
theta=-90:0.1:90;                                                                                                 
for   k=1:length(theta)
       Atheta=exp(-j*2*pi*n*d*sin(theta(k)/180*pi)/l);                                          
       P=Atheta'*C*Atheta;                                                                                    
       Pmusic(k)=abs(1/P);                                                                                    
end
Pmusic=20*log10(Pmusic);                                                                              

%画出空间谱的图形，找到它的p个峰值，它们就是待估计值，即外场信号的波达方向角
figure(1);
plot(theta,Pmusic);
axis([-100,100,-20,80]);
title('多重信号分类MUSIC方法');
grid on                                                                                                                  %显示分格线
xlabel('\theta')
ylabel('空间谱幅值的dB值')


%画方向性图

%当theta1=pi/4时的天线方向图
Btheta1=exp(-j*2*pi*n*d*sin(theta1)/l);                                                              %取权向量为Btheta1/m
theta=-90:0.1:90;                                                                                                 %在波达方向角的变化范围内画图
for   k=1:length(theta)
       Atheta=exp(-j*2*pi*n*d*sin(theta(k)/180*pi)/l);                                           
       P1=Btheta1'*Atheta/m;                                                                                 
       Ptheta1(k)=abs(P1);                                                                                      
end
Ptheta1=20*log10(Ptheta1);                                                                               
figure(2);
% subplot(3,1,1)
plot(theta,Ptheta1);
axis([-100,100,-90,10]);
title('\theta=\pi/4时的方向图');
grid on
xlabel('\theta')
ylabel('方向图幅值的dB值')

%当theta2=0时的天线方向图
Btheta2=exp(-j*2*pi*n*d*sin(theta2)/l);                                                              %取权向量为Btheta2/m
theta=-90:0.1:90;                                                                                                 %在波达方向角的变化范围内画图
for   k=1:length(theta)
       Atheta=exp(-j*2*pi*n*d*sin(theta(k)/180*pi)/l);                                           
       P2=Btheta2'*Atheta/m;                                                                                 
       Ptheta2(k)=abs(P2);                                                                                      
end
Ptheta2=20*log10(Ptheta2);                                                                               
figure(3);
% subplot(3,1,2)
plot(theta,Ptheta2);
axis([-100,100,-350,50]);
title('\theta=0时的方向图');
grid on
xlabel('\theta')
ylabel('方向图幅值的dB值')

%当theta3=-pi/4时的天线方向图
Btheta3=exp(-j*2*pi*n*d*sin(theta3)/l);                                                              %取权向量为Btheta3/m
theta=-90:0.1:90;                                                                                                 %在波达方向角的变化范围内画图
for   k=1:length(theta)
       Atheta=exp(-j*2*pi*n*d*sin(theta(k)/180*pi)/l);                                           
       P3=Btheta3'*Atheta/m;                                                                                 
       Ptheta3(k)=abs(P3);                                                                                      
end
Ptheta3=20*log10(Ptheta3);                                                                               
figure(4);
% subplot(3,1,3)
plot(theta,Ptheta3);
axis([-100,100,-90,10]);
title('\theta=-\pi/4时的方向图');
grid on
xlabel('\theta')
ylabel('方向图幅值的dB值')

%把方向性图画在同一坐标系中
figure(5);
plot(theta,Ptheta1,'b-');
axis([-100,100,-350,50]);
hold on
plot(theta,Ptheta2,'r:');
hold on
plot(theta,Ptheta3,'m-.');
title('\theta=-\pi/4,0,\pi/4时的方向图');
grid on
xlabel('\theta')
ylabel('方向图幅值的dB值')
legend('Ptheta1','Ptheta2','Ptheta3',4);