% 仿真程序
clear all
format long %将数据显示为长整型科学计数
K = 128; %快拍数
doa = [5 40]/180*pi; %信号到达角
gamma = [pi/6 pi/4]; %gamma
w = [pi/4 pi/3]'; %信号频率
N = 8; %阵元数
P = length(w); %信号个数2
lambda = 150; %波长
d = lambda/2; %阵元间距
snr = 18; %信噪比
As = zeros(P,N); %创建一个2行N列的0矩阵
for k = 1:P
As(k,:) = exp(-j*2*pi*d*sin(doa(k))/lambda*[0:N-1]); %矩阵赋值,每一行对应一个theta值 [as1;as2]
end
Ap = zeros(P,2); %创建一个2行2列的0矩阵
for k = 1:P
Ap(k,:) = [-cos(gamma(k)) j*cos(doa(k))*sin(gamma(k))]; %矩阵赋值,每一行对应一个gamma值
end
A = zeros(P,2*N); %创建一个2行2N列的0矩阵
for k = 1:P
A(k,:) = kron(As(k,:),Ap(k,:)); %矩阵赋值,As的第k行kronAp的第k行
end
B = A'; %转置矩阵 2N*2
xx = 2*exp(j*(w*[1:K])); %仿真信号,有2乘128个 p35 2.4.21
x = B*xx; % 2N * 128
x = x + awgn(x,snr); %加入高斯白噪声
R = x*x'; %数据协方差矩阵
[U, V] = eig(R); %求R的特征值和特征向量（特征值分解）
UU = U(:,1:2*N-P); %估计噪声子空间
tem_theta = [-20:0.5:60]*pi/180;
tem_gamma = [0:0.5:80]*pi/180;
dgr_theta = [-20:0.5:60];
dgr_gamma = [0:0.5:80];
Pmusic=zeros(length(tem_theta),length(tem_gamma));
%谱峰搜索
for ii = 1:length(tem_theta) %361个点
for kk =1:length(tem_gamma) %181个点
%1*N
AAs = exp(-j*2*[0:N-1]*pi*d*sin(tem_theta(ii))/lambda);
AAp=[-cos(tem_gamma(kk)) j*cos(tem_theta(ii))*sin(tem_gamma(kk))];
AA=kron(AAs,AAp);
WW = AA*UU*UU'*AA';
Pmusic(ii,kk) = abs(1/ WW);
end
end
fPmusic = 10*log10(Pmusic/max(max(Pmusic))); %空间谱函数
[X,Y]=meshgrid(dgr_theta,dgr_gamma);
h=mesh(X',Y',fPmusic);colorbar;
set(h,'Linewidth',2);
xlabel('theta');
ylabel('gamma');
zlabel('magnitude(dB)');
figure;
contour(X',Y',fPmusic,7);
评论0