music5.m

music 算法仿真 三维仿真图像

%--------------------------------------------
%适用于均匀圆阵
%盒子 于200803~200805期间研学，无版权
%--------------------------------------------
clc;
clear all;
%--------------------------------------------

%1：信源
%--------------------------------------------
Pd=2000;
Fd=1;
Fs=4*Fd;
R=0.5;
Delay=5;
M=4;
x1=randint(Pd,1,M); 
x2=randint(Pd,1,M);
x3=randint(Pd,1,M);
y1=modmap(x1,Fd,Fs,'qask',M);
y2=modmap(x2,Fd,Fs,'qask',M);
y3=modmap(x3,Fd,Fs,'qask',M);
[rcv_a1,ti]=rcosflt(y1,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a2,ti]=rcosflt(y2,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a3,ti]=rcosflt(y3,Fd,Fs,'fir/sqrt/Fs',R,Delay);
s1=amodce(rcv_a1,10,'qam'); 
s2=amodce(rcv_a2,10,'qam');
s3=amodce(rcv_a3,10,'qam');
save sig3 s1 s2 s3;
%--------------------------------------------

%2：全局变量定义
%--------------------------------------------
m=8; % 阵元数
p=3; % 信号数
%俯仰角为入射方向与圆阵轴线之间的夹角,方位角为x轴到入射线在圆阵上投影的夹角
st1=10; % 俯仰角
dt1=10; % 方位角
st2=50;
dt2=30;
st3=70;
dt3=30;

st=[st1;st2;st3];
dt=[dt1;dt2;dt3];

nn=4096; 
SN1=20; 
SN2=20;
SN3=20;
sn=[SN1;SN2;SN3];
degrad=pi/180;
%--------------------------------------------

%3：生成离散信号
%--------------------------------------------
load sig3
tt=1:nn;
S=[s1(tt),s2(tt),s3(tt)].';
noise=randn(m,nn)+j*randn(m,nn); 
Ps=S*S'/nn;
ps=diag(Ps);
refp=2*10.^(sn/10);
tmp=sqrt(refp./ps);
S2=diag(tmp)*S;
%--------------------------------------------

%4：圆半径为半波长
%--------------------------------------------
A=zeros(m,p);
k=[0:m-1]';
for t=1:p
A(:,t)=exp(j*pi*(sin(st(t)*degrad)*cos(dt(t)*degrad)*cos(2*pi*k/m)+sin(st(t)*degrad)*sin(dt(t)*degrad)*sin(2*pi*k/m)));
end
X=A*S2+noise;
%--------------------------------------------

%5：协方差矩阵，奇异值分解
%--------------------------------------------
Rxx=X*X'/nn;
[U,s,v]=svd(Rxx);
Vs=U(:,1:p);
Vu=U(:,p+1:m);
%--------------------------------------------

%6：做谱图
%--------------------------------------------
for st1=0:90
for dt1=0:90
k=[0:m-1]';
AA=exp(j*pi*(sin(st1*degrad)*cos(dt1*degrad)*cos(2*pi*k/8)+sin(st1*degrad)*sin(dt1*degrad)*sin(2*pi*k/8)));
WW=AA'*Vu*Vu'*AA;
Pmusic(st1+1,dt1+1)=abs(8./WW);%角谱
end 
end 
sst=0:90;
ddt=0:90;
mesh(sst,ddt,Pmusic);
title('MUSIC空间谱');
xlabel('方位角'); 
ylabel('俯仰角');
zlabel('角谱');
grid on;