f1.m

一个关于MUSIC算法的源程序

%构造信号源
n=-pi:0.01:pi;                           %范围    
x1=100*cos(50*sqrt(2.0)*pi*n+pi/3);      %构造三个余弦信号
x2=100*cos(50*sqrt(2.0)*pi*n+pi/3);
x3=100*cos(50*sqrt(2.0)*pi*n+pi/3);
save shy3 x1 x2 x3                        %保存信号源

%天线阵元及接收信号
clear
close all
i=sqrt(-1);    %复数单位
j=i;
m=8;      %天线阵元数
p=3;      %信号个数
angle1=45;angle2=60;angle3=75;    %信号来向
angle=[angle1;angle2;angle3];     
cys=600;   %采样数 
SN1=20;SN2=10;SN3=12;  %信噪比
SN=[SN1;SN2;SN3];
degrad=pi/180;  %将角度制化为弧度制    

%信号和噪声
load shy3    %载入信号源
TT=1:cys;    %采样数矩阵
S=[x1(TT).';x2(TT).';x3(TT).']; 
nr=randn(m,cys);  %建立高斯随机序列
ni=randn(m,cys);
E=nr+ni*j;        %建立高斯噪声
Ps=S*(S/cys)';    %方向向量
Ps1=diag(Ps);
refp=2*10.^(SN./10);
tmp=sqrt(refp./Ps1);
S1=diag(tmp)*S;

%计算协方差矩阵，特征值分解
tmp=-i*pi*sin(angle'.*degrad);  %相位差
tmp1=[0:m-1]';        %信号到达阵元的相位差组成的向量
A=tmp1*tmp;
A1=exp(A);            %m*p维方向矩阵
X=(A*S1)+E;           %第k个阵元上的接收信号
Rxx=X*(X/cys)';       %接收信号的协方差
[E,S,V]=svd(Rxx);     %计算特征值
Vs=E(:,1:p);          
Vu=E(:,p+1:m);

%计算空间谱函数
th2=[-90:90]';         %线性信号角度范围
tmp2=-i*pi*sin(th2'.*degrad);  %[-90 90]范围内的相位差组成的向量
tmp3=tmp1;
A2=tmp2*tmp3;       
A3=exp(A2);            %m*p维方向矩阵
num=diag(A3'*A3);      
Ena=Vu'*A3;
den=diag(Ena'*Ena);     
doa=num./den;          %MUSCI估计值

%做出空间谱图
semilogy(th2,doa);     %y轴对数坐标曲线
title('MUSIC算法');    %图表标题，x轴y轴标签及范围
xlable('角度（弧度制）');
ylable('值');
axis([-90 90 0.1 1e4]);