道尔夫-切比雪夫加权MATLAB实现:旁瓣抑制与主瓣展宽权衡
发布时间:2026/10/24 12:21:35 锦皓数字建站

简介面向阵列天线与阵列信号处理方向的学习者这份资源提供基于 MATLAB 的线阵道尔夫—切比雪夫加权综合实现用于解决在给定副瓣电平约束下求解各阵元激励幅度、并观察方向图变化的问题适合已具备天线基础与 MATLAB 编程能力、需要动手验证加权效果的本科生、研究生及工程人员。压缩包共 3 个文件以 2 个 m 脚本和 1 个 txt 说明为主脚本分别完成切比雪夫激励幅度计算与加权后方向图仿真说明文件交代使用方式包体约 2KB轻量易读。代码支持自行修改阵元个数、阵元间距、副瓣电平、波束指向角度与信号频率通过切换代码块即可对比不同参数下激励分布与方向图主瓣宽度、副瓣高度的差异关键步骤均附中文注释便于逐行理解综合流程。目前已有 1003 人学习下载可作为课堂实验、课程设计或自学阵列综合的参考脚本。1. 旁瓣压到 -30 dB 之后主瓣胖了一圈道尔夫-切比雪夫加权在做什么用 MATLAB 给线阵天线做道尔夫-切比雪夫加权第一次把 16 元阵的旁瓣从均匀加权的 -13.2 dB 硬压到 -30 dB 时方向图确实漂亮了可 -3 dB 波束宽度从 6.4° 涨到 9.1°测角精度掉了一截。这不是代码写错了而是这套加权方式的固有代价在阵元数和阵元间距给定的前提下它给出的是“等旁瓣电平约束下主瓣最窄”的那唯一一组电流分布。换句话说报一个旁瓣指标它就反算出唯一的激励系数没有第二组解可以挑选。这套 MATLAB 实现要落地的就是把这句话变成能跑的代码输入阵元数、间距、目标旁瓣电平输出加权系数、归一化方向图、实测旁瓣电平、-3 dB 波束宽度和零点位置。做阵列综合的、写雷达链路预算的、要给 DBF 通道加窗的都能直接把系数拿去用。2. 把旁瓣电平翻译成多项式系数x0、阶数与主瓣展宽的三个公式2.1 阵因子为什么能写成切比雪夫多项式N 元等间距线阵阵元位置取以阵列中心为原点第 n 个阵元的坐标是 x_n (n-(N-1)/2)dn 0,1,…,N-1。激励电流记作 I_n阵因子就是各阵元贡献的相干叠加AF(θ) Σ_{n0}^{N-1} I_n · exp(j·k·x_n·sinθ)把空间相位统一记成 ψ kd·sinθ 2π(d/λ)·sinθ上式可以写成 AF Σ I_n exp(j(n-(N-1)/2)ψ)。激励对称时 I_n I_{N-1-n}配对之后只剩余弦项AF(ψ) I_c 2Σ_{r0} I_r cos(rψ)对奇数 Nr 取 0,1,…,(N-1)/2对偶数 Nr 取 1/2,3/2,…,(N-1)/2因为阵元落在半整数倍的位置上。关键一步是变量替换 y cos(ψ/2)。代入 cos(rψ) 并用倍角公式展开AF 就变成 y 的 N-1 次多项式。道尔夫的做法是直接令这个多项式等于切比雪夫多项式 T_{N-1}(x0·y)其中 x0 是一个大于 1 的尺度因子。这样做的理由很实在T_m(x) 在 |x|≤1 区间里所有极值都等于 ±1等幅振荡在 |x|1 区间单调上升正好对应主瓣。把 |y|≤1 映射到旁瓣区、|y|1 映射到主瓣区就得到了“所有旁瓣等高”的效果这也是切比雪夫加权区别于汉明、汉宁这些窗函数的核心特征——后者旁瓣是递减的第一旁瓣最高越远越低。从多项式系数反推电流分布就是后面第 3 章脚本要做的事。这里先记住一点整条链路里唯一自由的参数就是 x0它由旁瓣电平决定跟阵元数、间距都无关。2.2 由旁瓣电平反解 x0 与主瓣展宽主瓣峰值出现在 ψ0此时 ycos(0)1T_{N-1}(x0) 就是主瓣电平第一旁瓣峰值出现在 ycos(ψ/2) 落到 ±1 的位置T_{N-1}(±1)±1。所以主瓣与旁瓣的电压比 R 满足T_{N-1}(x0) RR 10^(SLL/20)SLL 取正值例如目标 -30 dB 就代 30。切比雪夫多项式在 x1 时的显式表达是 T_m(x)cosh(m·arccosh(x))反解立刻得到x0 cosh( arccosh(R) / (N-1) )这个式子值得多看两眼。N 越大分母越大x0 越靠近 1意味着同样的旁瓣指标下多项式“挤”得更紧主瓣相对越窄——这就是大孔径带来的增益。反过来N 固定时想把 SLL 从 -20 dB 压到 -40 dBR 从 10 涨到 100arccosh 只按对数增长所以旁瓣指标的代价不是线性的压得越深越吃力。主瓣展宽没有闭式的简洁解工程上一般数值测量。定性规律是SLL 每降低 10 dB 左右-3 dB 波束宽度大约增加 10%~20%具体比例跟 N 有关。N 小的时候展宽更明显因为阵元数少、边缘电流被迫拉得很高。提示SLL 不要低于 -13.26 dB 去调用那个值附近切比雪夫解已经退化成均匀加权再低没有意义MATLAB 的 chebwin 会直接给警告。2.3 阵元间距和扫描角把边界卡死切比雪夫加权只管电流分布管不了栅瓣。栅瓣的出现条件是 ψ 在可见空间内出现多个 2π 周期即 d/λ ≥ 1/(1|sinθ_max|)。不扫描时 θ_max90°要求 d ≤ λ/2扫描到 60° 时要求 d ≤ λ/1.866 ≈ 0.536λ比半波长还宽松一点但如果扫描到 90° 端射d 必须小于 λ/2 才不出现栅瓣。实际工程里 d 通常取 0.5λ一是避开栅瓣二是让互耦可控。d 一旦超过 0.5λ方向图会在可见空间边缘冒出第二个大瓣切比雪夫加权压不住它——加权只改变各阵元的幅度不改变阵元间距决定的周期结构。扫描带来的第二个变化是波束展宽。阵列在 sinθ 空间的波束宽度是固定的换到 θ 空间要除以 cosθ0所以扫到 60° 时波束宽度大约是正扫时的 2 倍。这一点在做测角精度预算时必须算进去不能拿正扫的数字去套全扫描范围。3. 在 MATLAB 里生成加权系数chebwin 与自写递推两条路3.1 直接用 chebwin一行拿到系数但先核对三件事Signal Processing Toolbox 里的 chebwin 就是干这个的签名是 w chebwin(N, sll)N 是点数sll 是旁瓣相对主瓣的衰减单位 dB必须传正数。N 16; % 阵元数 sll 30; % 目标旁瓣 -30 dB这里传正数 d_lambda 0.5; % 阵元间距波长归一化 w chebwin(N, sll); % 返回 N×1 列向量峰值在中心 w w(:).; % 统一转行向量后面按行参与矩阵运算 w w / sum(w); % 归一化总激励为 1方便级联链路预算 fprintf(系数长度 %d最大值 %.4f最小值 %.4f\n, ... numel(w), max(w), min(w)); fprintf(动态范围比 %.2f (%.2f dB)\n, ... max(w)/min(w), 20*log10(max(w)/min(w)));这段代码的重点不在调用本身而在归一化。chebwin 返回的是未归一化的窗峰值在中点、取值约等于 1 的量级但总和不是 1。如果后面要跟发射功率、通道增益一起算链路预算统一归一化到 sum(w)1 更省事如果只是看方向图归不归一化都行因为画图前还会再除一次最大值。三件事要在调用前核对第一sll 必须是正数传 -30 会让内部算出 R1acosh 直接报错或返回复数第二N 为偶数时没有中心阵元chebwin 给的是对称的两半做相位中心对齐时要注意偏移 (N-1)/2 是半整数第三chebwin 是窗函数视角它的“主瓣”定义和阵列方向图的主瓣不完全是一回事两者数值上一致但概念要对齐别拿窗函数的等效噪声带宽去套阵列增益。3.2 自写递推脚本从 T_m 的升幂系数到阵元电流不想依赖工具箱或者想搞清楚系数到底怎么来的就自己写。思路分四步先用三项递推 T_{k1}2xT_k-T_{k-1} 求出 T_{N-1}(x) 的升幂系数再代入 xx0·y得到 y 的幂次系数然后用 cos^k(a)2^{-k}Σ C(k,j)cos((k-2j)a) 把 y^k 展成 cos(nψ/2) 的组合最后按阵元到中心的距离 r 取对应阶数的系数中心阵元取全值对称阵元各取一半。function w cheb_weight(N, sll_dB) % CHEB_WEIGHT N 元等距线阵的道尔夫-切比雪夫加权系数 % N 阵元数 2 % sll_dB 旁瓣电平正数例如 30 表示 -30 dB % 返回 1×N 行向量按阵元 1..N 排列和为 1 m N - 1; % 多项式阶数 阵元数 - 1 R 10^(sll_dB/20); % 主瓣/旁瓣电压比 x0 cosh(acosh(R)/m); % 峰值映射点 % ---- 1) 三项递推求 T_m(x) 的升幂系数 ---- Tkm1 1; % T_0 1 Tk [0; 1]; % T_1 x升幂排列 for k 1:(m-1) Tkp1 [0; 2*Tk] - [Tkm1; 0; 0]; % T_{k1} 2x*T_k - T_{k-1} Tkm1 Tk; Tk Tkp1; end cT Tk; % 长度 m1cT(k1) 是 x^k 的系数 % ---- 2) 代 x x0*cos(psi/2)得到 cos(psi/2) 的幂次系数 ---- kk (0:m).; cY cT .* (x0 .^ kk); % P(y) sum cY(k1)*y^k, y cos(psi/2) % ---- 3) 把 cos^k(psi/2) 展成 cos(n*psi/2) 的线性组合 ---- cCos zeros(1, m1); % cCos(n1) 对应 cos(n*psi/2) for k 0:m coef cY(k1) / 2^k; for j 0:k n abs(k - 2*j); % 倍角公式里的阶数 cCos(n1) cCos(n1) coef * nchoosek(k, j); end end % ---- 4) 按阵元到中心的距离取系数还原电流分布 ---- rlist abs((0:N-1) - (N-1)/2); % 每个阵元到中心的距离可为半整数 I zeros(1, N); for n 1:N A cCos(round(2*rlist(n))1); % cos(2r*psi/2) cos(r*psi) if rlist(n) 0 I(n) A; % 中心阵元独占 else I(n) A/2; % 左右对称阵元各分一半 end end w I / sum(I); % 归一化到总激励为 1 end逐步说明。第一步里 [Tkm1; 0; 0] 是在低位补两个零让长度和 [0; 2*Tk] 对齐这是升幂表示下做多项式加减最容易翻车的地方——补错位数不会报错只会算出一个看似合理但完全错误的多项式所以建议拿 N4、SLL30 手算一遍对照。第二步的 x0.^kk 在 N 较大时会比较大但还在 double 的范围内。第三步是整段代码的核心。cos^k(a) 的倍角展开用的是二项式系数nabs(k-2j) 会重复出现所以用累加而不是赋值这是“同阶合并”的关键漏掉累加的后果是方向图完全不成形。2^k 的除法必须放在累加之前做先累加后除会放大舍入误差。第四步的 rlist 对奇数 N 是整数序列对偶数 N 是半整数序列用 round(2*r)1 统一映射到 cCos 的下标这样奇偶两种情况用同一套索引逻辑不用分开写分支。中心阵元取全值、对称阵元取半值是因为展开式里 cos(rψ) 的系数对应的是“一对阵元合起来的贡献”。注意这个递推实现的上限大约在 N24 到 N28。超过之后 T_m 的系数和 x0^k 会一起变得很大正负抵消时 double 的 53 位有效位不够用算出来的系数会出现明显偏差甚至变号。N 超过这个量级就老老实实用 chebwin。3.3 归一化、动态范围与量化位宽系数算出来之后别急着画图先看两个数总和与动态范围比。总和归一化到 1 是为了让 AF(0)1方向图 dB 曲线的 0 dB 点就是主瓣峰值后面测量旁瓣时不用再折腾参考电平。动态范围比 DRR max(w)/min(w) 才是决定工程可实现性的指标。SLL-30 dB、N16 时DRR 大约在 12~15 dBSLL 压到 -40 dBDRR 会跳到 20 dB 以上-50 dB 时可能到 30 dB。这个数字直接对应功放的线性区要求中心阵元激励是边缘阵元的几百倍如果收发通道只有 30 dB 的线性动态范围边缘通道的信噪比会被压得很惨。做数字波束形成时还要多算一步量化。如果通道是 8 位 ADC满量程对应最大激励边缘阵元的量化台阶就是最大激励除以 256再乘上 1/DRR。用下面的小脚本快速估一下bits 8; % ADC 位宽 drr max(w)/min(w); % 动态范围比 lsb 1/(2^bits - 1); % 满量程归一化后的一个 LSB edge_level min(w)*drr*lsb; % 边缘通道的量化台阶占比 fprintf(DRR %.1f dB边缘通道量化台阶约为满量程的 %.4f%%\n, ... 20*log10(drr), 100*edge_level);如果这个百分比小到跟噪声底差不多说明旁瓣指标压得太狠、通道精度跟不上实测算出来的旁瓣根本到不了仿真值。这时候要么降指标要么换泰勒加权等旁瓣但边缘电流上升更缓要么给边缘通道单独加衰减器。4. 方向图仿真与指标测量从加权向量到 -3 dB 波束宽度4.1 阵因子计算与角度采样密度方向图计算本身只有三行但采样密度是个容易被忽略的坑。主瓣很窄的时候如果角度网格只取 0.1°波束宽度误差可能到 5% 以上旁瓣峰值会漏检——因为峰值点恰好落在两个采样点之间。N 16; d_lambda 0.5; sll 30; w chebwin(N, sll).; w w / sum(w); theta linspace(-90, 90, 18001); % 0.01° 网格兼顾精度和速度 psi 2*pi*d_lambda*sind(theta); % 空间相位 nvec (0:N-1).; AF w * exp(1j*nvec*psi); % 1×M 复阵因子 AF abs(AF); AF AF / max(AF); % 归一化到主瓣峰值 1 AF_dB 20*log10(max(AF, 1e-12)); % 加下限避免 log10(0) 出 -Infexp(1jnvecpsi) 这一步生成的是 N×M 的矩阵nvec 是列向量、psi 是行向量MATLAB 的隐式扩展会自动做外积。这里用的参考相位在阵元 0 上没有挪到阵列中心所以 AF 是个复数取 abs 之后相位因子 e^{j(N-1)ψ/2} 被消掉结果不受影响。但如果后面要做和差波束或者单脉冲测角相位必须保留那时候就得把参考点移到中心否则和差通道的相位关系会整体偏一个线性项。18 万点的网格跑起来大概几十毫秒N 再大也扛得住。如果做优化循环要反复调用可以降到 9001 点0.02°精度损失在 0.1% 量级速度翻倍。4.2 四种加权同图对比均匀、二项式、切比雪夫、泰勒单看一条切比雪夫曲线没有参照系把四种典型加权画在一张图上差异一眼就能看出来。% 四种加权系数 w_uni ones(1, N) / N; % 均匀加权 w_bin nchoosek(N-1, 0:N-1); w_bin w_bin/sum(w_bin); % 二项式 w_cheb chebwin(N, 30).; w_cheb w_cheb/sum(w_cheb); % 切比雪夫 -30 dB % 泰勒加权nbar4 控制近区旁瓣个数 if exist(taylorwin, file) w_tay taylorwin(N, 4, -30).; w_tay w_tay/sum(w_tay); else w_tay w_cheb; % 没装工具箱时退化处理仅保证脚本能跑 end figure; hold on; grid on; plot(theta, 20*log10(abs(w_uni * exp(1j*nvec*psi)) 1e-12), LineWidth, 1.1); plot(theta, 20*log10(abs(w_bin * exp(1j*nvec*psi)) 1e-12), LineWidth, 1.1); plot(theta, 20*log10(abs(w_cheb * exp(1j*nvec*psi)) 1e-12), LineWidth, 1.4); plot(theta, 20*log10(abs(w_tay * exp(1j*nvec*psi)) 1e-12), LineWidth, 1.1); xlabel(角度 \theta (deg)); ylabel(归一化幅度 (dB)); ylim([-80 0]); legend(均匀,二项式,切比雪夫 -30dB,泰勒 -30dB,Location,south);二项式系数用 nchoosek 一次生成它是所有 N 元加权里主瓣最宽、旁瓣为零的那一个常被拿来当“上界”参照。泰勒加权和切比雪夫的区别在于切比雪夫所有旁瓣严格等高泰勒让靠近主瓣的几个旁瓣保持等幅、远处旁瓣逐渐衰减代价是主瓣略宽一点点但边缘电流上升更缓DRR 通常小 2~4 dB。工程上如果指标只要求“第一旁瓣低于 -30 dB”而对远区旁瓣没硬性要求泰勒往往更划算。四种加权的典型指标对照如下N16、d0.5λ、正扫加权方式第一旁瓣 (dB)-3 dB 波束宽度系数 DRR (dB)特点均匀-13.26.4°0增益最高旁瓣最差二项式无单调下降10.8°约 24无旁瓣但主瓣最宽切比雪夫 -30 dB-30.09.1°约 13.5旁瓣严格等高泰勒 -30 dB-30.08.6°约 11远区旁瓣递减表里的数字跟 N 和 d 强相关换参数要重新测不能直接套。4.3 指标自动测量旁瓣、波束宽度与零点肉眼看图不够做参数扫描或者批量对比时必须自动测。三个指标各写一小段。第一旁瓣用局部极大值找排掉主瓣位置后取最大本文还有配套的精品资源点击获取
锦
锦皓数字建站
深耕本土企业品牌数字化升级,专注原创端正雅致商务官网,从视觉设计到稳定运维全程保驾护航。