www.pudn.com > music.rar > music.m

%-------------------------------------------- 
 
%相干信号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; 
%--------------------------------------------