music1.m

《统计与自适应信号处理》中的程序

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

%一段估计DOA的程序
%适用于均匀线阵
%没有修正的MUSIC算法
%盒子 于200803~200805期间研学，无版权
%--------------------------------------------

%1:构造信号源，对QASK数字调制信号做三个波
%--------------------------------------------
Pd=2000;
Fd= 1;
Fs= 4*Fd;
R= 0.5;
Delay= 5;
No=1;
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);
%y1=scatterplot(x1,Fd, Fs,'qask',M);
%y2=scatterplot(x2,Fd, Fs,'qask',M);
%y3=scatterplot(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：天线接收信号
%--------------------------------------------
%clear
i= sqrt(-1);
j= i;
m=8;%天线阵元数
p=2;%信号个数
angle1=-30;
angle2=50;
angle3=45;%信号来向
%th=[angle1;angle2;angle3];
th=[angle1;angle2];
nn=1024;%采样数
SN1=12; SN2=10; SN3=9;%信号的SNR
%sn=[SN1;SN2;SN3];
sn=[SN1;SN2];
degrad=pi/180;
%--------------------------------------------

%3：构造噪声和白噪声
%--------------------------------------------
%load sig3
tt=1:nn;
%S=[s1(tt).';s2(tt).';s3(tt).'];
S=[s1(tt).';s2(tt).'];
nr=randn(m,nn);
ni=randn(m,nn);
u=nr+j *ni;%构造噪声源,服从（0，2）分布

Ps=S*S'/nn;
ps=diag(Ps);
refp=2*10.^(sn/10);
tmp=sqrt(refp./ps);
S2=diag(tmp)* S;
%--------------------------------------------

%4：计算协方差矩正解特征向量
%--------------------------------------------
tmp =-i*pi*sin(th'*degrad);
tmp2=[0:m-1 ]';
a2=tmp2*tmp;
A=exp(a2);%天线阵子

X=A*S2+u;%阵列接收信号
Rxx=X*X'/nn;%接收信号协方差
[U,S,V]=svd(Rxx);% （SVD）奇异值分解
Vs=U(:,1:p);
Vu=U(:,p+1:m);
%--------------------------------------------

%5：空间谱函数
%--------------------------------------------
th2=[-90:90]';
tmp=-i* pi* sin(th2'* degrad);
tmp2=[0:m-1 ]';
a2=tmp2*tmp;
A2=exp(a2);%A2为8*181矩阵，使角度变化180次

num=diag(A2'*A2);%分子为181*1的矩阵

Ena=Vu'* A2;%信号向量与噪声向量乘积
den=diag(Ena'* Ena);%分母也为181*1的矩阵
doa=(num./den);%空间谱函数
%--------------------------------------------

%6：作出谱图
%--------------------------------------------
%plot(th2, doa);
semilogy(th2, abs(doa));
title('MUSIC Spectrum');
xlabel('Angle(deg)');%横轴
ylabel('Spectrum');%纵轴
axis([-90 90 0.1 1e5]);
grid;