文件頭:
%構(gòu)造信號源

n=-pi:0.01:pi;                           %范圍    

x1=100*cos(50*sqrt(2.0)*pi*n+pi/3);      %構(gòu)造三個余弦信號

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);    %復(fù)數(shù)單位

j=i;

m=8;      %天線陣元數(shù)

p=3;      %信號個數(shù)

angle1=45;angle2=60;angle3=75;    %信號來向

angle=[angle1;angle2;angle3];     

cys=600;   %采樣數(shù) 

SN1=20;SN2=10;SN3=12;  %信噪比

SN=[SN1;SN2;SN3];

degrad=pi/180;  %將角度制化為弧度制    



%信號和噪聲

load shy3    %載入信號源

TT=1:cys;    %采樣數(shù)矩陣

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;



%計算協(xié)方差矩陣，特征值分解

tmp=-i*pi*sin(angle'.*degrad);  %相位差

tmp1=[0:m-1]';        %信號到達(dá)陣元的相位差組成的向量

A=tmp1*tmp;

A
文件尾:
1=exp(A);            %m*p維方向矩陣

X=(A*S1)+E;           %第k個陣元上的接收信號

Rxx=X*(X/cys)';       %接收信號的協(xié)方差

[E,S,V]=svd(Rxx);     %計算特征值

Vs=E(:,1:p);          

Vu=E(:,p+1:m);



%計算空間譜函數(shù)

th2=[-90:90]';         %線性信號角度范圍

tmp2=-i*pi*sin(th2'.*degrad);  %[-90 90]范圍內(nèi)的相位差組成的向量

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軸對數(shù)坐標(biāo)曲線

title('MUSIC算法');    %圖表標(biāo)題，x軸y軸標(biāo)簽及范圍

xlable('角度（弧度制）');

ylable('值');

axis([-90 90 0.1 1e4]);