clear all;
derad = pi/180; % 角度->弧度
radeg = 180/pi; % 弧度->角度
k = 16; % 阵元数
n = 500; % 采样次数
M = 500; % 蒙特卡罗次数
snr = [-4:12]; % 信噪比
ns = 1; % 信源数
d = 1; % 阵元间距
lambda = 2; % 波长
theta = [-30 ]; % 入射角
A = exp(-j*2*pi*[0:15].'*d*sin(theta*derad)/lambda); % 方向矢量
thetam = zeros(1,M);
thetab = zeros(1,M);
thetac = zeros(1,M);
RMSEm = zeros(1,length(snr));
RMSEb = zeros(1,length(snr));
RMSEc = zeros(1,length(snr));
for is = 1:length(snr)
for ii = 1:M
S = randn(ns,n); % 信源信号
X = A*S; % 接收信号
X1 = awgn(X,snr(is),'measured'); % 添加高斯白噪声
Rxx = X1*X1'/n; % 构建接收信号的协方差矩阵
IRxx = inv(Rxx); % 求逆矩阵
[P,V] = eig(Rxx); % 特征值分解
EV = diag(V)';
[EV,I] = sort(EV); % 将特征值从小到大排列
EVA = fliplr(EV); % 翻转，将特征值从大到小排列
EVV = fliplr(P(:,I)); % 对应特征矢量排序
beta = [-pi/2:pi/500:pi/2];
% 构造MUSIC谱函数
music = zeros(1,length(beta));
bartlet = zeros(1,length(beta));
capon = zeros(1,length(beta));
for i = 1:length(beta)
a = exp(-j*2*pi*[0:15].'*sin(beta(i))/lambda);
En = EVV(:,ns+1:k);
music(i) = 1/(a'*En*En'*a);
bartlet(i) = (a'*Rxx*a)/(a'*a);
capon(i) = inv(a'*IRxx*a);
end
music_abs = abs(music);
music_max = max(music_abs);
wm = find(music_abs==music_max);
thetam(ii) = beta(wm)*radeg;
Pmusic = 10*log10(music_abs/music_max);
bartlet_abs = abs(bartlet);
bartlet_max = max(bartlet_abs);
bm = find(bartlet_abs==bartlet_max);
thetab(ii) = beta(bm)*radeg;
Pbartlet = 10*log10(bartlet_abs/bartlet_max);
capon_abs = abs(capon);
capon_max = max(capon_abs);
cm = find(capon_abs==capon_max);
thetac(ii) = beta(cm)*radeg;
Pcapon = 10*log10(capon_abs/capon_max);
end
% 均方误差
RMSEm(is) = sqrt(sum((thetam-theta*ones(1,M)).^2)/M);
RMSEb(is) = sqrt(sum((thetab-theta*ones(1,M)).^2)/M);
RMSEc(is) = sqrt(sum((thetac-theta*ones(1,M)).^2)/M);
end
% 绘图
figure(1);
t = beta*radeg;
plot(t,Pmusic,'k');
xlabel('Angle in Degress(theta)');
ylabel('P(dB)');
title('Spectrum plot(snr=4dB)');
grid on
hold on;
plot(t,Pbartlet,'b');
plot(t,Pcapon,'r');
legend('Music','Bartlet','Capon');
figure(2);
plot(snr,RMSEm,snr,RMSEb,snr,RMSEc);
grid on;
legend('Music','Bartlet','Capon');
