雷达模糊函数MATLAB仿真:矩形、高斯与LFM波形对比

发布时间:2026/10/26 11:40:13
雷达模糊函数MATLAB仿真:矩形、高斯与LFM波形对比

简介:这套 MATLAB 仿真代码聚焦雷达与通信信号分析中的模糊函数计算,面向学习信号处理、雷达波形设计的本科生、研究生及工程技术人员,帮助其理解矩形脉冲、高斯脉冲与线性调频(LFM)信号的时频耦合特性及分辨性能。压缩包共 6 个文件,均为 .m 脚本与函数,体积约 4KB,分别完成三类信号的模糊函数仿真、高斯脉冲与矩形脉冲的独立函数实现,以及将模糊函数计算封装为 ambiguity() 以便直接调用;配套脚本可绘制模糊函数图、模糊度图、距离模糊函数图和速度模糊函数图,便于对照比较不同波形的距离—速度分辨能力。已有 4217 人学习,说明其在小规模实验复现中具有一定参考价值。读者可据此快速搭建波形模糊函数分析框架,理解封装函数的输入输出与调用方式,并在此基础上替换参数或扩展其他调制信号,用于课程实验、专题报告或算法预研中的波形对比与结果验证。

1. 两个目标差 15 米、速度差 30 米每秒,模糊函数就是回答“分不分得开”

某实验室做近程探测波形选型时遇到一个具体问题:两个目标在距离上相差 15 m,径向速度相差 30 m/s,接收机用匹配滤波处理,能不能把这两个回波分成两个峰?距离差 15 m 折合成双程时延 0.1 μs,速度差 30 m/s 在 3 cm 波长下折合成多普勒 2 kHz。把这对 (τ, ξ) 代入候选波形的模糊函数,如果归一化幅度掉到 −3 dB 以下,说明这两个目标在匹配滤波输出上会互相顶开、形成两个可分辨的峰;如果还挂在 −1 dB 上,那就是一片糊,换什么检测算法都救不回来。

模糊函数 χ(τ, ξ) 就是把“时延差 + 多普勒差”这组二维坐标映射成匹配滤波输出复包络的那个函数。它一次性说清了三件事:分辨率有多细、旁瓣有多高、距离和速度之间会不会互相串扰。矩形脉冲、高斯脉冲、LFM 线性调频这三类基带波形,是几乎所有雷达波形设计教材的起点,也是最容易在 MATLAB 里用三十行代码跑通、跑对、跑出结构差异的三个样本。下面从离散化开始,一路写到能直接量出主瓣宽度和峰值旁瓣比的脚本。适合正在做波形选型、或者第一次上手模糊函数仿真的人照着复现。

2. 把 χ(τ,ξ) 拆成矩阵乘:离散化、网格选取与三种波形的生成

2.1 时延用整数移位、多普勒用复指数:离散化绕不开的两个硬约束

连续形式的模糊函数有几种等价约定,工程上最常用的是这一种:

χ(τ, ξ) = ∫ s(t) · s*(t − τ) · exp(j2πξt) dt

其中 s(t) 是复基带信号,τ 是时延,ξ 是多普勒频率。取 τ = 0、ξ = 0 时 χ(0,0) = ∫|s(t)|² dt,就是信号能量。后面所有图都按这个值归一化到 0 dB。

落到离散域,这个积分变成三重求和:对每个时延格点 k、每个多普勒格点 m,做一次长度为 N 的内积。直接三重循环写出来能跑,但在 MATLAB 里的速度会让人怀疑人生。正确的做法是把多普勒那一维提出来做成矩阵,让内积变成一次矩阵乘法。

这里有两个绕不开的硬约束,先讲清楚,后面代码才看得懂。

第一个约束是时延分辨率被采样率限制住。s(t − τ) 在离散域只能是 s(n − ns),ns 必须是整数。也就是说,τ 只能取 1/fs 的整数倍,没有中间值。想让时延轴更细,只有提高采样率或者做频域插值上采样,没有第三条路。很多人在这一步用linspace生成一个看起来很密的时延轴,然后在循环里对 τ·fs 做round,结果图上会出现周期性的阶梯条纹,这就是典型的翻车点。

第二个约束是多普勒轴的周期性。exp(j2πξt) 在 t = n/fs 上采样后,对 ξ 的周期是 fs。所以离散计算出来的模糊函数在多普勒方向上以 fs 为周期重复。物理上有意义的区间只有 [−fs/2, fs/2],超过这个范围算出来的东西是混叠出来的镜像,看着像真的,其实是鬼影。

把这两个约束合起来看,一幅可信的模糊函数图,前提条件是 fs 明显大于信号带宽 B(时延方向采样够密),同时多普勒轴范围不超过 fs/2。

2.2 时延轴取几个 1/B、多普勒轴取几个 1/T:轴范围定错的图没有信息量

轴范围比代码更容易出错。一个常见现象是:图画出来了,花花绿绿一大片,但完全看不出三种波形的差别。八成是轴范围取错了——要么太宽,主瓣缩成一条线;要么太窄,只看到中心一个点。

定轴范围的依据是两个特征尺度:时延方向的特征宽度是 1/B(B 为信号带宽),多普勒方向的特征宽度是 1/T(T 为信号有效时宽)。这两个量分别是距离分辨率和速度分辨率的量纲来源。

波形时延轴建议范围多普勒轴建议范围图上应该看到什么
矩形脉冲±1.1·T±3/T三角形时延脊,多普勒方向 sinc 旁瓣
高斯脉冲±3σ±3/(2πσ)一片光滑的二维高斯,无可见旁瓣
LFM(B·T = 100)±2/B±2/T一条斜率 K 的斜刃,沿刃线有 sinc 起伏

LFM 那一行的轴范围值得多说一句。斜刃的脊线方程是 ξ = Kτ,K = B/T。如果时延轴取到 2/B,那么脊线在这段时延上扫过的多普勒范围正好是 K·(2/B) = 2/T,跟多普勒轴范围对上了。这就是为什么 LFM 的时延轴要用 1/B 而不是 T 来定标——用 T 定标的话,你看到的只是斜刃最中间那一小段,几乎是一条水平线,看不出倾斜。

提示:三个波形的轴范围不要强行统一。模糊函数本身是各自归一化的,轴各自按特征尺度取,图才有可比性。硬要放在同一组坐标下,LFM 的斜刃会被压扁成一条线。

2.3 矩形、高斯、LFM 三种基带信号的生成与能量归一化

信号生成这一步有个隐藏的坑:能量归一化。如果不做归一化,矩形脉冲的能量是 T,LFM 也是 T,但高斯脉冲取决于 σ,三个波形的 χ(0,0) 完全不同,画出来的 dB 图没法横向比较。统一的做法是让每个信号的离散能量 Σ|s(n)|²·dt 等于 1,这样 χ(0,0) 恒为 1,任何一张图上的 −10 dB 都是同一个含义。

function [s, t, fs] = make_waveform(type, T, fs, param) % 生成复基带波形,能量归一化到 sum(|s|^2)*dt = 1 % type : 'rect' | 'gauss' | 'lfm' % T : 信号窗长(秒) % fs : 采样率(Hz) % param: gauss 时为 sigma(秒);lfm 时为带宽 B(Hz) N = 2*floor(T*fs/2) + 1; % 取奇数点,保证有严格的零时刻样点 t = ((0:N-1).' - (N-1)/2) / fs; % 时间轴以中心为零点 switch lower(type) case 'rect' s = ones(N,1); case 'gauss' sigma = param; s = exp(-t.^2 / (2*sigma^2)); case 'lfm' B = param; K = B / T; % 调频斜率 s = exp(1j*pi*K*t.^2); % 复基带 LFM otherwise error('unknown waveform type'); end s = s / sqrt(sum(abs(s).^2)/fs); % 能量归一化 end

三个分支里只有 LFM 是复数信号,矩形和高斯生成出来是实向量,计算时会自动按复数处理,不用额外转换。

N取奇数这一步不是洁癖。偶数点时时间轴没有严格落在 0 上的样点,最接近的两个是 ±dt/2,归一化后 χ(0,0) 会略微偏小,而且零多普勒切面会左右不对称,肉眼能看出来。

t的零点放在中心而不是起点,是为了让 exp(j2πξt) 的相位在信号两端对称。严格来说多普勒相位的绝对时间原点只影响一个常数相位因子,不影响幅度,但对称取法能让后面画出来的相位图好看一些,也便于和解析式对齐。

参数上,本文一律取 T = 10 μs、fs = 100 MHz,此时 N = 1001,每个时延格点间隔 10 ns。LFM 取 B = 10 MHz,时宽带宽积 D = B·T = 100。高斯取 σ = T/8 = 1.25 μs,窗口 ±5 μs 对应 ±4σ,截断处的幅度是 exp(−8) ≈ −69 dB,远低于后面 −40 dB 的显示门限,不会污染旁瓣。

2.4 一个 30 行的 ambg 函数:一次矩阵乘算出整张模糊面

核心函数长这样。思路是把多普勒相位做成一个 N×Nfd 的矩阵,对每个时延做一次sd' * SE,一次矩阵乘把整行多普勒算完。

function [AF, tau, fd] = ambg(s, fs, tauMax, fdMax, Nfd) % AMBG 离散模糊函数 chi(tau,fd) = sum s(n)*conj(s(n-ns))*exp(j2*pi*fd*t)*dt % s : 复基带样本(列向量),已能量归一化 % fs : 采样率 % tauMax : 时延轴单边最大值(秒),实际会被截到 1/fs 的整数倍 % fdMax : 多普勒轴单边最大值(Hz) % Nfd : 多普勒格点数,内部强制为奇数以包含零频 s = s(:); N = numel(s); dt = 1/fs; t = ((0:N-1).' - (N-1)/2) / fs; % 与生成时一致的中心化时间轴 if mod(Nfd,2) == 0, Nfd = Nfd + 1; end fd = linspace(-fdMax, fdMax, Nfd).'; % 多普勒网格 K = floor(tauMax*fs); % 时延只能取整数个样点 ns = (-K:K).'; % 移位样点数,含 0 tau = ns / fs; % 对应的时延轴(秒) E = exp(1j*2*pi*(t*fd.')); % N x Nfd 多普勒相位矩阵 SE = s .* E; % 先把信号调制到各多普勒频点 AF = zeros(numel(ns), Nfd); for k = 1:numel(ns) n0 = ns(k); sd = zeros(N,1); % 补零移位,绝不用 circshift if n0 >= 0 sd(n0+1:end) = s(1:N-n0); % s(t - tau),延迟 else sd(1:N+n0) = s(1-n0:end); % 超前 end AF(k,:) = (sd' * SE) * dt; % 1xN 乘以 NxNfd end end

几个参数和实现细节值得逐条交代。

tauMax传进来之后被floor(tauMax*fs)截断,实际时延轴范围是 ±K/fs。这意味着如果你传 10 μs、fs = 100 MHz,K = 1000,正好对上;如果传 10.005 μs,会被截成 10 μs。这个“静默截断”是有意为之——与其让你以为轴范围是按传入值算的,不如强制对齐到物理可实现的分辨率上。

sd用补零移位而不是circshift。circshift会把移出窗口的样本从另一端绕回来,等效于把信号当成周期信号处理,结果是在 τ = ±T 附近凭空长出旁瓣,而且 τ 越大旁瓣越强,看起来像是某种物理效应,实际是纯粹的数值伪影。这个坑后面还会专门讲。

(sd' * SE) * dt里的dt不能省。省掉之后 χ(0,0) 等于 N 而不是 1,不同 fs 下算出来的数值不可比,而且和解析式对不上。乘上 dt 之后,离散求和就是积分的一致逼近,fs 加大时数值会稳定收敛。

内存和耗时方面,矩形脉冲那一组 N = 1001、Ntau = 2001、Nfd = 201,SE占 1001×201 个复数,约 3 MB;循环 2001 次,每次做 1001×201 的内积,总运算量约 4×10⁸ 次复乘,普通笔记本上一两秒能出结果。如果嫌慢,把 Ntau 降到 401(轴范围不变、格点变稀)也不影响主瓣形状,只是等高线会粗糙一点。

3. 三种波形的模糊图长什么样:矩形、高斯与 LFM 的结构差异

3.1 矩形脉冲:三角形时延脊与 ±1/T 处的 sinc 旁瓣

矩形脉冲的模糊函数有闭式解,写成归一化形式是:

|χ(τ, ξ)| = |(1 − |τ|/T) · sinc(ξ(T − |τ|))|, |τ| ≤ T

MATLAB 的sinc(x)定义为 sin(πx)/(πx),零点在非零整数处,所以多普勒方向的第一零点出现在 ξ = 1/(T − |τ|),在 τ = 0 处就是 1/T。时延方向则是一个严格的三角形,τ = ±T 处降到零。

这两个特征直接给出了矩形脉冲的两条设计结论。距离分辨率由时延三角形的半宽决定,量级是 T——10 μs 的脉冲,双程时延分辨率 10 μs,折合单程距离 1.5 km。这个数非常糟糕,正因为如此,单纯用矩形脉冲做雷达测距在现代系统里基本见不到。速度分辨率由 sinc 主瓣宽度决定,量级是 1/T = 100 kHz,反倒还不错。

旁瓣方面,sinc 函数的第一个旁瓣在 −13.2 dB,之后按 1/ξ 衰减。这一点和高斯脉冲形成鲜明对比。更要命的是,时延方向虽然主瓣之外严格为零(三角形在 |τ| > T 处截断),但多普勒方向的旁瓣在整个时延范围内都存在,只是幅度被三角形因子压小。表现在等高线图上,就是一条沿 τ 轴延伸的三角形亮脊,两侧挂着一条条平行的旁瓣条纹。

算这张图的时候要注意 K 的取值。τ 到了 ±T 边缘时,(1 − |τ|/T) 因子归零,sinc的参数也归零,0 乘上一个有限的 sinc 值仍然是 0。数值上不会出问题,但如果你用解析式做对照,记得在 |τ| ≥ T 时直接置零,别让sinc的振荡漏出来。

3.2 高斯脉冲:二维高斯、全程无旁瓣,代价藏在哪

高斯脉冲 s(t) = exp(−t²/(2σ²)) 的模糊函数是可以分离变量的:

|χ(τ, ξ)| / χ(0,0) = exp(−τ²/(4σ²)) · exp(−π²σ²ξ²)

时延方向是标准差为 √2·σ 的高斯,多普勒方向是标准差为 1/(πσ√2) 的高斯。两个方向各自单调衰减,永远不穿过零点,也就永远没有旁瓣。在 −40 dB 的动态范围图上,你看到的是一团干净的椭圆光斑,四周什么都没有。

这个性质的稀缺程度值得强调一下。模糊函数有一个“体积不变”约束:∫∫|χ(τ,ξ)|² dτ dξ = |χ(0,0)|²。归一化之后,整张曲面下的体积恒等于 1。这意味着你没法让所有地方都小——压低了旁瓣,主瓣必然变宽,或者能量跑到别的什么地方去。高斯脉冲是唯一能在不产生任何旁瓣的前提下把这个体积塞进主瓣里的波形。

但代价同样明确。高斯脉冲没有锐截止,时延方向的高斯尾巴拖得很长,两个相距 3σ 的目标,第二个的峰值会叠在第一个的尾巴上,虽然不会被误判成旁瓣,但会造成峰值位置偏移。另外,真正的雷达发射机做不出无限长的高斯包络,总要截断,一截断就重新长出旁瓣。本文的 σ = T/8,窗口 ±4σ,截断处 −69 dB,如果换成 σ = T/4、窗口 ±2σ,截断处只有 −8.7 dB,画出来的图会完全失去“无旁瓣”这个卖点。这一点在第 5 章还会再提一次。

还有一点容易被忽略:高斯脉冲的时宽带宽积有一个理论下界。它的时间标准差和频率标准差之积恒为 1/(4π),无论 σ 怎么调都变不了。这意味着高斯脉冲的分辨率乘积是一个固定值,想同时提高距离和速度分辨率,高斯脉冲给不了你自由度。

3.3 LFM:斜率为 K 的刀刃与距离-多普勒耦合

LFM 的模糊函数近似式是:

|χ(τ, ξ)| ≈ |(1 − |τ|/T) · sinc((ξ − Kτ)(T − |τ|))|, |τ| ≤ T

跟矩形脉冲的式子只差一个− Kτ,但这一项改变了整张图的拓扑结构。矩形脉冲的脊线在 ξ = 0 上,是一条水平的亮带;LFM 的脊线在 ξ = Kτ 上,是一条斜线,斜率就是调频斜率 K。

这条斜线的物理含义是距离-多普勒耦合:一个速度不为零的目标,回波的多普勒频移会让匹配滤波输出的峰值在时延轴上偏移 Δτ = ξ/K。用 B = 10 MHz、T = 10 μs 的参数,K = 10¹² Hz/s,多普勒 100 kHz 对应的时延偏移是 0.1 μs,也就是一个 1/B。这正是 LFM 在多目标场景下需要额外做速度补偿的原因——不补偿的话,测出来的距离带着一个正比于径向速度的偏差。

从图上看,这条斜刃在归一化坐标下是 45° 的。取时延轴 ±2/B、多普勒轴 ±2/T,脊线 ξ = Kτ 在这段时延上扫过多普勒 ±2/T,正好覆盖整个显示范围。沿着刃线走,幅度被 sinc 调制,出现一串亮点和暗点,间隔是 1/(T − |τ|)——这就是 LFM 在“距离-速度联合平面”上的旁瓣结构。垂直于刃线的方向衰减极快,因为 sinc 的自变量里 (T − |τ|) 很大,等于把带宽压缩比 D = 100 全部用上了。

这份“刃很薄、但刃本身是斜的”的结构,是 LFM 最核心的取舍:它用 D = B·T 的压缩比换来了 1/B 量级的距离分辨率和 1/T 量级的速度分辨率同时成立,矩形脉冲做不到这一点。代价就是那条斜刃——耦合带来的距离偏差必须在信号处理链路里补掉。

4. 出图、切面与读数:把 dB 等高线变成分辨率、PSLR 和主瓣宽度

4.1 归一化到 0 dB 再截 −40 dB:动态范围选错会把旁瓣藏起来

abs(AF)的绝对值没有任何可比性,先除以全局最大值归一化,再取 20log10。截断门限选多少,直接决定你能看到什么。

AFn = abs(AF) / max(abs(AF(:))); % 归一化到 0 dB L = 20*log10(AFn + 1e-12); % 加 eps 防止 log10(0) L(L < -40) = -40; % 低于 -40 dB 的全部压平 figure('Position',[100 100 720 520]); contourf(tau*1e6, fd*1e-3, L.', ... [-40 -35 -30 -25 -20 -15 -10 -6 -3 -1 0], 'LineStyle','none'); colormap(flipud(gray)); colorbar; caxis([-40 0]); xlabel('\tau (\mus)'); ylabel('\xi (kHz)'); title(sprintf('%s |\\chi| 归一化 dB', name));

几点说明。contourf的 X、Y 长度要分别等于数据的列数和行数,而AF是 Ntau×Nfd,所以传入的是L.',X 传tau(长度 Ntau),Y 传fd(长度 Nfd)。这个转置漏掉的话,要么报错,要么画出来的图轴标签是反的——而且不报错的概率更大,因为 Ntau 和 Nfd 有时候恰好相等。

门限取 −40 dB 是个工程折中。低于 −40 dB 的旁瓣对检测性能的影响已经可以忽略(噪声基底通常在这个量级以上),压平之后图的层次更清楚。如果做的是低截获概率波形的旁瓣分析,需要看到 −60 dB,那就改成L(L < -60) = -60,同时把caxis和contourf的等值线向量一起改掉,三处不一致的话会出现一块纯白或者纯黑。

colormap(flipud(gray))让高值变黑、低值变白,打印出来更清楚。想要彩色就用parula或jet,但jet在等高线图上容易产生视觉上的假边界,用在论文插图里会被审稿人挑。

4.2 零多普勒切面与零时延切面的两种读法

二维等高线看结构,一维切面读数。两个切面各有分工:零多普勒切面(ξ = 0)量距离分辨率,零时延切面(τ = 0)量速度分辨率。

[~, i0] = min(abs(tau)); % tau = 0 的行号 [~, j0] = min(abs(fd)); % fd = 0 的列号 cut_tau = 20*log10(abs(AF(:,j0)) + 1e-12); % 沿 tau,看距离分辨 cut_fd = 20*log10(abs(AF(i0,:)) + 1e-12); % 沿 fd,看速度分辨 % -3 dB 主瓣宽度(以零多普勒切面为例) c = <p> <a href="https://download.csdn.net/download/zrh86/85303218" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源,点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>