doaestimate2.m

通信信号处理中的七阵元四信号元DOA估计

%music.m 
clear 
clear all 
%构造信号源 
n=-pi:2.0/2000*pi:pi; %范围 
fs=2000;
f1=0.015;
f2=0.05;
f3=0.02;
f4=0.035;
%s1=100*randn(1)*cos(100*randn(1)*sqrt(2.0)*pi*n+pi*randn(1)); %构造四个余弦信号 
%s2=100*randn(1)*cos(100*randn(1)*sqrt(2.0)*pi*n+pi*randn(1)); 
%s3=100*randn(1)*cos(100*randn(1)*sqrt(2.0)*pi*n+pi*randn(1)); 
%s4=100*randn(1)*cos(100*randn(1)*sqrt(2.0)*pi*n+pi*randn(1)); 

s1=cos(100*f1*sqrt(2.0)*pi*n); %构造四个余弦信号 
s2=cos(100*f2*sqrt(2.0)*pi*n); 
s3=cos(100*f3*sqrt(2.0)*pi*n); 
s4=cos(100*f4*sqrt(2.0)*pi*n); 
save sig4 s1 s2 s3 s4

%天线阵元及接收信号 
i=sqrt(-1); 
j=i; 
m=7;%天线阵元数 
p=3; 
angle1=45;%从-60 0 60入射的信号源 
angle2=60; 
angle3=30; 
angle4=-45;
th=[angle1;angle2;angle3;angle4]; 
nn=2000;%采样点数 
SN1=10;%信噪比 
SN2=10; 
SN3=10;
SN4=10;
sn=[SN1;SN2;SN3;SN4]; 
degrad=pi/180; 

load sig4 
tt=1:nn; 
S=[s1(tt);s2(tt);s3(tt);s4(tt)]; 
nr=randn(m,nn); 
ni=randn(m,nn); 
U=nr+j*ni; %建立高斯噪声 

Ps=S*S'/nn; 
ps=diag(Ps); 
refp=2*10.^(sn/10); 
tmp=sqrt(refp./ps); %refp,tmp,[;;] 
S2=diag(tmp)*S;%diag() then3*3 ,diag^ 

%计算协方差矩阵，特征值分解 
tmp22=-i*pi*sin(th'*degrad); %相位差 
tmp2=[0:m-1]'; %信号到达阵元的相位差组成的向量 
a2=tmp2*tmp22; 
A=exp(a2); %m*p维方向矩阵 
X=A*S2+U; %第k个阵元上的接收信号
%X=A*S2;
Rxx=X*X'/nn; %接收信号的协方差 

%计算空间谱函数 
th2=[-90:90]'; %线性信号角度范围 
tmp33=-i*pi*sin(th2'*degrad); %[-90 90]范围内的相位差组成的向量 
tmp2=[0:m-1]'; 
a2=tmp2*tmp33; 
A2=exp(a2); %m*p维方向矩阵 

[V,D]=eig(Rxx); 
G=V(:,1:m-p); %信号与噪声分离 
for n=1:length(th2) 
doa(n)=1.0/(A2(:,n)'*G*G'*A2(:,n)); 
end 

%做出空间谱图 
semilogy(th2,doa); %y轴对数坐标曲线 
title('MUSIC Spectyum'); %图表标题，x轴y轴标签及范围 
xlabel('Angle(deg)'); 
ylabel('Spectrum'); 
axis([-90 90 0.1 1e3]); 
grid;