波束赋形MATLAB代码实现

一、波束赋形基本原理

波束赋形(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; end
    if isempty(right_idx), right_idx = length(theta); end
    width = 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_UE
    a = 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_pulses
    range_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/contentcsn/83082.html

五、总结与扩展

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

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

2. 扩展方向

1. 混合波束赋形：结合模拟和数字波束赋形，降低硬件复杂度
2. 深度学习波束赋形：使用神经网络优化波束形成权重
3. 毫米波波束赋形：针对高频段的大带宽、窄波束特性优化
4. 全息波束赋形：利用超材料实现连续孔径波束赋形

3. 工程实践建议

1. 阵列校准：定期进行阵列校准，补偿阵元位置误差
2. 通道均衡：对每个接收通道进行幅相一致性校准
3. 鲁棒设计：考虑导向矢量误差，采用鲁棒波束形成算法
4. 硬件实现：FPGA实现时注意定点量化效应

注意：实际系统中还需考虑互耦效应、阵元方向图、热噪声等非理想因素。本代码提供了波束赋形的基本实现框架，可根据具体需求进行修改和优化。