波束赋形MATLAB代码实现

news/2025/12/1 11:17:32/文章来源:https://www.cnblogs.com/I989809/p/19292068

一、波束赋形基本原理

波束赋形(Beamforming)是一种通过控制天线阵列中各阵元的幅度和相位，在特定方向上增强信号强度、抑制干扰的空间滤波技术。核心原理是利用相控阵原理，通过调整各阵元信号的加权系数，使阵列在目标方向形成相干叠加，在其他方向形成相消干涉。

关键技术参数：

阵列因子(AF)：描述阵列方向特性的数学表达式
导向矢量(Steering Vector)：描述信号到达各阵元的相位差
加权向量(Weight Vector)：控制各阵元信号的幅度和相位
空间谱(Spatial Spectrum)：表示不同方向的信号强度

二、MATLAB代码实现

1. 均匀线阵(ULA)波束赋形

clear; clc; close all;%% ==================== 参数设置 ====================
c = 3e8;             % 光速 (m/s)
fc = 2.4e9;          % 载波频率 (Hz)
lambda = c/fc;       % 波长 (m)
d = lambda/2;        % 阵元间距 (m)
N = 8;               % 阵元数量
theta_desired = 30;  % 期望信号方向 (度)
theta_interf = [60, -20]; % 干扰信号方向 (度)
SNR = 20;            % 信噪比 (dB)
INR = 30;            % 干噪比 (dB)
fs = 10e6;           % 采样频率 (Hz)
t = 0:1/fs:1e-3;      % 时间向量 (1ms)
A_desired = 1;       % 期望信号幅度
A_interf = [0.8, 0.6]; % 干扰信号幅度%% ==================== 阵列流形向量 ====================
% 导向矢量函数
steering_vector = @(theta, N, d, lambda) exp(-1j*2*pi*d*(0:N-1)'*sind(theta)/lambda);% 期望信号导向矢量
a_desired = steering_vector(theta_desired, N, d, lambda);% 干扰信号导向矢量
a_interf1 = steering_vector(theta_interf(1), N, d, lambda);
a_interf2 = steering_vector(theta_interf(2), N, d, lambda);%% ==================== 信号模型 ====================
% 期望信号 (正弦波)
s_desired = A_desired * exp(1j*2*pi*fc*0.1*t); % 载波调制% 干扰信号
s_interf1 = A_interf(1) * exp(1j*2*pi*fc*0.2*t);
s_interf2 = A_interf(2) * exp(1j*2*pi*fc*0.3*t);% 接收信号
X = a_desired*s_desired + a_interf1*s_interf1 + a_interf2*s_interf2;% 添加噪声
noise_power = 1/(10^(SNR/10));
noise = sqrt(noise_power/2)*(randn(N, length(t)) + 1j*randn(N, length(t)));
X = X + noise;%% ==================== 波束赋形算法 ====================
% 1. 延迟求和波束形成器 (Conventional Beamformer)
w_ds = a_desired / N;  % 均匀加权% 2. 最小方差无失真响应波束形成器 (MVDR)
R = X*X'/length(t);    % 接收信号协方差矩阵
w_mvdr = inv(R)*a_desired/(a_desired'*inv(R)*a_desired); % MVDR权重% 3. 最大信噪比波束形成器 (Max-SNR)
w_snr = inv(R)*a_desired;% 4. 线性约束最小方差波束形成器 (LCMV)
C = [a_desired, a_interf1, a_interf2]; % 约束矩阵
f = [1; 0; 0];                       % 约束向量
w_lcmv = inv(R)*C*(inv(C'*inv(R)*C))*f;%% ==================== 波束形成输出 ====================
% 延迟求和输出
y_ds = w_ds'*X;% MVDR输出
y_mvdr = w_mvdr'*X;% Max-SNR输出
y_snr = w_snr'*X;% LCMV输出
y_lcmv = w_lcmv'*X;%% ==================== 波束方向图绘制 ====================
% 扫描角度范围
theta_scan = -90:0.1:90;
AF_ds = zeros(size(theta_scan));
AF_mvdr = zeros(size(theta_scan));
AF_snr = zeros(size(theta_scan));
AF_lcmv = zeros(size(theta_scan));for i = 1:length(theta_scan)a = steering_vector(theta_scan(i), N, d, lambda);AF_ds(i) = abs(w_ds'*a);AF_mvdr(i) = abs(w_mvdr'*a);AF_snr(i) = abs(w_snr'*a);AF_lcmv(i) = abs(w_lcmv'*a);
end% 归一化
AF_ds = 20*log10(AF_ds/max(AF_ds));
AF_mvdr = 20*log10(AF_mvdr/max(AF_mvdr));
AF_snr = 20*log10(AF_snr/max(AF_snr));
AF_lcmv = 20*log10(AF_lcmv/max(AF_lcmv));% 绘制波束方向图
figure;
plot(theta_scan, AF_ds, 'b', 'LineWidth', 1.5); hold on;
plot(theta_scan, AF_mvdr, 'r', 'LineWidth', 1.5);
plot(theta_scan, AF_snr, 'g', 'LineWidth', 1.5);
plot(theta_scan, AF_lcmv, 'm', 'LineWidth', 1.5);
xline(theta_desired, '--k', 'LineWidth', 1.5);
xline(theta_interf(1), '--k', 'LineWidth', 1.5);
xline(theta_interf(2), '--k', 'LineWidth', 1.5);
xlabel('角度 (度)');
ylabel('归一化增益 (dB)');
title('波束方向图比较');
legend('延迟求和', 'MVDR', 'Max-SNR', 'LCMV', '期望方向', '干扰方向');
grid on;
axis([-90 90 -40 0]);

2. 平面阵波束赋形（2D方向图）

%% ==================== 平面阵参数设置 ====================
M = 8;               % 行数
N = 8;               % 列数
dx = lambda/2;       % 水平间距
dy = lambda/2;       % 垂直间距
theta_desired = [30, 20]; % 期望方位角和俯仰角 (度)
theta_interf = [60, 30; -20, 10]; % 干扰方向% 生成平面阵导向矢量
[az_grid, el_grid] = meshgrid(-90:1:90, -90:1:90);
steering_vector_2d = @(az, el) exp(-1j*2*pi*(dx*(0:M-1)'*sind(az)*cosd(el) + dy*(0:N-1)'*sind(az)*sind(el))/lambda);% 期望信号导向矢量
a_desired = steering_vector_2d(theta_desired(1), theta_desired(2));
a_desired = kron(a_desired(:), ones(N,1)); % 展开为MN×1向量% 干扰信号导向矢量
a_interf1 = steering_vector_2d(theta_interf(1,1), theta_interf(1,2));
a_interf1 = kron(a_interf1(:), ones(N,1));
a_interf2 = steering_vector_2d(theta_interf(2,1), theta_interf(2,2));
a_interf2 = kron(a_interf2(:), ones(N,1));%% ==================== 平面阵波束方向图 ====================
AF_2d = zeros(size(az_grid));for i = 1:size(az_grid,1)for j = 1:size(az_grid,2)a = steering_vector_2d(az_grid(i,j), el_grid(i,j));a = kron(a(:), ones(N,1));AF_2d(i,j) = abs(w_mvdr'*a); % 使用MVDR权重end
end% 归一化并转换为dB
AF_2d = 20*log10(AF_2d/max(AF_2d(:)));% 绘制3D波束方向图
figure;
surf(az_grid, el_grid, AF_2d, 'EdgeColor', 'none');
view(45, 30);
xlabel('方位角 (度)');
ylabel('俯仰角 (度)');
zlabel('增益 (dB)');
title('平面阵波束方向图');
colorbar;
colormap jet;
shading interp;

3. 自适应波束赋形（LMS算法）

%% ==================== LMS自适应波束形成 ====================
mu = 0.01;            % 步长参数
w_lms = zeros(N, 1);  % 初始化权重向量% 期望信号（已知）
d = a_desired'*X;     % 期望输出% LMS算法迭代
y_lms = zeros(1, length(t));
e = zeros(1, length(t));for n = 1:length(t)x = X(:, n);      % 当前输入向量y_lms(n) = w_lms'*x; % 波束形成器输出e(n) = d(n) - y_lms(n); % 误差信号w_lms = w_lms + mu*e(n)*conj(x); % 权重更新
end% 绘制LMS收敛曲线
figure;
subplot(2,1,1);
plot(t, real(d), 'b', t, real(y_lms), 'r');
legend('期望信号', 'LMS输出');
xlabel('时间 (s)');
ylabel('幅度');
title('LMS波束形成器输出');subplot(2,1,2);
plot(t, e);
xlabel('时间 (s)');
ylabel('误差');
title('LMS误差信号');

三、波束赋形性能评估

1. 方向图指标计算

%% ==================== 方向图性能指标 ====================
% 主瓣宽度
main_lobe_width_ds = calculate_beamwidth(theta_scan, AF_ds, theta_desired);
main_lobe_width_mvdr = calculate_beamwidth(theta_scan, AF_mvdr, theta_desired);% 旁瓣电平
sidelobe_level_ds = calculate_sidelobe_level(theta_scan, AF_ds, theta_desired);
sidelobe_level_mvdr = calculate_sidelobe_level(theta_scan, AF_mvdr, theta_desired);% 零陷深度
null_depth_mvdr = calculate_null_depth(theta_scan, AF_mvdr, theta_interf);fprintf('延迟求和波束形成器:\n');
fprintf('  主瓣宽度: %.2f°\n', main_lobe_width_ds);
fprintf('  旁瓣电平: %.2f dB\n', sidelobe_level_ds);fprintf('\nMVDR波束形成器:\n');
fprintf('  主瓣宽度: %.2f°\n', main_lobe_width_mvdr);
fprintf('  旁瓣电平: %.2f dB\n', sidelobe_level_mvdr);
fprintf('  零陷深度: %.2f dB (%.1f°), %.2f dB (%.1f°)\n', ...null_depth_mvdr(1), theta_interf(1), ...null_depth_mvdr(2), theta_interf(2));% 计算主瓣宽度函数
function width = calculate_beamwidth(theta, AF, theta_desired)[max_val, idx] = max(AF);left_idx = find(AF(1:idx) <= max_val-3, 1, 'last');right_idx = find(AF(idx:end) <= max_val-3, 1, 'first') + idx - 1;if isempty(left_idx), left_idx = 1; endif isempty(right_idx), right_idx = length(theta); endwidth = theta(right_idx) - theta(left_idx);
end% 计算旁瓣电平函数
function level = calculate_sidelobe_level(theta, AF, theta_desired)[max_val, idx] = max(AF);main_lobe_mask = (theta > theta_desired-10) & (theta < theta_desired+10);sidelobe_region = AF;sidelobe_region(main_lobe_mask) = -Inf;level = max(sidelobe_region);
end% 计算零陷深度函数
function depth = calculate_null_depth(theta, AF, theta_interf)depth = zeros(size(theta_interf));for i = 1:length(theta_interf)[~, idx] = min(abs(theta - theta_interf(i)));depth(i) = AF(idx);end
end

2. 输出信干噪比(SINR)分析

%% ==================== SINR分析 ====================
% 计算输出功率
P_out_ds = mean(abs(y_ds).^2);
P_out_mvdr = mean(abs(y_mvdr).^2);
P_out_snr = mean(abs(y_snr).^2);
P_out_lcmv = mean(abs(y_lcmv).^2);% 计算干扰加噪声功率
interf_noise_power = noise_power + sum(A_interf.^2);% 计算输入SINR
SINR_in = 10*log10(A_desired^2 / interf_noise_power);% 计算输出SINR
SINR_out_ds = 10*log10(A_desired^2 / (P_out_ds - A_desired^2));
SINR_out_mvdr = 10*log10(A_desired^2 / (P_out_mvdr - A_desired^2));
SINR_out_snr = 10*log10(A_desired^2 / (P_out_snr - A_desired^2));
SINR_out_lcmv = 10*log10(A_desired^2 / (P_out_lcmv - A_desired^2));fprintf('\n输入SINR: %.2f dB\n', SINR_in);
fprintf('延迟求和输出SINR: %.2f dB (改善: %.2f dB)\n', SINR_out_ds, SINR_out_ds - SINR_in);
fprintf('MVDR输出SINR: %.2f dB (改善: %.2f dB)\n', SINR_out_mvdr, SINR_out_mvdr - SINR_in);
fprintf('Max-SNR输出SINR: %.2f dB (改善: %.2f dB)\n', SINR_out_snr, SINR_out_snr - SINR_in);
fprintf('LCMV输出SINR: %.2f dB (改善: %.2f dB)\n', SINR_out_lcmv, SINR_out_lcmv - SINR_in);

四、波束赋形应用场景

1. 5G Massive MIMO系统

%% ==================== 5G Massive MIMO波束赋形 ====================
% 参数设置
N_BS = 64;           % 基站天线数
N_UE = 8;            % 用户设备天线数
fc = 28e9;           % 毫米波频率 (GHz)
lambda = c/fc;
d = lambda/2;
theta_ue = rand(1, N_UE)*60-30; % 用户方位角 (-30°~30°)% 生成用户信道矩阵
H = zeros(N_BS, N_UE);
for k = 1:N_UEa = steering_vector(theta_ue(k), N_BS, d, lambda);H(:,k) = (randn(N_BS,1) + 1j*randn(N_BS,1))/sqrt(2) .* a; % Rayleigh衰落
end% ZF预编码
W_zf = H*(H'*H)\eye(N_UE); % 迫零预编码% MMSE预编码
W_mmse = H*(H'*H + 0.1*eye(N_UE))\eye(N_UE); % MMSE预编码% 波束方向图
theta_scan = -90:0.1:90;
AF_zf = zeros(size(theta_scan));
AF_mmse = zeros(size(theta_scan));for i = 1:length(theta_scan)a = steering_vector(theta_scan(i), N_BS, d, lambda);AF_zf(i) = 20*log10(abs(a'*W_zf(:,1))); % 第一个用户的波束AF_mmse(i) = 20*log10(abs(a'*W_mmse(:,1)));
end% 绘制波束方向图
figure;
plot(theta_scan, AF_zf, 'b', theta_scan, AF_mmse, 'r');
xlabel('角度 (度)');
ylabel('增益 (dB)');
title('5G Massive MIMO波束赋形');
legend('ZF预编码', 'MMSE预编码');
grid on;
axis([-90 90 -40 0]);

2. 雷达系统波束赋形

%% ==================== 雷达系统波束赋形 ====================
% 参数设置
N = 16;              % 阵元数量
theta_target = 10;   % 目标方向 (度)
theta_clutter = 30;  % 杂波方向 (度)
PRF = 2e3;           % 脉冲重复频率 (Hz)
T = 1/PRF;           % 脉冲间隔
num_pulses = 64;     % 脉冲数% 生成雷达回波信号
target_sig = exp(1j*2*pi*1e6*(0:num_pulses-1)*T); % 目标信号
clutter_sig = exp(1j*2*pi*0.5e6*(0:num_pulses-1)*T); % 杂波信号% 接收信号矩阵
X = steering_vector(theta_target, N, d, lambda)*target_sig + ...steering_vector(theta_clutter, N, d, lambda)*clutter_sig + ...0.1*(randn(N, num_pulses) + 1j*randn(N, num_pulses)); % 噪声% 动目标显示(MTI)波束形成
w_mti = ones(N,1) - [1; zeros(N-1,1)]; % 相邻脉冲相减% 多普勒处理
range_profile = zeros(N, num_pulses);
for n = 1:num_pulsesrange_profile(:,n) = fftshift(fft(X(:,n)));
end% 绘制距离-多普勒图
figure;
imagesc((0:num_pulses-1)*PRF/num_pulses, linspace(-1/(2*T), 1/(2*T), N), 20*log10(abs(range_profile)));
xlabel('多普勒频率 (Hz)');
ylabel('距离门');
title('雷达距离-多普勒图');
colorbar;

参考代码波束赋形的matlab代码 www.youwenfan.com/contentcnm/83082.html

五、总结与扩展

1. 不同波束赋形算法对比

算法	优点	缺点	适用场景
延迟求和	简单易实现，计算量小	旁瓣高，干扰抑制差	简单定向通信
MVDR	高分辨率，干扰抑制强	需要准确协方差矩阵估计	雷达、抗干扰通信
LCMV	多约束，灵活性强	计算复杂，需矩阵求逆	多目标跟踪
LMS/RLS	自适应，无需统计先验	收敛速度慢，稳态误差	时变环境
ZF/MMSE	高容量，抗多用户干扰	需要信道状态信息	大规模MIMO