clc
clear all
close all
signalnum=3;%信号源
SNR1=0;
SNR2=0;%选择三个信号源,定义其信噪比。
SNR3=0;%信噪比
jay=sqrt(-1);%定义虚数分量的定义,就是j
angle1=30;
angle2=20;
angle3=-10;%到达角度,三个,可以自己根据需要 设置。
angle1=angle1*pi/180;
angle2=angle2*pi/180;
angle3=angle3*pi/180;%角度转换,将角度转换成弧度制。
SA1=sqrt(2*10^(SNR1/10));
SA2=sqrt(2*10^(SNR2/10));
SA3=sqrt(2*10^(SNR3/10));%信号的振幅,主要控制信号的能量
c=1500;%定义的波速,这里按照水声的速度先定义吧,无所谓的,雷达波速度不影响这个结论
wavelength=1;
M = 8;%阵元数
fa=c/wavelength;%定义的波频率
d=wavelength/2;%阵元距
len=2;%采样时长
fs=3000;
f1=700;
f2=600;
f3=500;%三个达到波的频率
t=[0:1/fs:len];%采样拍数
ss1=sin(2*pi*f1*t);
s1=SA1*ss1;
ss2=sin(2*pi*f2*t);%这个主要写的是信号的采样,其实就是经过数字信号转换以后的程序。
s2=SA2*ss2;
ss3=sin(2*pi*f3*t);
s3=SA3*ss3;
signal=[s1;s2;s3];%信号矩阵
%阵列流型
ta1=([0:d:(M-1)*d]*sin(angle1))/c;
ta2=([0:d:(M-1)*d]*sin(angle2))/c;
ta3=([0:d:(M-1)*d]*sin(angle3))/c;
a1=exp(-jay*2*pi*fa*ta1).';
a2=exp(-jay*2*pi*fa*ta2).';
a3=exp(-jay*2*pi*fa*ta3).';
direction=[a1 a2 a3];%阵列方向矢量
SRr=direction*signal+10*randn(M,length(s1));%要加上噪声,在做程序时候,要仔细观察噪声的分量
Rr=SRr*SRr'/length(t);
[u v]=eig(Rr)%特征值分解
for i=1:M-signalnum
e1(:,i)=u(:,i);
end
Enr=e1;%噪声空间
theta=[-90:0.1:90];
for i=1:length(theta)
angle4=theta(i)*pi/180;
for j=1:M
ae(j)=exp(-jay*2*pi*fa*(j-1)*d*sin(angle4)/c);
end
a=ae.';%
Prmusic(i)=1/(a'*Enr*Enr'*a);%music判断函数,这个要看空间谱估计那本书。
end
PMUSIC=20*log10(abs(Prmusic));
figure
plot(theta,PMUSIC,'-g');%画图程序,简单看一下。
grid on
title('music算法');