music.m

经行direction of arrival(DOA)估计的MUSIC算法。

%--------------------------------------------

%相干信号DOA估计MD（矩阵分解）算法 
%适用于均匀线阵
%修正后的music算法
%程序分为信源、全局变量定义、信号产生、预处理算法、DF算法、做图等几部分
%盒子 于200803~200805期间研学，无版权
%--------------------------------------------

clc;
%1：信源
%--------------------------------------------
Pd=2000;
Fd=1;
Fs=4*Fd;
R=0.5;
Delay=5;
M=4;
x1=randint(Pd,1,M); %生成一个随机数为2000*1的矩阵，输出范围为0~3
x2=randint(Pd,1,M);
y1=modmap(x1,Fd,Fs,'qask',M);%数字信号调制模拟信号，modmap可用scatterplot代替
y2=modmap(x2,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);
s1=amodce(rcv_a1,10,'qam');%输出复包络，amodce可用ammod来代替,s1为8040*1
% s2=amodce(rcv_a2,10,'qam');
s2=s1; %假设信号s1和s2相关
save sig3 s1 s2 %把变量s1,s2保存为sig3.mat文件
%--------------------------------------------

%2：全局变量定义
%--------------------------------------------
% clear
m=8; %阵元数
angle1=-50;%信号角度
angle2=60;
th=[angle1;angle2];%th为2*1矩阵
nn=1024;%采样数
SN1=10;%信号的SNR
SN2=10;
SN3=12;
sn=[SN1;SN2];%sn为2*1矩阵
degrad=pi/180;
%--------------------------------------------

%3：生成离散信号
%--------------------------------------------
load sig3 %把sig3.mat文件中的s1,s2变量装入内存
tt=1:nn;
S=[s1(tt).';s2(tt).'];%S为2*1024矩阵
nr=randn(m,nn);%nr为8*1024矩阵,为正态分布随机阵
ni=randn(m,nn);
U=nr+j*ni; %构造噪声源，U为8*1024复数矩阵
Ps=S*S'/nn;%信号的协方差矩阵，Ps为2*2矩阵
ps=diag(Ps);%取主对角线上值，ps为2*1矩阵
refp=2*10.^(sn/10);
tmp=sqrt(refp./ps);%tmp为2*1矩阵
S2=diag(tmp)*S;%S2为2*1024矩阵
%--------------------------------------------

%4：协方差矩阵及前后向平滑
%--------------------------------------------
tmp2=[0:m-1]'; %阵元位置
A=exp(-i*pi*tmp2*sin(th'*degrad)); %A为8*2矩阵，方向矩阵由公式得出
X=A*S2+U;%X为8*1024矩阵，1024次采集后得到的数据，阵列接受信号
Rxx=X*X'/nn;%Rxx为8*8矩阵，阵列接受信号协方差

a=eye(8);
b=a(:,8:-1:1);%构造交换矩阵
RXX=(Rxx+b*Rxx.'*b)/2;
%--------------------------------------------

%5：奇异值分解
%--------------------------------------------
[Q,SS,W]=svd(RXX); 
Vu=Q(:,3:8 ); %小特征值对应的特征向量，即噪声子空间
%--------------------------------------------

%6：做图
%--------------------------------------------
th2=[-90:90]';
tmp2=[0:m-1]';
a2=-i*tmp2*pi*sin(th2'*degrad);
A2=exp(a2);
num=diag(A2'*A2);
Ena=Vu'*A2;
den=diag(Ena'*Ena);
doa=num./den;
semilogy(th2,abs(doa));
axis([-90 90 0.1 1e5]);
grid on;
hold on;
%--------------------------------------------