简介:面向MATLAB用户与天线设计人员,这份压缩包演示了如何在MATLAB中实现三维极坐标绘图。由于surf、mesh等原生函数基于笛卡尔坐标,难以直接绘制三维方向图,资源以天线三维方向图为例,给出先由极坐标转换至直角坐标、再调用surf渲染的完整流程,并与HFSS仿真结果对比验证,帮助读者掌握等效绘图思路。压缩包仅2个文件,包含1个m脚本与1个csv数据源,合计约320KB,脚本可直接运行并修改,csv存放绘图所需的场强数据。已有16174人学习下载,适合需要绘制球面方向图、功率辐射图等场景的MATLAB中高级用户。通过学习可获得坐标转换与surf绘制的可复用模板,以及从数据导入到图形渲染的完整代码路径,免去自行摸索。1. MATLAB 三维极坐标绘图:为什么方向图不能直接 plot3在 MATLAB 里画天线的三维方向图,多数人第一个念头是把(theta, phi, 增益)这组球坐标数据丢给 plot3,结果画出来要么是两极密、赤道疏的“刺猬”,要么根本看不出主瓣指向。原因在于 plot3 不做球面拓扑,它只按点的顺序连折线;而三维极坐标绘图的本质,是把球坐标数据映射成一张带拓扑关系的曲面。这篇文章用手头就有的 theta/phi 网格数据,按“预处理 → 坐标变换 → surf 绘制 → 信息增强”的顺序,把天线三维方向图做成能交给报告评审的图,最后给出方向图绘制里最常踩的五个坑。适合做天线仿真后处理、链路预算、波束覆盖评估的工程师,也适合拿三维极坐标做课程设计的同学。2. 方向图数据预处理:先弄清 theta/phi 网格和 dB 转线性的次序我拿到方向图数据时从来不急着画,先问三个问题:theta 的起点在哪,phi 的采样是 0~360 还是 -180~180,数据是 dB 还是线性。这三个问题决定后面会不会翻车。方向图数据的来源五花八门,但预处理的目标只有一个:把它们归一成“行对应 theta、列对应 phi”的完整矩阵。这一章把预处理拆成三步,每一步都给出可以直接复制的代码。2.1 三种来源的方向图数据长什么样常见的方向图数据有三种来源。第一种是电磁仿真软件导出的表格,通常每行是(theta, phi, gain_dB),theta 和 phi 按均匀网格遍历;第二种是实测暗室数据,结构类似,但很多系统只给 E 面和 H 面两个主切面,没有完整的三维场数据;第三种是解析表达式,比如偶极子方向图 cos(theta) 的模型,用来验证代码流程最方便。这三种来源的处理路径不同,但最终都要落到同一个接口上:一个 size 为 [numTheta, numPhi] 的数值矩阵。如果手里只有两个主切面的数据,我一般会先说明“三维图只能做近似”,再用旋转对称插值补出完整曲面。仿真导出的数据则要留意导出顺序,这在第 3 章会展开。2.2 meshgrid 的方向:行是 theta,列是 phi生成网格是预处理的第一步,很多人在这里就埋了雷。meshgrid 的返回尺寸和直觉相反,两个输入参数谁先生成行、谁先生成列,直接决定后续 reshape 要不要转置。先用一段代码把网格定下来:% 极角 theta:从 z 轴向下量,0~180 度,换算成弧度 theta = linspace(0, 180, 91) * pi/180; % 方位角 phi:x-y 平面内从 x 轴起,0~360 度 phi = linspace(0, 360, 181) * pi/180; % 生成角网格:meshgrid 的第一个参数是列方向,第二个参数是行方向 [Phi, Theta] = meshgrid(phi, theta); % 验证一下尺寸:行数 = length(theta),列数 = length(phi) fprintf('Theta 矩阵尺寸: %d 行 x %d 列\n', size(Theta));这段代码里,linspace 的采样点数决定了方向图的分辨率:theta 取 91 点意味着每 2 度一个点,phi 取 181 点对应每 2 度一个方位采样。meshgrid(phi, theta) 返回的 Theta 矩阵每一行相同,每一列随 theta 变化,也就是说 theta 沿行方向、phi 沿列方向。这个约定和后面 reshape 仿真数据时必须保持一致。有了网格,就能构造一个用于验证流程的模拟方向图。以下代码生成一个带旁瓣的轴对称波束,主瓣在 z 轴方向,旁瓣在低仰角区:% 模拟方向图:主瓣沿 z 轴,附带三个旁瓣结构 r_lin = abs(cos(Theta)).^1.5 + 0.06 * abs(cos(3 * Theta)); % 转成 dB 相对值;加 eps 防止 log10(0) 产生 -Inf r_dB = 20 * log10(r_lin / max(r_lin(:)) + eps);函数 f = cos(θ)^1.5 描述一个向 z 轴收拢的主波束,指数 1.5 控制主瓣宽度,指数越大主瓣越窄;0.06 * cos(3θ) 在低仰角区制造三个幅度为 -24dB 左右的旁瓣,用来检验第 4 章的动态范围设置。dB 转换里加 eps 是必须的,因为方向图在某些角度上可能是数学上严格的零,直接取对数会出现 -Inf,surf 会在对应位置留下空洞。2.3 仿真导出数据的 reshape 与统一接口实际工作中拿到的往往不是数学函数,而是仿真软件导出的表格。假设导入的表格有三列,依次是 theta(deg)、phi(deg)、gain(dB),预处理代码是这样:% 假设导入的表格三列分别是 theta(deg)、phi(deg)、gain(dB) raw = readmatrix('pattern.csv'); % 提取升序的 theta、phi 列表,并统计各有多少个采样点 tList = unique(raw(:,1)); pList = unique(raw(:,2)); numT = numel(tList); numP = numel(pList); % 按行=theta、列=phi 的结构重排增益数据 gain_dB = reshape(raw(:,3), numT, numP); % dB 转线性幅度:场强方向图用 /20,功率方向图用 /10 gain_lin = 10.^(gain_dB / 20);这里最隐蔽的坑是遍历顺序。如果导出文件是按 theta 外循环、phi 内循环排列的,reshape(gain_dB, numT, numP) 直接得到行对应 theta、列对应 phi 的矩阵;但如果导出顺序是 phi 外循环、theta 内循环,必须改成 reshape(gain_dB, numP, numT).',否则画出来的方向图会发生转置错位,波瓣都跑到了错误的角度上。拿到数据后我一般先对 raw(:,1) 和 raw(:,2) 做 unique,对比两者长度,再打印 reshape 后的矩阵前几行看一眼规律。这个习惯能省掉后面一整轮的排查。另外,dB 转线性时要分清是场强方向图还是功率方向图。场强用 10^(dB/20),功率用 10^(dB/10),用错的话所有数值等比缩小或放大,方向图形状看起来差不多,但颜色映射的尺度就全错了。判断方法很简单:看数据的表头或仿真报告里写的单位是 dB 还是 dBm,如果是 dBm,那通常是功率量纲。3. 球坐标转直角坐标:sph2cart 的约定差与 surf 画法预处理完数据,下一步是把球坐标映射到直角坐标。MATLAB 自带 sph2cart 函数,但它的第二个参数不是极角 theta,而是仰角 elevation。这个差异不报错、不警告,画出来的方向图照样是圆润的曲面,只是波束指向整体偏转,不仔细校验根本发现不了。这一章把两套坐标定义讲透,再给出 surf 的完整画法。3.1 极角与仰角:MATLAB 的 sph2cart 为什么容易用错天线方向图里的 theta 通常采用物理学的极角约定:从 z 轴正方向起算,0 在天顶,pi 在天底。在这个约定下,天顶方向的直角坐标应该是 (0, 0, r)。而 MATLAB 的 sph2cart 采用方位角加仰角的约定,elevation 从 x-y 平面向上量,天顶对应 pi/2。两者相差 90 度,换算关系是:% 极角 theta 转仰角 elevation Elev = pi/2 - Theta; % sph2cart 的第一个参数是方位角 Phi,第二个参数是仰角 Elev [X1, Y1, Z1] = sph2cart(Phi, Elev, r_lin);这段代码里 Phi 是方位角矩阵,Elev 是换算后的仰角矩阵,r_lin 是线性幅度。很多人在这一步直接把 Theta 填进去,结果方向图的“北极”被掰到了 x-y 平面上,原本垂直指向天顶的主波束看起来变成了水平指向。更麻烦的是这个错误不会产生任何报错,图形结构依然完整,只有和已知的方向图基准对照时才能发现。所以我每次处理新数据时,都会先画一个解析偶极子方向图作为基准,确认 z 轴确实是天顶方向,再换真实数据。3.2 手写坐标变换:不依赖 sph2cart 的版本如果不放心 sph2cart 的约定,可以直接用球坐标到直角坐标的标准公式手写变换,这样每个变量的含义都摆在明面上:% 标准数学约定的球坐标变换:theta 从 z 轴起算 X = r_lin .* sin(Theta) .* cos(Phi); Y = r_lin .* sin(Theta) .* sin(Phi); Z = r_lin .* cos(Theta);这个公式里,sin(Theta) 把幅度投影到 x-y 平面,cos(Theta) 决定 z 轴分量;Phi 再决定 x-y 平面内的朝向。从几何上看,它相当于把每一条从原点出发的射线,按照它的 theta、phi 方向拉长到 r_lin 的距离。与 sph2cart 版本相比,手写版更适合检查坐标约定,sph2cart 版更适合从现成代码里搬用,两个我都放在自己的工具箱里,画完互相验证一次。3.3 surf 绘制与 edgealpha 参数坐标变换完成后,用 surf 绘制球面。这一步的关键参数是 EdgeAlpha 和 axis equal,少了任何一个,图形都会失真:figure; surf(X, Y, Z, r_lin, 'EdgeAlpha', 0.25, 'FaceAlpha', 1.0); shading interp; % 颜色在网格间做插值,避免色块边界 axis equal; % 关键:不加这一步,z 轴会被拉伸成椭圆 view(3); % 切换到三维视角 colormap(jet); colorbar; xlabel('x'); ylabel('y'); zlabel('z'); title('天线三维方向图 (线性幅度)');surf 的前三个参数是 X、Y、Z 坐标矩阵,第四个参数 r_lin 作为 CData 控制颜色映射。EdgeAlpha 控制网格线的透明度,0.25 是一个常用取值,既能保留曲面结构又不会让网格线抢掉主视觉;想要纯色曲面可以把 EdgeAlpha 设为 'none',但代价是近距离观察时不容易分辨曲面起伏。shading interp 让颜色在网格点之间平滑过渡,避免亮暗分明的色块感。axis equal 在这里极其关键,不加的话 MATLAB 会自动拉伸坐标轴,一个标准球体会被压成椭球,主瓣宽度看起来全错了。4. 让方向图可读:颜色映射、光照与主平面切片surf 画出来的裸图能看,但离“能进报告”还差几步。这一章处理三个问题:用线性幅度还是 dB 幅度、光照怎么打、怎么把 E 面和 H 面切片叠进三维图。4.1 线性幅度还是 dB:两张图表达的信息不同用线性幅度直接画 surf,主瓣非常突出,但旁瓣会被压到几乎贴地,颜色映射的整个动态范围都被主瓣占据;用 dB 值画,旁瓣结构能清晰显现,但主瓣看起来相对“瘦”。我的做法是:看波束覆盖和主瓣形状用线性,评估旁瓣电平、前后比时用 dB,并在 dB 图上做动态范围截断:% 用 dB 值绘制并截断动态范围:从峰值往下 30dB surf(X, Y, Z, r_dB, 'EdgeAlpha', 0.25); clim([max(r_dB(:)) - 30, max(r_dB(:))]); colormap(parula); colorbar;clim 把颜色映射的下限锁定在“峰值减 30dB”,低于这个值的区域全部用色标最底端的颜色表示。天线工程里 30dB、40dB、50dB 是几种常用的截断范围,看蜂窝基站方向图时我习惯用 30dB,因为它能同时展示主瓣和第一旁瓣的关系,又不至于把远旁瓣的细节过度放大。需要注意的是,旧版 MATLAB 里这个函数叫 caxis,新版统一为 clim;如果代码要在不同版本间迁移,建议加个判断或直接升级到支持 clim 的版本。4.2 光照与材质:让曲面有立体感surf 默认的平面着色在三维视图里看起来非常平,尤其是大角度旋转视角时,很难分辨曲面的凹凸。加几行光照命令,立体感立刻出来:camlight('headlight'); % 在相机位置加一个点光源 lighting gouraud; % 对颜色做插值光照,曲面更平滑 material dull; % 使用漫反射材质,避免高光掩盖数据camlight('headlight') 把光源放在观察者身后,保证照亮的部分永远正对着视角;lighting gouraud 对每个顶点的法线做插值,比 lighting flat 的块状着色平滑得多;material dull 降低镜面反射,防止方向图表面出现一块块白色高光。这里需要注意的是,光照命令必须在 surf 之后执行,否则不会生效;而且如果后面还要叠加切片曲线,光照可能会改变曲线的观感,所以我会先画切片再打光。4.3 把 E 面和 H 面切片叠进三维图纯三维曲面在报告里虽然直观,但工程师评审时往往还想看到两个主平面的精确轮廓。把 E 面和 H 面切片叠加到三维图上,一张图就能同时表达三维波束和二维主平面数据:hold on; % E 面:phi = 0°,对应矩阵第一列 idxE = 1; plot3(X(:,idxE), Y(:,idxE), Z(:,idxE), 'k-', 'LineWidth', 2); % H 面:phi = 90°,在 0~360 采样下位于行数的四分之一处 idxH = round(numel(phi) / 4); plot3(X(:,idxH), Y(:,idxH), Z(:,idxH), 'w-', 'LineWidth', 1.5); hold off;这条切片线的原理是沿着球面在指定 phi 角度上“切一刀”,把对应列的 (X, Y, Z) 坐标连成一条三维曲线。idxE 取第一列对应 phi = 0°,idxH 取四分之一列数对应 phi = 90°,前提是 phi 从 0 开始均匀采样到 360。如果数据里 phi 的取值范围是 -180~180,那么 90° 对应的就不是四分之一处,需要先用 find 找到 phi 列表中离 pi/2 最近的索引。黑色粗线标 E 面、白色细线标 H 面,在 jet 或 parula 色标上都有足够的对比度。5. 避坑排查:五个方向图绘制的翻车现场前面几章把完整流程走通了,但真正到了自己动手的时候,翻车点往往集中在几个固定的位置。这一章写五个我实际遇到过的问题,每条按“现象 → 原因 → 解决”的顺序记录,省得你再踩一遍。5.1 翻车一:方向图整体旋转了 90 度现象:标准偶极子的方向图应该是一个沿 z 轴竖起来的“甜甜圈”,画出来却变成躺在 x-y 平面上的形态,主瓣方向指向了水平方向。原因:sph2cart 的第二个参数要求传仰角 elevation,从 x-y 平面起算;而方向图数据里的 theta 是从 z 轴起算的极角。直接把 Theta 传给 sph2cart,等于把天顶方向硬掰到了水平面上,而且整个过程不报错。解决:传入之前先做一次换算,Elev = pi/2 - Theta。更稳妥的做法是用手写坐标变换 X = r.*sin(Theta).*cos(Phi) 那套公式,两套代码各跑一遍,对比峰值坐标是否一致。现在我每画一张方向图,都会先用解析偶极子跑一次基准校验,再上真实数据。5.2 翻车二:方向图变成刺猬或出现空洞现象:曲面布满从原点向外散射的尖刺,或者某些角度区域出现 NaN 空白,图形像被虫啃过。原因:实测数据或某些仿真导出数据并不是严格等间隔网格。surf 假设输入是规则网格,一旦 theta 或 phi 的采样间隔不均匀,网格点之间的连线会交叉,产生刺猬状伪影;数据里有 NaN 时则直接形成空洞。解决:先用 scatteredInterpolant 把散点插值到均匀网格上:% thetaVec、phiVec、rVec 是从数据表提取的三列 F = scatteredInterpolant(thetaVec, phiVec, rVec, 'linear', 'none'); % 在规则网格上采样 [PhiG, ThetaG] = meshgrid(phi, theta); rInterp = F(ThetaG(:), PhiG(:)); rInterp = reshape(rInterp, size(ThetaG));scatteredInterpolant 的 'linear' 表示线性插值,适合方向图这种变化平缓的数据;'none' 表示不对网格范围外的点做外插,避免边缘出现无法解释的数值。插值完成后还要检查是否还有 NaN,如果存在,多半是数据覆盖范围不够。5.3 翻车三:颜色和数值对不上现象:曲面形状看起来正常,但颜色表达的信息和实际数值不一致,比如 colorbar 上标的是负的 dB 值,曲面中心却显示为最亮的黄色。原因:surf 的 CData 参数传错了矩阵。很多人把坐标变换后的 Z 矩阵当作 CData 传入,导致颜色表现的是 z 轴高度而不是方向图幅度;或者内部做了线性到 dB 的转换,但忘了同步更新 colorbar 的量纲。解决:显式给 surf 传入需要显示的幅度矩阵,并且保证这个矩阵和坐标变换用的是同一份数据;colorbar 上额外标注单位。如果截断了动态范围,用 clim 控制上下限,避免颜色分布被远旁瓣带偏。绘制完成后我一般会抽查几个特定角度的数值,和 colorbar 上的读数对照一遍。5.4 翻车四:网格线把表面全盖住现象:图是画出来了,但打印到报告里黑乎乎一片,主瓣轮廓完全被网格线遮住,根本看不清波束形状。原因:surf 默认的 EdgeColor 是黑色,并且边缘透明度为 1,网格一多整个曲面就像罩了一张黑网。解决:设置 EdgeAlpha 为 0.15~0.3,或干脆设成 'none' 关闭网格线。如果还需要突出曲面起伏,可以用 contourslice 在指定平面上画等高线,而不是依赖网格线。上面第 3 章完整代码里的 0.25 就是我常用的值,网格密度高的时候可以再降一档。5.5 翻车五:第三方极坐标绘图脚本不生效现象:从网上找到的 polarplot3d 脚本,第一参数传入方向图矩阵后报维度错误,或者画出来图形方向不对。原因:第三方三维极坐标绘图函数内部定义的坐标约定和天线方向图的常用约定往往不一致,有的把 0 度放在 x 轴方向,有的把极角从 x-y 平面起算,和你的 theta/phi 定义差出 90 度或 180 度。解决:既然要画三维极坐标方向图,干脆用 sph2cart 加 surf 自己搭一个十行代码的小模板,坐标约定自己完全掌控。第三方脚本通常只适合快速预览,不适合直接进入正式报告流程。网上脚本再有名的,也要先用偶极子基准函数验证坐标方向再上真实数据。6. 进阶技巧:模拟下倾角并验证波束指向实际工程里,评估基站天线覆盖时经常要回答一个问题:把方向图往下倾 10 度,主瓣扫到哪里?与其在仿真软件里重算一遍,不如对已经画好的球面做一次刚体旋转变换,几行代码就能看效果。6.1 旋转矩阵与下倾视角把方向图绕 y 轴旋转一个角度,即可模拟电下倾或机械下倾的效果:% 下倾角 10 度,转弧度 beta = 10 * pi/180; % 绕 y 轴旋转的旋转矩阵 Ry = [cos(beta) 0 sin(beta); 0 1 0 ; -sin(beta) 0 cos(beta)]; % 对每个散点做旋转:行向量右乘旋转矩阵 XYZ = [X(:), Y(:), Z(:)] * Ry.'; Xr = reshape(XYZ(:,1), size(X)); Yr = reshape(XYZ(:,2), size(Y)); Zr = reshape(XYZ(:,3), size(Z));旋转矩阵里,beta 为正时 z 轴分量减小、x 轴分量增大,对应波束向下倾;如果发现方向反了,把 beta 取负即可。这里的“下倾”是纯几何旋转,适合快速评估,如果要做带地面反射的精确覆盖仿真,还是要回到电磁仿真软件里重新求解。6.2 用峰值索引自检波束指向做完旋转变换,要验证主瓣确实落到了预期的角度上。最便宜的方式是找幅度最大点对应的 theta 和 phi:% 在旋转后的几何里找到最高的点的索引 [rMax, idx] = max(Zr(:)); % 一维索引转回矩阵行列 [ti, pi_] = ind2sub(size(Zr), idx); % 反查角度并转成度数 thetaPeak = theta(ti) * 180/pi; phiPeak = phi(pi_) * 180/pi; fprintf('波束峰值: theta=%.1f deg, phi=%.1f deg\n', thetaPeak, phiPeak);这段代码把旋转后的曲面里最高的点找出来,反推它在原始 theta/phi 网格里的位置。理想情况下,下倾 10 度后峰值 theta 应该从 0 度变到 10 度;如果反查出来的角度对不上,说明旋转矩阵的方向或者索引和网格的对应关系有误。这个自检步骤看起来简单,但能挡住绝大多数坐标约定问题。从那以后,我每次画三维方向图都会先跑一遍解析偶极子做基准校验,确认坐标约定没有反转,再上真实数据;旋转变换做完也一定用峰值反查确认指向。这个习惯帮我省掉了好几次报告返工,希望帮到你。本文还有配套的精品资源,点击获取
简介:面向MATLAB用户与天线设计人员,这份压缩包演示了如何在MATLAB中实现三维极坐标绘图。由于surf、mesh等原生函数基于笛卡尔坐标,难以直接绘制三维方向图,资源以天线三维方向图为例,给出先由极坐标转换至直角坐标、再调用surf渲染的完整流程,并与HFSS仿真结果对比验证,帮助读者掌握等效绘图思路。压缩包仅2个文件,包含1个m脚本与1个csv数据源,合计约320KB,脚本可直接运行并修改,csv存放绘图所需的场强数据。已有16174人学习下载,适合需要绘制球面方向图、功率辐射图等场景的MATLAB中高级用户。通过学习可获得坐标转换与surf绘制的可复用模板,以及从数据导入到图形渲染的完整代码路径,免去自行摸索。1. MATLAB 三维极坐标绘图:为什么方向图不能直接 plot3在 MATLAB 里画天线的三维方向图,多数人第一个念头是把(theta, phi, 增益)这组球坐标数据丢给 plot3,结果画出来要么是两极密、赤道疏的“刺猬”,要么根本看不出主瓣指向。原因在于 plot3 不做球面拓扑,它只按点的顺序连折线;而三维极坐标绘图的本质,是把球坐标数据映射成一张带拓扑关系的曲面。这篇文章用手头就有的 theta/phi 网格数据,按“预处理 → 坐标变换 → surf 绘制 → 信息增强”的顺序,把天线三维方向图做成能交给报告评审的图,最后给出方向图绘制里最常踩的五个坑。适合做天线仿真后处理、链路预算、波束覆盖评估的工程师,也适合拿三维极坐标做课程设计的同学。2. 方向图数据预处理:先弄清 theta/phi 网格和 dB 转线性的次序我拿到方向图数据时从来不急着画,先问三个问题:theta 的起点在哪,phi 的采样是 0~360 还是 -180~180,数据是 dB 还是线性。这三个问题决定后面会不会翻车。方向图数据的来源五花八门,但预处理的目标只有一个:把它们归一成“行对应 theta、列对应 phi”的完整矩阵。这一章把预处理拆成三步,每一步都给出可以直接复制的代码。2.1 三种来源的方向图数据长什么样常见的方向图数据有三种来源。第一种是电磁仿真软件导出的表格,通常每行是(theta, phi, gain_dB),theta 和 phi 按均匀网格遍历;第二种是实测暗室数据,结构类似,但很多系统只给 E 面和 H 面两个主切面,没有完整的三维场数据;第三种是解析表达式,比如偶极子方向图 cos(theta) 的模型,用来验证代码流程最方便。这三种来源的处理路径不同,但最终都要落到同一个接口上:一个 size 为 [numTheta, numPhi] 的数值矩阵。如果手里只有两个主切面的数据,我一般会先说明“三维图只能做近似”,再用旋转对称插值补出完整曲面。仿真导出的数据则要留意导出顺序,这在第 3 章会展开。2.2 meshgrid 的方向:行是 theta,列是 phi生成网格是预处理的第一步,很多人在这里就埋了雷。meshgrid 的返回尺寸和直觉相反,两个输入参数谁先生成行、谁先生成列,直接决定后续 reshape 要不要转置。先用一段代码把网格定下来:% 极角 theta:从 z 轴向下量,0~180 度,换算成弧度 theta = linspace(0, 180, 91) * pi/180; % 方位角 phi:x-y 平面内从 x 轴起,0~360 度 phi = linspace(0, 360, 181) * pi/180; % 生成角网格:meshgrid 的第一个参数是列方向,第二个参数是行方向 [Phi, Theta] = meshgrid(phi, theta); % 验证一下尺寸:行数 = length(theta),列数 = length(phi) fprintf('Theta 矩阵尺寸: %d 行 x %d 列\n', size(Theta));这段代码里,linspace 的采样点数决定了方向图的分辨率:theta 取 91 点意味着每 2 度一个点,phi 取 181 点对应每 2 度一个方位采样。meshgrid(phi, theta) 返回的 Theta 矩阵每一行相同,每一列随 theta 变化,也就是说 theta 沿行方向、phi 沿列方向。这个约定和后面 reshape 仿真数据时必须保持一致。有了网格,就能构造一个用于验证流程的模拟方向图。以下代码生成一个带旁瓣的轴对称波束,主瓣在 z 轴方向,旁瓣在低仰角区:% 模拟方向图:主瓣沿 z 轴,附带三个旁瓣结构 r_lin = abs(cos(Theta)).^1.5 + 0.06 * abs(cos(3 * Theta)); % 转成 dB 相对值;加 eps 防止 log10(0) 产生 -Inf r_dB = 20 * log10(r_lin / max(r_lin(:)) + eps);函数 f = cos(θ)^1.5 描述一个向 z 轴收拢的主波束,指数 1.5 控制主瓣宽度,指数越大主瓣越窄;0.06 * cos(3θ) 在低仰角区制造三个幅度为 -24dB 左右的旁瓣,用来检验第 4 章的动态范围设置。dB 转换里加 eps 是必须的,因为方向图在某些角度上可能是数学上严格的零,直接取对数会出现 -Inf,surf 会在对应位置留下空洞。2.3 仿真导出数据的 reshape 与统一接口实际工作中拿到的往往不是数学函数,而是仿真软件导出的表格。假设导入的表格有三列,依次是 theta(deg)、phi(deg)、gain(dB),预处理代码是这样:% 假设导入的表格三列分别是 theta(deg)、phi(deg)、gain(dB) raw = readmatrix('pattern.csv'); % 提取升序的 theta、phi 列表,并统计各有多少个采样点 tList = unique(raw(:,1)); pList = unique(raw(:,2)); numT = numel(tList); numP = numel(pList); % 按行=theta、列=phi 的结构重排增益数据 gain_dB = reshape(raw(:,3), numT, numP); % dB 转线性幅度:场强方向图用 /20,功率方向图用 /10 gain_lin = 10.^(gain_dB / 20);这里最隐蔽的坑是遍历顺序。如果导出文件是按 theta 外循环、phi 内循环排列的,reshape(gain_dB, numT, numP) 直接得到行对应 theta、列对应 phi 的矩阵;但如果导出顺序是 phi 外循环、theta 内循环,必须改成 reshape(gain_dB, numP, numT).',否则画出来的方向图会发生转置错位,波瓣都跑到了错误的角度上。拿到数据后我一般先对 raw(:,1) 和 raw(:,2) 做 unique,对比两者长度,再打印 reshape 后的矩阵前几行看一眼规律。这个习惯能省掉后面一整轮的排查。另外,dB 转线性时要分清是场强方向图还是功率方向图。场强用 10^(dB/20),功率用 10^(dB/10),用错的话所有数值等比缩小或放大,方向图形状看起来差不多,但颜色映射的尺度就全错了。判断方法很简单:看数据的表头或仿真报告里写的单位是 dB 还是 dBm,如果是 dBm,那通常是功率量纲。3. 球坐标转直角坐标:sph2cart 的约定差与 surf 画法预处理完数据,下一步是把球坐标映射到直角坐标。MATLAB 自带 sph2cart 函数,但它的第二个参数不是极角 theta,而是仰角 elevation。这个差异不报错、不警告,画出来的方向图照样是圆润的曲面,只是波束指向整体偏转,不仔细校验根本发现不了。这一章把两套坐标定义讲透,再给出 surf 的完整画法。3.1 极角与仰角:MATLAB 的 sph2cart 为什么容易用错天线方向图里的 theta 通常采用物理学的极角约定:从 z 轴正方向起算,0 在天顶,pi 在天底。在这个约定下,天顶方向的直角坐标应该是 (0, 0, r)。而 MATLAB 的 sph2cart 采用方位角加仰角的约定,elevation 从 x-y 平面向上量,天顶对应 pi/2。两者相差 90 度,换算关系是:% 极角 theta 转仰角 elevation Elev = pi/2 - Theta; % sph2cart 的第一个参数是方位角 Phi,第二个参数是仰角 Elev [X1, Y1, Z1] = sph2cart(Phi, Elev, r_lin);这段代码里 Phi 是方位角矩阵,Elev 是换算后的仰角矩阵,r_lin 是线性幅度。很多人在这一步直接把 Theta 填进去,结果方向图的“北极”被掰到了 x-y 平面上,原本垂直指向天顶的主波束看起来变成了水平指向。更麻烦的是这个错误不会产生任何报错,图形结构依然完整,只有和已知的方向图基准对照时才能发现。所以我每次处理新数据时,都会先画一个解析偶极子方向图作为基准,确认 z 轴确实是天顶方向,再换真实数据。3.2 手写坐标变换:不依赖 sph2cart 的版本如果不放心 sph2cart 的约定,可以直接用球坐标到直角坐标的标准公式手写变换,这样每个变量的含义都摆在明面上:% 标准数学约定的球坐标变换:theta 从 z 轴起算 X = r_lin .* sin(Theta) .* cos(Phi); Y = r_lin .* sin(Theta) .* sin(Phi); Z = r_lin .* cos(Theta);这个公式里,sin(Theta) 把幅度投影到 x-y 平面,cos(Theta) 决定 z 轴分量;Phi 再决定 x-y 平面内的朝向。从几何上看,它相当于把每一条从原点出发的射线,按照它的 theta、phi 方向拉长到 r_lin 的距离。与 sph2cart 版本相比,手写版更适合检查坐标约定,sph2cart 版更适合从现成代码里搬用,两个我都放在自己的工具箱里,画完互相验证一次。3.3 surf 绘制与 edgealpha 参数坐标变换完成后,用 surf 绘制球面。这一步的关键参数是 EdgeAlpha 和 axis equal,少了任何一个,图形都会失真:figure; surf(X, Y, Z, r_lin, 'EdgeAlpha', 0.25, 'FaceAlpha', 1.0); shading interp; % 颜色在网格间做插值,避免色块边界 axis equal; % 关键:不加这一步,z 轴会被拉伸成椭圆 view(3); % 切换到三维视角 colormap(jet); colorbar; xlabel('x'); ylabel('y'); zlabel('z'); title('天线三维方向图 (线性幅度)');surf 的前三个参数是 X、Y、Z 坐标矩阵,第四个参数 r_lin 作为 CData 控制颜色映射。EdgeAlpha 控制网格线的透明度,0.25 是一个常用取值,既能保留曲面结构又不会让网格线抢掉主视觉;想要纯色曲面可以把 EdgeAlpha 设为 'none',但代价是近距离观察时不容易分辨曲面起伏。shading interp 让颜色在网格点之间平滑过渡,避免亮暗分明的色块感。axis equal 在这里极其关键,不加的话 MATLAB 会自动拉伸坐标轴,一个标准球体会被压成椭球,主瓣宽度看起来全错了。4. 让方向图可读:颜色映射、光照与主平面切片surf 画出来的裸图能看,但离“能进报告”还差几步。这一章处理三个问题:用线性幅度还是 dB 幅度、光照怎么打、怎么把 E 面和 H 面切片叠进三维图。4.1 线性幅度还是 dB:两张图表达的信息不同用线性幅度直接画 surf,主瓣非常突出,但旁瓣会被压到几乎贴地,颜色映射的整个动态范围都被主瓣占据;用 dB 值画,旁瓣结构能清晰显现,但主瓣看起来相对“瘦”。我的做法是:看波束覆盖和主瓣形状用线性,评估旁瓣电平、前后比时用 dB,并在 dB 图上做动态范围截断:% 用 dB 值绘制并截断动态范围:从峰值往下 30dB surf(X, Y, Z, r_dB, 'EdgeAlpha', 0.25); clim([max(r_dB(:)) - 30, max(r_dB(:))]); colormap(parula); colorbar;clim 把颜色映射的下限锁定在“峰值减 30dB”,低于这个值的区域全部用色标最底端的颜色表示。天线工程里 30dB、40dB、50dB 是几种常用的截断范围,看蜂窝基站方向图时我习惯用 30dB,因为它能同时展示主瓣和第一旁瓣的关系,又不至于把远旁瓣的细节过度放大。需要注意的是,旧版 MATLAB 里这个函数叫 caxis,新版统一为 clim;如果代码要在不同版本间迁移,建议加个判断或直接升级到支持 clim 的版本。4.2 光照与材质:让曲面有立体感surf 默认的平面着色在三维视图里看起来非常平,尤其是大角度旋转视角时,很难分辨曲面的凹凸。加几行光照命令,立体感立刻出来:camlight('headlight'); % 在相机位置加一个点光源 lighting gouraud; % 对颜色做插值光照,曲面更平滑 material dull; % 使用漫反射材质,避免高光掩盖数据camlight('headlight') 把光源放在观察者身后,保证照亮的部分永远正对着视角;lighting gouraud 对每个顶点的法线做插值,比 lighting flat 的块状着色平滑得多;material dull 降低镜面反射,防止方向图表面出现一块块白色高光。这里需要注意的是,光照命令必须在 surf 之后执行,否则不会生效;而且如果后面还要叠加切片曲线,光照可能会改变曲线的观感,所以我会先画切片再打光。4.3 把 E 面和 H 面切片叠进三维图纯三维曲面在报告里虽然直观,但工程师评审时往往还想看到两个主平面的精确轮廓。把 E 面和 H 面切片叠加到三维图上,一张图就能同时表达三维波束和二维主平面数据:hold on; % E 面:phi = 0°,对应矩阵第一列 idxE = 1; plot3(X(:,idxE), Y(:,idxE), Z(:,idxE), 'k-', 'LineWidth', 2); % H 面:phi = 90°,在 0~360 采样下位于行数的四分之一处 idxH = round(numel(phi) / 4); plot3(X(:,idxH), Y(:,idxH), Z(:,idxH), 'w-', 'LineWidth', 1.5); hold off;这条切片线的原理是沿着球面在指定 phi 角度上“切一刀”,把对应列的 (X, Y, Z) 坐标连成一条三维曲线。idxE 取第一列对应 phi = 0°,idxH 取四分之一列数对应 phi = 90°,前提是 phi 从 0 开始均匀采样到 360。如果数据里 phi 的取值范围是 -180~180,那么 90° 对应的就不是四分之一处,需要先用 find 找到 phi 列表中离 pi/2 最近的索引。黑色粗线标 E 面、白色细线标 H 面,在 jet 或 parula 色标上都有足够的对比度。5. 避坑排查:五个方向图绘制的翻车现场前面几章把完整流程走通了,但真正到了自己动手的时候,翻车点往往集中在几个固定的位置。这一章写五个我实际遇到过的问题,每条按“现象 → 原因 → 解决”的顺序记录,省得你再踩一遍。5.1 翻车一:方向图整体旋转了 90 度现象:标准偶极子的方向图应该是一个沿 z 轴竖起来的“甜甜圈”,画出来却变成躺在 x-y 平面上的形态,主瓣方向指向了水平方向。原因:sph2cart 的第二个参数要求传仰角 elevation,从 x-y 平面起算;而方向图数据里的 theta 是从 z 轴起算的极角。直接把 Theta 传给 sph2cart,等于把天顶方向硬掰到了水平面上,而且整个过程不报错。解决:传入之前先做一次换算,Elev = pi/2 - Theta。更稳妥的做法是用手写坐标变换 X = r.*sin(Theta).*cos(Phi) 那套公式,两套代码各跑一遍,对比峰值坐标是否一致。现在我每画一张方向图,都会先用解析偶极子跑一次基准校验,再上真实数据。5.2 翻车二:方向图变成刺猬或出现空洞现象:曲面布满从原点向外散射的尖刺,或者某些角度区域出现 NaN 空白,图形像被虫啃过。原因:实测数据或某些仿真导出数据并不是严格等间隔网格。surf 假设输入是规则网格,一旦 theta 或 phi 的采样间隔不均匀,网格点之间的连线会交叉,产生刺猬状伪影;数据里有 NaN 时则直接形成空洞。解决:先用 scatteredInterpolant 把散点插值到均匀网格上:% thetaVec、phiVec、rVec 是从数据表提取的三列 F = scatteredInterpolant(thetaVec, phiVec, rVec, 'linear', 'none'); % 在规则网格上采样 [PhiG, ThetaG] = meshgrid(phi, theta); rInterp = F(ThetaG(:), PhiG(:)); rInterp = reshape(rInterp, size(ThetaG));scatteredInterpolant 的 'linear' 表示线性插值,适合方向图这种变化平缓的数据;'none' 表示不对网格范围外的点做外插,避免边缘出现无法解释的数值。插值完成后还要检查是否还有 NaN,如果存在,多半是数据覆盖范围不够。5.3 翻车三:颜色和数值对不上现象:曲面形状看起来正常,但颜色表达的信息和实际数值不一致,比如 colorbar 上标的是负的 dB 值,曲面中心却显示为最亮的黄色。原因:surf 的 CData 参数传错了矩阵。很多人把坐标变换后的 Z 矩阵当作 CData 传入,导致颜色表现的是 z 轴高度而不是方向图幅度;或者内部做了线性到 dB 的转换,但忘了同步更新 colorbar 的量纲。解决:显式给 surf 传入需要显示的幅度矩阵,并且保证这个矩阵和坐标变换用的是同一份数据;colorbar 上额外标注单位。如果截断了动态范围,用 clim 控制上下限,避免颜色分布被远旁瓣带偏。绘制完成后我一般会抽查几个特定角度的数值,和 colorbar 上的读数对照一遍。5.4 翻车四:网格线把表面全盖住现象:图是画出来了,但打印到报告里黑乎乎一片,主瓣轮廓完全被网格线遮住,根本看不清波束形状。原因:surf 默认的 EdgeColor 是黑色,并且边缘透明度为 1,网格一多整个曲面就像罩了一张黑网。解决:设置 EdgeAlpha 为 0.15~0.3,或干脆设成 'none' 关闭网格线。如果还需要突出曲面起伏,可以用 contourslice 在指定平面上画等高线,而不是依赖网格线。上面第 3 章完整代码里的 0.25 就是我常用的值,网格密度高的时候可以再降一档。5.5 翻车五:第三方极坐标绘图脚本不生效现象:从网上找到的 polarplot3d 脚本,第一参数传入方向图矩阵后报维度错误,或者画出来图形方向不对。原因:第三方三维极坐标绘图函数内部定义的坐标约定和天线方向图的常用约定往往不一致,有的把 0 度放在 x 轴方向,有的把极角从 x-y 平面起算,和你的 theta/phi 定义差出 90 度或 180 度。解决:既然要画三维极坐标方向图,干脆用 sph2cart 加 surf 自己搭一个十行代码的小模板,坐标约定自己完全掌控。第三方脚本通常只适合快速预览,不适合直接进入正式报告流程。网上脚本再有名的,也要先用偶极子基准函数验证坐标方向再上真实数据。6. 进阶技巧:模拟下倾角并验证波束指向实际工程里,评估基站天线覆盖时经常要回答一个问题:把方向图往下倾 10 度,主瓣扫到哪里?与其在仿真软件里重算一遍,不如对已经画好的球面做一次刚体旋转变换,几行代码就能看效果。6.1 旋转矩阵与下倾视角把方向图绕 y 轴旋转一个角度,即可模拟电下倾或机械下倾的效果:% 下倾角 10 度,转弧度 beta = 10 * pi/180; % 绕 y 轴旋转的旋转矩阵 Ry = [cos(beta) 0 sin(beta); 0 1 0 ; -sin(beta) 0 cos(beta)]; % 对每个散点做旋转:行向量右乘旋转矩阵 XYZ = [X(:), Y(:), Z(:)] * Ry.'; Xr = reshape(XYZ(:,1), size(X)); Yr = reshape(XYZ(:,2), size(Y)); Zr = reshape(XYZ(:,3), size(Z));旋转矩阵里,beta 为正时 z 轴分量减小、x 轴分量增大,对应波束向下倾;如果发现方向反了,把 beta 取负即可。这里的“下倾”是纯几何旋转,适合快速评估,如果要做带地面反射的精确覆盖仿真,还是要回到电磁仿真软件里重新求解。6.2 用峰值索引自检波束指向做完旋转变换,要验证主瓣确实落到了预期的角度上。最便宜的方式是找幅度最大点对应的 theta 和 phi:% 在旋转后的几何里找到最高的点的索引 [rMax, idx] = max(Zr(:)); % 一维索引转回矩阵行列 [ti, pi_] = ind2sub(size(Zr), idx); % 反查角度并转成度数 thetaPeak = theta(ti) * 180/pi; phiPeak = phi(pi_) * 180/pi; fprintf('波束峰值: theta=%.1f deg, phi=%.1f deg\n', thetaPeak, phiPeak);这段代码把旋转后的曲面里最高的点找出来,反推它在原始 theta/phi 网格里的位置。理想情况下,下倾 10 度后峰值 theta 应该从 0 度变到 10 度;如果反查出来的角度对不上,说明旋转矩阵的方向或者索引和网格的对应关系有误。这个自检步骤看起来简单,但能挡住绝大多数坐标约定问题。从那以后,我每次画三维方向图都会先跑一遍解析偶极子做基准校验,确认坐标约定没有反转,再上真实数据;旋转变换做完也一定用峰值反查确认指向。这个习惯帮我省掉了好几次报告返工,希望帮到你。本文还有配套的精品资源,点击获取