ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

MATLAB天线仿真实战:从方向图建模到阵列波束形成

MATLAB天线仿真实战:从方向图建模到阵列波束形成 简介本资源是一份面向通信工程、电磁场与微波技术方向初学者及课程设计者的MATLAB阵列天线仿真轻量工具包聚焦天线方向图绘制、波束形成原理验证与阵列参数快速分析等核心教学与实践需求。压缩包为RAR格式仅含1个MATLAB脚本文件.m体积仅620B代码简洁可读适用于Matlab R2018a及以上版本便于直接运行、修改与教学演示。已有273人学习下载反映出其在高校实验课、课程设计及自学入门阶段的实用价值。用户可直接调用该脚本完成均匀线阵建模、激励相位/幅度设置、空间方向图实时绘制并直观观察主瓣宽度、副瓣电平及波束指向变化规律是理解阵列天线基本原理与Matlab电磁仿真入门的高效辅助工具。1. 从“gbo.rar”到天线仿真一个MATLAB工程师的实战复盘最近在整理一个老项目时翻出了一个名为“gbo.rar”的压缩包。解压开来里面是几年前用MATLAB写的一堆关于天线方向图、波束形成和阵列天线仿真的脚本和函数。看着这些代码当时为了搞懂一个公式或者调通一个仿真模型在实验室熬到深夜的场景又历历在目。天线仿真尤其是阵列天线仿真对于通信、雷达、电子对抗等领域的朋友来说是绕不开的基本功。它不像纯理论推导那么抽象也不像硬件调试那么“玄学”它恰恰是连接理论与实践的桥梁——在电脑上用数学和算法去预测和设计一个真实天线的辐射特性。今天我就以这个尘封的“gbo.rar”项目为引子和大家系统地聊聊如何用MATLAB这把“瑞士军刀”从零开始构建一套属于自己的天线仿真工具箱。我们会聚焦于三个核心天线方向图、波束形成与阵列天线仿真。无论你是正在学习《天线原理》的学生还是需要快速验证天线设计方案的工程师这篇文章都将带你绕过我当年踩过的那些坑直击要害手把手教你如何用MATLAB实现从单个天线单元到复杂阵列的完整仿真链路。我们会从最基础的辐射模型讲起逐步深入到波束扫描、自适应调零等高级话题并提供可直接运行的代码片段和避坑指南。2. 天线方向图辐射特性的“指纹”与MATLAB建模天线方向图通俗讲就是天线向空间各个方向辐射或接收电磁波能力的图形化表示。它是天线的“指纹”决定了天线“看”向哪里、能“看”多清楚。在MATLAB里仿真方向图核心是建立正确的数学模型。2.1 理论基础从点源到实际天线模型仿真始于模型。最简单的模型是各向同性点源它在所有方向均匀辐射其方向图函数F(θ, φ) 1。这是一个理想的参考基准但实际天线不存在。更常用的是**电基本振子赫兹偶极子**模型其远场方向图函数为F(θ) sin(θ)在垂直于振子轴的方向θ90°辐射最强沿振子轴方向θ0°, 180°辐射为零。这个模型是理解方向性概念的起点。对于更实际的天线如半波偶极子、微带贴片天线等其方向图函数更为复杂往往由理论公式或实测数据拟合得到。在MATLAB中我们通常用一个关于球坐标角度θ俯仰角和φ方位角的函数F(θ, φ)来表示它。这个函数通常是复数包含了幅度和相位信息。幅度图表示功率密度分布相位图表示波前的空间变化。2.2 MATLAB实现绘制2D与3D方向图有了方向图函数绘图就水到渠成。MATLAB的polarplot、plot和patternCustom函数是得力工具。1. 绘制二维方向图如E面或H面假设我们有一个半波偶极子天线其E面包含振子轴和最大辐射方向的平面方向图近似为F(θ) cos((pi/2)*cos(θ))/sin(θ)。我们可以这样可视化% 定义角度范围弧度 theta linspace(0, 2*pi, 361); % 计算方向图函数值避免除零错误 F zeros(size(theta)); for i 1:length(theta) if abs(sin(theta(i))) 1e-10 % 避免sin(0)或sin(pi)导致的除零 F(i) abs(cos((pi/2)*cos(theta(i))) / sin(theta(i))); else F(i) 0; % 在理论上为零的点手动赋值 end end % 归一化 F F / max(F); % 使用极坐标绘图 figure; polarplot(theta, F, ‘LineWidth‘, 2); title(‘半波偶极子E面方向图极坐标‘); % 或者使用直角坐标绘图 figure; plot(theta*180/pi, 20*log10(F)); % 转换为度和对数尺度dB xlabel(‘角度 (度)‘); ylabel(‘增益 (dBi)‘); title(‘半波偶极子E面方向图直角坐标‘); grid on; xlim([0, 360]);注意直接计算cos((pi/2)*cos(θ))/sin(θ)在θ0或π时会出现NaN非数。在实际编程中必须加入判断语句处理这些奇异点或者使用MATLAB的sinc函数等数值稳定的写法。这是我早期代码里最常见的bug之一。2. 绘制三维方向图三维方向图能全面展示空间辐射特性。patternCustom函数非常便捷但它需要的是直角坐标系下的点数据。通常我们先在球坐标系下生成网格计算方向图值再转换。% 生成球坐标网格 [theta, phi] meshgrid(linspace(0, pi, 181), linspace(0, 2*pi, 361)); % 计算一个示例方向图例如一个具有方向性的贴片天线模型 % 假设在俯仰面和方位面都有余弦类分布 F abs(cos(theta).^2 .* cos(phi).^2); % 示例函数非真实物理公式 % 将球坐标转换为直角坐标用于绘图 [x, y, z] sph2cart(phi, pi/2 - theta, F); % 注意MATLAB球坐标定义azimuth, elevation, r % 使用patternCustom绘制需要将矩阵展平 figure; patternCustom(F(:), theta(:)*180/pi, phi(:)*180/pi); title(‘三维天线方向图示例‘);对于更复杂的自定义绘图可以使用mesh或surf函数并将F作为高度或颜色值映射到球面上。2.3 关键参数提取波束宽度与旁瓣电平方向图不仅是看着好看更要从中提取关键工程参数。最核心的两个是半功率波束宽度HPBW方向图主瓣功率下降到最大值一半-3dB时所对应的角度宽度。第一旁瓣电平SLL主瓣旁边第一个副瓣的峰值电平相对于主瓣峰值的差值通常为负的dB值。手动从图形上读取这些值既不准又麻烦。我们需要编写算法自动提取。% 假设我们已经有了一个一维方向图数组 F_dB单位dB对应角度数组 ang单位度 % 步骤1找到主瓣最大值及其位置 [max_gain, max_idx] max(F_dB); % 步骤2找到主瓣两侧-3dB点的位置 half_power max_gain - 3; % 向左搜索 idx_left max_idx; while idx_left 1 F_dB(idx_left) half_power idx_left idx_left - 1; end % 线性插值以获得更精确的左-3dB点角度 if idx_left max_idx ang_left interp1(F_dB(idx_left:idx_left1), ang(idx_left:idx_left1), half_power); else ang_left ang(idx_left); end % 向右搜索 idx_right max_idx; while idx_right length(F_dB) F_dB(idx_right) half_power idx_right idx_right 1; end if idx_right max_idx ang_right interp1(F_dB(idx_right-1:idx_right), ang(idx_right-1:idx_right), half_power); else ang_right ang(idx_right); end % 计算HPBW HPBW ang_right - ang_left; fprintf(‘半功率波束宽度HPBW为%.2f 度\n‘, HPBW); % 步骤3寻找第一旁瓣电平SLL % 需要排除主瓣区域。一个简单方法是定义主瓣区域为最大值左右各扩展HPBW/2的范围 mainlobe_start max_idx - round(HPBW/(ang(2)-ang(1))/2); mainlobe_end max_idx round(HPBW/(ang(2)-ang(1))/2); mainlobe_region false(size(F_dB)); mainlobe_region(max(1,mainlobe_start):min(length(F_dB),mainlobe_end)) true; % 在主瓣区域外寻找最大值即为第一旁瓣峰值 F_dB_excluding_main F_dB; F_dB_excluding_main(mainlobe_region) -inf; [first_sidelobe_gain, ~] max(F_dB_excluding_main); SLL first_sidelobe_gain - max_gain; fprintf(‘第一旁瓣电平SLL为%.2f dB\n‘, SLL);实操心得自动提取算法的鲁棒性至关重要。实际仿真或测量得到的方向图数据可能有噪声主瓣不一定在0度旁瓣可能不止一个。上述简单算法在理想情况下工作良好但对于复杂情况可能需要更稳健的方法如使用findpeaks函数配合参数来识别主瓣和旁瓣或者对数据进行平滑处理。在“gbo.rar”的后期版本中我专门写了一个analyzePattern函数来封装这些逻辑并加入了异常处理。3. 阵列天线基础从单元到“军团”的合成法则单个天线的能力是有限的。阵列天线通过将多个完全相同的天线单元阵元按一定规则排列并通过控制各阵元的馈电幅度和相位能够合成出远优于单个单元的方向图特性实现波束扫描、赋形、零点控制等功能。这是现代相控阵雷达和5G Massive MIMO技术的基石。3.1 阵列因子理解方向图合成的钥匙阵列天线的总方向图根据方向图乘积原理等于单元方向图Element Pattern和阵列因子Array Factor, AF的乘积。单元方向图描述了单个阵元自身的辐射特性而阵列因子完全由阵元的几何排列和馈电激励决定与单元本身无关。这为我们分析阵列特性提供了极大的便利我们可以先研究理想点源阵列的阵列因子再乘以实际单元方向图。对于一个由N个阵元组成的直线阵假设阵元等间距d排列在x轴上第n个阵元的激励电流为复数I_n A_n * exp(1j * φ_n)其中A_n是幅度φ_n是相位。那么该直线阵的阵列因子为AF(θ, φ) Σ_{n0}^{N-1} [ I_n * exp(j * k * d * n * sinθ cosφ) ]其中k 2π/λ是波数λ是波长。对于更复杂的二维或三维面阵公式会扩展为双重或三重求和但核心思想不变阵列因子是各阵元在远场观察点贡献的复振幅的相干叠加。3.2 MATLAB仿真均匀直线阵与方向图特性让我们从最简单的均匀直线阵Uniform Linear Array, ULA开始。所有阵元等幅同相馈电A_n1, φ_n0间距为d。% 参数设置 f 2.4e9; % 频率 2.4 GHz c 3e8; % 光速 lambda c / f; % 波长 d 0.5 * lambda; % 阵元间距通常取半波长以避免栅瓣 N 8; % 阵元数量 % 生成角度向量这里只看x-z平面即φ0的情况 theta_deg linspace(-90, 90, 1801); theta deg2rad(theta_deg); % 计算波数 k 2 * pi / lambda; % 计算阵列因子 AF zeros(size(theta)); for n 0:N-1 % 阵元位置x_n n * d % 波程差引起的相位差k * d * n * sin(theta) 因为φ0, cosφ1 AF AF exp(1j * k * d * n * sin(theta)); end AF_mag abs(AF); % 幅度 AF_dB 20*log10(AF_mag / max(AF_mag)); % 归一化dB值 % 绘图 figure; plot(theta_deg, AF_dB, ‘b-‘, ‘LineWidth‘, 2); xlabel(‘角度 (度)‘); ylabel(‘归一化阵列因子 (dB)‘); title([‘均匀直线阵方向图 (N‘, num2str(N), ‘, d‘, num2str(d/lambda), ‘λ)‘]); grid on; xlim([-90, 90]); ylim([-40, 0]);运行这段代码你会看到一个典型的多波束方向图。主瓣位于0度法线方向两侧对称分布着多个旁瓣。阵元数N主要影响主瓣宽度N越大主瓣越窄阵元间距d主要影响栅瓣grating lobe的出现。当d λ时在某些角度会出现与主瓣幅度相当的栅瓣这是要极力避免的。3.3 波束扫描原理相位控制的艺术阵列天线最迷人的特性之一是电扫描——不动天线只改变馈电相位就能让波束指向任意方向。假设我们希望主瓣指向θ0方向。根据前面阵列因子的公式为了使所有阵元在θ0方向同相叠加需要补偿因空间位置不同带来的波程差。这要求第n个阵元的激励相位为φ_n -k * d * n * sin(θ0)在MATLAB中实现扫描非常简单只需在计算阵列因子时为每个阵元乘上这个相位补偿项即可。% 接上段代码参数 theta0_deg 30; % 期望的波束指向角度 theta0 deg2rad(theta0_deg); % 计算扫描所需的激励相位 phase_shift -k * d * (0:N-1)‘ * sin(theta0); % 列向量 % 重新计算阵列因子带扫描相位 AF_scan zeros(size(theta)); for n 0:N-1 AF_scan AF_scan exp(1j * (k * d * n * sin(theta) phase_shift(n1))); end AF_scan_dB 20*log10(abs(AF_scan) / max(abs(AF_scan))); % 绘图对比 figure; hold on; plot(theta_deg, AF_dB, ‘b-‘, ‘DisplayName‘, ‘法向波束 (θ00°)‘); plot(theta_deg, AF_scan_dB, ‘r--‘, ‘LineWidth‘, 2, ‘DisplayName‘, [‘扫描波束 (θ0‘, num2str(theta0_deg), ‘°)‘]); xlabel(‘角度 (度)‘); ylabel(‘归一化阵列因子 (dB)‘); title(‘均匀直线阵波束扫描效果‘); legend(‘show‘); grid on; xlim([-90, 90]); ylim([-40, 0]); hold off;关键点剖析波束扫描的本质是通过线性相位梯度来“扳直”波前。在θ0方向来自各阵元的信号经过空间波程差和馈电相位补偿后到达远场观察点时相位完全一致实现同相叠加产生主瓣。改变θ0就改变了相位梯度从而改变了波束指向。这就是相控阵Phased Array的核心思想。在仿真中务必注意相位φ_n的计算符号一个正负号的错误会导致波束指向完全相反的方向。我早期的代码就曾因此调试了半天。4. 波束形成算法从固定权值到自适应优化上一节的均匀加权和固定相位扫描是最基础的波束形成Beamforming。但实际应用中我们往往对方向图有更复杂的要求比如更低的旁瓣、在干扰方向形成零点、或者根据环境自适应调整波束。这就引出了波束形成权值的设计问题。4.1 窗函数法抑制旁瓣的经典手段均匀加权矩形窗的阵列因子旁瓣较高第一旁瓣约-13dB。借鉴数字信号处理中的窗函数法对阵列各阵元的激励幅度进行加权即采用非均匀幅度分布可以有效抑制旁瓣代价是主瓣会略微展宽。常用的窗函数有汉宁Hanning、汉明Hamming、切比雪夫Chebyshev、泰勒Taylor等。% 参数设置 (沿用之前的ULA) N 16; d 0.5; theta deg2rad(linspace(-90, 90, 1801)); % 1. 均匀加权矩形窗 w_uniform ones(N, 1); % 2. 汉宁窗加权 w_hann hann(N); % MATLAB函数返回N点汉宁窗 % 3. 道尔夫-切比雪夫加权指定旁瓣电平 SLL_dB -30; % 期望的旁瓣电平为-30dB w_cheb chebwin(N, -SLL_dB); % chebwin需要旁瓣衰减值正数 % 计算并绘制方向图 figure; hold on; colors {‘b‘, ‘r‘, ‘g‘, ‘m‘}; windows {w_uniform, w_hann, w_cheb}; window_names {‘均匀窗‘, ‘汉宁窗‘, [‘切比雪夫窗(‘, num2str(-SLL_dB), ‘dB)‘]}; for i 1:length(windows) w windows{i}; AF zeros(size(theta)); for n 0:N-1 AF AF w(n1) * exp(1j * 2*pi/lambda * d * n * sin(theta)); end AF_dB 20*log10(abs(AF)/max(abs(AF))); plot(rad2deg(theta), AF_dB, colors{i}, ‘DisplayName‘, window_names{i}); end xlabel(‘角度 (度)‘); ylabel(‘归一化方向图 (dB)‘); title(‘不同窗函数加权对ULA方向图的影响 (N16)‘); legend(‘show‘, ‘Location‘, ‘best‘); grid on; xlim([-90, 90]); ylim([-80, 0]); hold off;运行代码可以清晰看到汉宁窗大幅降低了旁瓣但主瓣变宽切比雪夫窗则在指定旁瓣电平下实现了最窄的主瓣宽度是一种最优设计。4.2 自适应波束形成MVDR算法实例当存在已知方向的强干扰时我们希望在天线方向图的干扰方向形成很深的零点同时保持期望信号方向的增益。这需要自适应地计算权向量。最小方差无失真响应MVDR波束形成器是一个经典算法其目标是在保证期望方向增益为1的前提下使阵列输出的总功率最小从而抑制干扰和噪声。假设我们有一个N元阵列期望信号来自方向θ_s有K个干扰来自方向θ_i1, θ_i2, ...。阵列接收数据的协方差矩阵为R。MVDR的权向量w由下式给出w (R^{-1} * a(θ_s)) / (a(θ_s)^H * R^{-1} * a(θ_s))其中a(θ_s)是期望信号方向的导向矢量Steering Vector对于ULAa(θ) [1, exp(j*k*d*sinθ), ..., exp(j*k*d*(N-1)*sinθ)]^T。% 参数设置 N 10; % 阵元数 d 0.5; theta_deg linspace(-90, 90, 361); theta deg2rad(theta_deg); % 定义信号方向 theta_s_deg 20; % 期望信号方向 theta_s deg2rad(theta_s_deg); % 定义干扰方向 theta_i_deg [-30, 40]; % 两个干扰方向 theta_i deg2rad(theta_i_deg); % 生成导向矢量函数 steering_vec (ang) exp(1j * 2*pi/lambda * d * (0:N-1)‘ * sin(ang)); a_s steering_vec(theta_s); % 期望信号导向矢量 % 模拟接收数据协方差矩阵 R % 假设期望信号功率为1干扰功率各为100强干扰噪声功率为1单位阵 R a_s * a_s‘; % 信号部分 for i 1:length(theta_i) a_i steering_vec(theta_i(i)); R R 100 * (a_i * a_i‘); % 干扰部分 end R R eye(N); % 加性白噪声部分 % 计算MVDR权向量 w_mvdr (R \ a_s) / (a_s‘ * (R \ a_s)); % 等价于 inv(R)*a_s / (a_s‘*inv(R)*a_s) % 计算MVDR方向图 AF_mvdr zeros(size(theta)); for idx 1:length(theta) a steering_vec(theta(idx)); AF_mvdr(idx) w_mvdr‘ * a; end AF_mvdr_dB 20*log10(abs(AF_mvdr)/max(abs(AF_mvdr))); % 作为对比计算常规波束形成CBF方向图即均匀加权指向期望方向 w_cbf a_s / N; % 常规波束形成权值相位对齐即可 AF_cbf zeros(size(theta)); for idx 1:length(theta) a steering_vec(theta(idx)); AF_cbf(idx) w_cbf‘ * a; end AF_cbf_dB 20*log10(abs(AF_cbf)/max(abs(AF_cbf))); % 绘图 figure; hold on; plot(theta_deg, AF_cbf_dB, ‘b--‘, ‘LineWidth‘, 1.5, ‘DisplayName‘, ‘常规波束形成 (CBF)‘); plot(theta_deg, AF_mvdr_dB, ‘r-‘, ‘LineWidth‘, 2, ‘DisplayName‘, ‘MVDR波束形成‘); % 标记信号和干扰方向 xline(theta_s_deg, ‘k:‘, ‘LineWidth‘, 1, ‘DisplayName‘, [‘信号: ‘, num2str(theta_s_deg), ‘°‘]); for i 1:length(theta_i_deg) xline(theta_i_deg(i), ‘g:‘, ‘LineWidth‘, 1, ‘DisplayName‘, [‘干扰: ‘, num2str(theta_i_deg(i)), ‘°‘]); end xlabel(‘角度 (度)‘); ylabel(‘归一化方向图 (dB)‘); title(‘MVDR与常规波束形成方向图对比‘); legend(‘show‘, ‘Location‘, ‘best‘); grid on; xlim([-90, 90]); ylim([-50, 0]); hold off;从结果图中可以明显看到MVDR算法在期望的20度方向保持了高增益同时在-30度和40度干扰方向形成了非常深的零点低于-40dB而常规波束形成在这两个方向仍有较高的旁瓣响应。这展示了自适应波束形成强大的抗干扰能力。重要提醒与避坑MVDR算法对导向矢量失配即实际的θ_s与预设的有偏差和协方差矩阵R的估计误差非常敏感。在实际应用中R通常由采样数据估计得到R_hat (1/L) * X * X^H其中X是N x L的数据快拍矩阵。如果L不够大经验上需要L 2N估计误差会导致性能严重下降甚至出现信号相消。因此产生了对角加载Diagonal Loading、稳健自适应波束形成等改进算法。在仿真中如果发现MVDR方向图异常如主瓣分裂、增益极低首先检查协方差矩阵是否满秩、条件数是否过大。5. 综合仿真实战一个完整的平面阵设计与分析案例现在我们将前面所有的知识点串联起来完成一个更接近实际工程的仿真案例设计一个均匀矩形平面阵列Uniform Rectangular Array, URA实现波束扫描和方向图赋形。5.1 URA建模与三维方向图绘制一个M x N的URA阵元在x-y平面上等间距排列。其阵列因子是x和y两个维度直线阵列因子的乘积可分离性。% URA参数设置 f 10e9; % 10 GHz典型雷达频段 c 3e8; lambda c / f; dx 0.5 * lambda; % x方向间距 dy 0.5 * lambda; % y方向间距 M 8; % x方向阵元数 N 8; % y方向阵元数 % 定义扫描角度 (方位角az, 俯仰角el) az_scan_deg 30; % 方位扫描30度 el_scan_deg 20; % 俯仰扫描20度 az_scan deg2rad(az_scan_deg); el_scan deg2rad(el_scan_deg); % 生成角度网格 az_deg linspace(-90, 90, 181); el_deg linspace(0, 90, 91); [Az, El] meshgrid(deg2rad(az_deg), deg2rad(el_deg)); % 计算扫描所需的相位补偿 % 对于URA第(m,n)个阵元的相位补偿为 -k * (m*dx*sin(el)*cos(az) n*dy*sin(el)*sin(az)) % 其中 m0:M-1, n0:N-1 [m_idx, n_idx] meshgrid(0:M-1, 0:N-1); phase_comp -2*pi/lambda * (m_idx(:)*dx*sin(el_scan)*cos(az_scan) n_idx(:)*dy*sin(el_scan)*sin(az_scan)); % 初始化方向图矩阵 pattern zeros(size(Az)); % 计算每个角度点的阵列响应 for i 1:size(Az, 1) for j 1:size(Az, 2) az Az(i, j); el El(i, j); % 计算该方向所有阵元的空间相位 spatial_phase 2*pi/lambda * (m_idx(:)*dx*sin(el)*cos(az) n_idx(:)*dy*sin(el)*sin(az)); % 总相位 空间相位 扫描补偿相位 total_phase spatial_phase phase_comp; % 阵列响应等幅激励 pattern(i, j) sum(exp(1j * total_phase)); end end pattern_dB 20*log10(abs(pattern) / max(abs(pattern(:)))); % 三维方向图可视化 figure; surf(az_deg, el_deg, pattern_dB, ‘EdgeColor‘, ‘none‘); xlabel(‘方位角 Azimuth (度)‘); ylabel(‘俯仰角 Elevation (度)‘); zlabel(‘增益 (dB)‘); title([‘URA (8x8) 三维方向图 - 扫描至 (Az‘, num2str(az_scan_deg), ‘°, El‘, num2str(el_scan_deg), ‘°)‘]); colormap(‘jet‘); colorbar; view(45, 30); % 调整视角这段代码生成了一个8x8的URA扫描至(30°, 20°)的三维方向图。你可以通过旋转图形观察主瓣的位置和形状。5.2 方向图赋形实现余割平方波束在某些应用如地面雷达对空搜索时希望波束在俯仰面上按余割平方规律变化使得对等高度、不同距离的目标回波强度一致。这需要对阵列各单元的幅度和相位进行综合优化。这里我们使用一种简单的方法在俯仰维y方向采用特定的幅度分布来近似实现。% 接上部分URA参数现在只关注俯仰面固定方位角az0 el_deg linspace(0, 90, 901); el deg2rad(el_deg); % 目标赋形方向图余割平方从el_min到el_max el_min_deg 5; el_max_deg 60; el_min deg2rad(el_min_deg); el_max deg2rad(el_max_deg); target_gain_dB zeros(size(el)); % 在赋形区间内增益与csc^2(el)成正比即与sin^2(el)成反比 idx_shape (el el_min) (el el_max); target_gain_dB(idx_shape) 20*log10(sin(el_min) ./ sin(el(idx_shape))); % 归一化到el_min处为0dB target_gain_dB(el el_min) 0; % 低于el_min增益快速下降 target_gain_dB(el el_max) target_gain_dB(find(el el_max, 1)); % 高于el_max保持常数 % 将目标增益转换为线性值 target_gain_lin 10.^(target_gain_dB/20); % 使用傅里叶变换法Orchard-Elliott方法简化版综合幅度分布 % 对于直线阵阵列因子AF(theta)是激励电流I_n的离散时间傅里叶变换(DTFT)。 % 我们可以通过逆DTFT即IFFT来近似计算激励。 % 注意这是一个近似方法对于大阵效果较好。 N_y N; % 使用y方向的阵元 % 将目标方向图从角度域转换到u域u sin(el) u sin(el); % 对目标方向图在u域进行插值使其采样点数为阵元数的倍数用于IFFT N_fft 256; % FFT点数远大于阵元数以提高分辨率 u_uniform linspace(-1, 1, N_fft); % u的范围是[-1, 1] % 由于目标方向图只在正半空间定义我们需要构造一个对称的或共轭对称的频域响应。 % 简单起见我们假设为对称实函数。 target_lin_uniform interp1(u, target_gain_lin, abs(u_uniform), ‘linear‘, 0); % 映射到对称区间 % 执行逆FFT得到“激励分布” currents_raw ifft(ifftshift(target_lin_uniform)); % ifftshift将零频移到中心 % 取中间的N_y个点作为阵元激励因为激励是空间有限的 center_idx floor(N_fft/2); start_idx center_idx - floor(N_y/2); end_idx start_idx N_y - 1; I_y currents_raw(start_idx:end_idx); I_y I_y / max(abs(I_y)); % 归一化 % 确保激励为实数幅度分布 I_y abs(I_y); % 绘制综合得到的幅度分布 figure; subplot(1,2,1); stem(0:N_y-1, I_y, ‘filled‘); xlabel(‘阵元序号 (y方向)‘); ylabel(‘归一化激励幅度‘); title(‘综合得到的幅度分布‘); grid on; % 计算采用此幅度分布后的实际方向图 AF_shape zeros(size(el)); for n 0:N_y-1 AF_shape AF_shape I_y(n1) * exp(1j * 2*pi/lambda * n * dy * sin(el)); end AF_shape_dB 20*log10(abs(AF_shape)/max(abs(AF_shape))); % 绘制对比图 subplot(1,2,2); hold on; plot(el_deg, target_gain_dB, ‘b--‘, ‘LineWidth‘, 2, ‘DisplayName‘, ‘目标方向图 (Csc^2)‘); plot(el_deg, AF_shape_dB, ‘r-‘, ‘LineWidth‘, 1.5, ‘DisplayName‘, ‘综合得到的方向图‘); xlabel(‘俯仰角 (度)‘); ylabel(‘归一化增益 (dB)‘); title(‘俯仰面方向图赋形效果‘); legend(‘show‘, ‘Location‘, ‘best‘); grid on; xlim([0, 90]); ylim([-40, 5]); hold off;这个例子展示了通过逆FFT方法进行方向图综合的基本流程。实际工程中会使用更专业的综合算法如凸优化、遗传算法等但原理相通根据期望的方向图反推阵列的激励分布。5.3 性能评估与可视化技巧一个完整的仿真报告离不开性能评估和专业的可视化。方向图切割三维方向图虽全面但二维切割图更利于定量分析。可以固定方位角或俯仰角查看方向图剖面。% 接5.1的三维方向图数据 pattern_dB, az_deg, el_deg fixed_az 30; % 查看方位角30度的俯仰面切割 [~, idx_az] min(abs(az_deg - fixed_az)); pattern_cut_el pattern_dB(:, idx_az); figure; plot(el_deg, pattern_cut_el, ‘LineWidth‘, 2); xlabel(‘俯仰角 (度)‘); ylabel(‘增益 (dB)‘); title([‘URA方向图切割 (方位角固定为 ‘, num2str(fixed_az), ‘°)‘]); grid on;等高线图与二维色彩图用于观察主瓣和旁瓣的二维分布。figure; contour(az_deg, el_deg, pattern_dB, -30:3:0); % 绘制-30dB到0dB的等高线间隔3dB xlabel(‘方位角 (度)‘); ylabel(‘俯仰角 (度)‘); title(‘URA方向图等高线 (扫描波束)‘); colorbar; grid on; figure; imagesc(az_deg, el_deg, pattern_dB); xlabel(‘方位角 (度)‘); ylabel(‘俯仰角 (度)‘); title(‘URA方向图二维色彩图‘); colormap(‘jet‘); colorbar; axis xy; % 确保y轴方向正确关键指标计算将第2.3节的参数提取函数扩展至二维自动计算扫描后波束的指向精度、3dB波束宽度方位和俯仰、峰值旁瓣电平等。项目复盘经验在“gbo.rar”这类仿真项目中代码的模块化至关重要。我将核心功能封装成了独立的函数例如generate_steering_vector(): 生成任意阵列结构的导向矢量。calc_array_pattern(): 计算给定权向量和阵列结构的方向图。beamforming_mvdr(): 实现MVDR等自适应算法。pattern_metrics(): 从方向图数据中提取HPBW, SLL, 指向误差等指标。plot_pattern_2d/3d(): 标准化的绘图函数。 这样主脚本变得非常简洁类似于“搭积木”便于调试和复用。另外务必养成写注释和文档的习惯否则几个月后自己都可能看不懂某段复杂的相位计算是干什么的。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进