clear?all;
close?all;
clc;

source_number=3;%信元數
sensor_number=8;%陣元數
N_x=1024;?%信號長度
snapshot_number=N_x;%快拍數
w=[pi/6?pi/6?pi/6].‘;%信號頻率
l=0.15;%信號波長??
d=l/2;%陣元間距
snr=20;%信噪比

source_doa=[20?40?60];%三個信號的入射角度
A=[exp（-j*（0:sensor_number-1）*d*2*pi*sin（source_doa（1）*pi/180）/l）;exp（-j*（0:sensor_number-1）*d*2*pi*sin（source_doa（2）*pi/180）/l）;exp（-j*（0:sensor_number-1）*d*2*pi*sin（source_doa（3）*pi/180）/l）].‘;%陣列流型

s=sqrt（10.^（snr/10））*exp（j*w*[0:N_x-1]）;%仿真信號
%x=awgn（ssnr）;
x=A*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪聲后的陣列接收信號
%求四階累積量C4
temp1?=?kron（xconj（x））;
C4=temp1*temp1‘/（snapshot_number.^2）-temp1*temp1‘-kron（x*x‘/snapshot_numberconj（x*x‘）/snapshot_number）;
disp（C4）;
[VD]=eig（-C4）;%對C4進行特征分解
[lambda?index]=sort（diag（D））;
En=V（:index（1:（sensor_number.^2-source_number^2）））;
disp（En）;
searching_doa=-90:0.5:90;%線陣的搜索范圍為-90~90度
?for?i=1:length（searching_doa）
???a_theta=exp（-j*（0:sensor_number-1）‘*2*pi*d*sin（pi*searching_doa（i）/180）/l）;
???b_theta=kron（a_thetaconj（a_theta））;
???Pmusic（i）=1./（norm（b_theta‘*En））.^2;
?end
?figure;
plot（searching_doa10*log（Pmusic））;
xlabel（‘入射角/degree‘）;
ylabel（‘功率譜/dB‘）;
title（‘四階累積量MUSIC估計‘）;
grid?;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????1057??2018-09-23?19:30??music.m

?????文件???????1368??2012-08-27?19:25??PM.m

?????文件???????1317??2018-09-23?21:09??CUM_MUSIC.m

-----------?---------??----------?-----??----

?????????????????3742????????????????????3