麦克风阵列波束形成（8阵元圆阵）

热爱学习的任同学

于 2024-03-07 15:24:49 发布

阅读量627

点赞数 17

文章标签： matlab 人工智能

本文链接：https://blog.csdn.net/weixin_44583910/article/details/136531918

版权

麦克风阵列波束形成（圆阵）

引言
MATLAB波束形成程序
结果图
参考文献

引言

最近需要使用麦克风阵列，但是网上很少找到关于平面阵麦克风阵列的代码，对于我这种时间紧张的初学者而言很难受，吃了不少亏，现在把阶段性学习成果分享给大家，希望对后面刚接触这方面的同学有点启发。

MATLAB波束形成程序

创建基本阵列参数

c = 340;                        				%声速
R = 0.04354;                    				%阵元半径
M = 8;                          				%阵元个数
phin = 360/M;                   				%阵元位置角度
phin = 0:phin:360-phin;
zero_rho_z = zeros(1,M);
rho = [R*cosd(phin);R*sind(phin);zero_rho_z];	%直角坐标系下阵元坐标
snapshot_num = 400;             				%快拍数
f = 10000;										%目标频率
lambda = c/f;                   				%波长
k = 2*pi/lambda;                				%波数值

设定期望值及观察范围

theta0 = 0;										%目标俯仰角
phi0 = 180;										%目标方向角
theta = linspace(-90,90,361);   				%俯仰角范围
phi = linspace(0,360,361);      				%水平角范围

%***********若是需要抑制一定方向的干扰则加入下方代码***********%
%theta1 = 10;									%干扰1俯仰角
%phi1 = 80;										%干扰1水平角
%theta2 = -30;									%干扰2俯仰角
%phi2 = 200;									%干扰2水平角
%************************************************************%

创建数据协方差矩阵

kv0 = k*(-[sind(theta0)*cosd(phi0);sind(theta0)*sind(phi0);cosd(theta0)]);	%目标方向波数向量										
P0 = exp(-1j*kv0'*rho).';												  	%目标方向响应向量
tar_sig = wgn(1,snapshot_num,0);                                        	%目标信号
noise = (randn(M,snapshot_num)+randn(M,snapshot_num)*1i)/sqrt(2);       	%环境噪声
rec_sig = P0*tar_sig+noise;													%无干扰的输入信号

%***********若是需要抑制一定方向的干扰则加入下方代码***********%
%kv1 = k*(-[sind(theta1)*cosd(phi1);sind(theta1)*sind(phi1);cosd(theta1)]);	%干扰1方向波数向量
%kv2 = k*(-[sind(theta2)*cosd(phi2);sind(theta2)*sind(phi2);cosd(theta2)]);	%干扰2方向波数向量
%P1 = exp(-1j*kv1'*rho).';													%干扰1方向响应向量
%P2 = exp(-1j*kv2'*rho).';													%干扰2方向响应向量
%inf1 = wgn(1,snapshot_num,20);                                          	%干扰信号1
%inf2 = wgn(1,snapshot_num,15);                                          	%干扰信号2
%rec_sig = P0*tar_sig+P1*inf1+P2*inf2+noise;								%有干扰的输入信号
%************************************************************%

Rx = rec_sig*rec_sig'/snapshot_num;											%数据协方差矩阵

计算方向响应矩阵（或称阵列流形向量）

kv_all = k*(-[kron(sind(theta),cosd(phi));kron(sind(theta),sind(phi));reshape(repmat(cosd(theta)',1,361)',[1,130321])]);
P_all = exp(-1j*kv_all'*rho).';

RCB算法计算各传感器权值

epsilon0 = 2;                     			%导向向量误差范围上界
R_inv = inv(Rx); 							%求逆矩阵
cvx_begin quiet
    variable a(M) complex;
    minimize norm(sqrtm(R_inv)*a)
subject to
   norm(P0-a) <= epsilon0;
cvx_end
w = R_inv*a/(a'*R_inv*a);

绘制波束响应图

BF = reshape(w'*P_all,[361,361]).';			%计算波束响应并整理格式
BF = 20*log10(abs(BF)/max(max(abs(BF))));	%转成dB单位

figure(1)									%水平角为目标值时的图3刨面
theta0_index = find(theta == theta0);
plot(phi,BF(theta0_index,:));
xlabel('水平角 \phi °');
ylabel('增益大小 (dB)');
title('圆阵RCB');
xlim([0,360]);
ylim([-25,0]);
grid;

figure(2)									%俯仰角为目标值时的图3刨面
plot(theta,BF(:,phi0));
xlabel('俯仰角 \theta °');
ylabel('增益大小 (dB)');
title('圆阵RCB');
xlim([-90,90]);
ylim([-25,0]);
grid;

figure(3)									%直角坐标系中显示
mesh(phi,theta,BF)
xlabel('水平角 \phi °');
ylabel('俯仰角 \theta ° ');
zlabel('增益大小 (dB)');
title('圆环阵RCB三维视图');

figure(4)									%转换为极坐标系中显示
[ttt,ppp] = meshgrid(theta,phi);
BF = (BF-min(BF(:)))/(max(BF(:))-min(BF(:)));
x = BF.*sind(ttt').*cosd(ppp');
y = BF.*sind(ttt').*sind(ppp');
z = BF.*cosd(ttt');
xlabel('x');
ylabel('y');
zlabel('z');
mesh(x,y,z)