tic
%代码块
%% initial matlab workspace
clc
clear
close all
%% 坐标参数设置
H = 7; % 圆柱中轴线的高度
r = 3; % 圆柱半径
minx = -10; % 圆柱开端
maxx = 10; % 圆柱终端
% 信号参数设置
SN = 1024; % 快拍数
SNR = 15; % 信噪比,最低到-5db
f = 80; % 工作频率,最低降到100Hz(频率取极端低失效)
c = 1500; % 声速
% 信源设置
ph = [ pi, pi ];
doas_x = [-5, 5];
doas_y = r .* cos(ph);
doas_z = r .* sin(ph) + H;
% 提取信源数目
M = 2;
% 阵列参数设置
N = 16; % 阵元数量
%d = 0.5*c/f; % 单元间距
d = 2; %可调大调小2
% 阵列坐标生成
cx = (0:N-1)*d - (N-1)/2*d;
cy = zeros( 1, N );
cz = zeros( 1, N );
% 曲面生成
ph = linspace(0,2*pi,201);
x = linspace(minx,maxx,101);
[ph,xx] = meshgrid( ph, x );
% 生成对应的yy和zz
yy = r .* cos(ph);
zz = r .* sin(ph) + H;
% % 画图 check
% figure
% mesh( xx, yy, zz, -( xx.^2 ) )
%% 接收信号生成
% 生成随机信源信号
st = rand( M, SN) + 1i * rand(M,SN);
% 生成信源坐标至所有位置的距离rall
rall = pdist2( [cx',cy',cz'], [doas_x',doas_y',doas_z'] );
% 导向矩阵生成
k = 2*pi/(c/f); % 波数
A = exp( -1i*k.*rall )./rall;
% 生成接收信号x
x = A * st;
% 接收信号添加噪声
x = awgn( x, SNR, 'measured');
%% MUSIC DOA估计
% 生成协方差矩阵R
R = x*x'/SN;
R_inv=inv(R);
% 初始化空间谱
P = zeros( size(xx) );
% 迭代求空间谱
for i1 = 1:size(P,1)
for i2 = 1:size(P,2)
% 提取待估计点的坐标cxyz
xyz = [ xx(i1,i2), yy(i1,i2), zz(i1,i2) ];
% 生成导向矢量a
ra = pdist2( [cx',cy',cz'], xyz );
a = exp( -1i.*k.*ra )./ra;
%a = exp( -1i.*k.*ra )
% 计算P
P(i1,i2) = 1/(a'*R_inv*a);
end
end
% 空间谱归一化 并化为dB
% P = 1./real(P);
% P = P./max(max(P));
% P = 10.*log10( P );
P = 10*log10(abs(P)./max(max(abs(P))));
% 画图 空间谱
figure
Pmin = min(P(:));
for i1 = 1:size(P,1)
for i2 = 1:size(P,2)
if yy(i1,i2) > 0
P(i1,i2) = Pmin;
end
end
end
mesh( xx, yy, zz, P )
xlabel('x/m')
ylabel('y/m')
zlabel('z/m');
set(gca,'ylim',[0,4]);%y坐标轴范围
set(gca,'FontSize',13)
colorbar
colormap(jet)
toc
%disp(['运行时间: ',num2str(toc)]);
【matlab三维绘图-半圆柱(空间谱)】
猜你喜欢
转载自blog.csdn.net/weixin_44312889/article/details/128182865
今日推荐
周排行