简介资源包为LFM线性调频信号仿真与短时傅里叶变换STFT时频分析专题代码面向雷达信号处理、通信系统检测的初学者及需要理解非平稳信号分析的研究人员。内含1个MATLAB脚本.m约1KB配套实现LFM信号生成、加窗STFT计算、D倍抽取下采样提升时频分辨率以及时频图可视化等完整流程便于直接运行和二次修改。LFM频率线性变化常规STFT受窗长限制易出现时频模糊代码通过D倍抽取缓解该问题同时提示需防范混叠并合理选择下采样因子具有较强实操参考价值。已有880人学习可结合理论公式对照代码逐步掌握调频斜率、窗函数与分辨率的权衡。1. LFM 信号与 STFT 时频分析这份仿真代码到底能帮你什么做雷达信号处理或者通信系统仿真的人大概率都跟线性调频LFM信号打过交道。LFM 信号频率随时间线性变化宽频带、大时宽带宽积、抗干扰能力强几乎成了脉冲压缩体制的标准配置。但真正上手做仿真时有一个问题很磨人信号生成很简单cos一把就能造出来可当你需要观察它的时频特性尤其是瞬时频率随时间的变化趋势时直接做 FFT 得到的频谱完全看不出频率随时间的变化关系——因为频谱是全局的时间信息被积分掉了。这就是 STFT短时傅里叶变换要解决的问题给信号加一个滑动的窗得到频率成分随时间变化的二维时频分布。这份资源的核心就两样东西LFM.rar压缩包和里面的LFM.m脚本它完整实现了 LFM 信号生成、STFT 时频分析、D 倍抽取改善时频分辨率、可视化出图这一整条链路。适合雷达专业的学生做课程仿真也适合刚接触时频分析的一线工程师拿来当作 baseline直接改参数跑自己的数据。2. LFM 信号生成数学模型、参数选型与第一版代码2.1 从理想模型到可仿真信号载频、扫频率和初始相位的关系LFM 信号的数学表达式在多数教科书里长这样[ s(t) A \cos(2\pi(f_c t \frac{\beta}{2}t^2)) ]其中 ( A ) 是幅度( f_c ) 是初始载频t0 时刻的频率( \beta ) 是扫频率单位 Hz/s。注意相位是 ( 2\pi(f_c t \frac{\beta}{2}t^2) )对时间求导得到瞬时频率[ f_{inst}(t) f_c \beta t ]这地方有个容易犯迷糊的点瞬时频率是相位的导数不是相位直接除以 ( 2\pi t )。很多初学者拿cos(2*pi*fc*t pi*beta*t.^2)去仿真跟拿cos(2*pi*(fc*t beta*t.^2))去仿真得到的信号频率变化率差了一倍。写代码时相位表达式里t.^2的系数应该是pi*beta不是2*pi*beta原因就是求导后要乘出那个 2 来。这块搞错了后面 STFT 图上的频率范围会对不上理论值。再就是带宽。LFM 信号在持续时间 T 内扫过的频带宽度为 ( B |\beta| T )。工程上更常用的设计顺序是先定带宽 B 和脉冲宽度 T再反推扫频率 ( \beta B / T )。比如要设计一个带宽 20 MHz、脉宽 10 us 的雷达脉冲扫频率就是 ( 2 \times 10^{12} ) Hz/s。这种设计方式的好处是参数语义清晰带宽和脉宽直接跟系统指标挂钩。2.2 生成 LFM 信号的 MATLAB 代码与参数含义LFM.m里的信号生成部分常见的做法是这样写的%% 参数定义 fs 100e6; % 采样率 100 MHz T 10e-6; % 脉冲宽度 10 us B 20e6; % 信号带宽 20 MHz fc 10e6; % 初始载频 10 MHz A 1; % 信号幅度 beta B / T; % 扫频率 带宽 / 脉宽 N round(T * fs); % 总采样点数 t (0 : N-1) / fs; % 时间序列长度为 N % LFM 信号相位 2*pi*(fc*t 0.5*beta*t.^2) s A * cos(2*pi*fc*t pi*beta*t.^2); % 画个时域波形看眼 figure; plot(t*1e6, s); xlabel(时间 (us)); ylabel(幅度); title(LFM 时域波形);这里N round(T * fs)是总采样点数t (0:N-1)/fs生成 0 到 T 之间的 N 个时间点。pi*beta*t.^2这个写法对应前面讲的相位系数关系——如果在这里写成2*pi*beta*t.^2瞬时频率就会变成 ( f_c 2\beta t )扫频速率翻倍时频图上的斜率直接偏掉。fc选 10 MHz 是故意的让信号频率范围落在 10 MHz 到 30 MHz 之间全部低于奈奎斯特频率 50 MHz这样后面 STFT 出的时频图频率轴不会出现混叠。实际工程里fs的选择要满足带通采样定理如果信号带宽是 B中心频率很高可以用带通采样降低采样率。但如果只是仿真 LFM 的分析方法本身直接用低通等效模型更省事——把信号搬到零中频附近再做时频分析频率轴含义更清晰。2.3 信号参数与仿真参数对照参数符号本代码取值设计依据采样率fs100 MHz高于信号最高频率 30 MHz 的 2 倍以上脉冲宽度T10 us对应 1000 个采样点STFT 时窗长可调范围大信号带宽B20 MHz扫频范围 10~30 MHz带宽适中扫频率beta2e12 Hz/s等于 B/T采样点数N1000时频图时间轴分辨率的基础这套参数组合起来时频图上的瞬时频率线会从 10 MHz 斜着爬到 30 MHz斜率是 2 MHz/us。看 STFT 图时这条线的斜率就是验证代码正确性的第一道关卡——如果斜率对不上或者方向反了优先检查相位表达式里t.^2的系数。3. STFT 实现与窗函数选型短时傅里叶变换的运行细节3.1 为什么 LFM 信号需要 STFT 而不是直接做 FFTLFM 信号是非平稳信号它的频率随时间线性变化。直接做 FFT 得到的是整个时间窗口内的平均频谱结果是 20 MHz 带宽内所有频率分量的叠加完全看不出频率随时间的变化轨迹。STFT 的思路是把长信号切成一段一段的短信号每一段近似看成平稳信号做 FFT然后把各段的结果按时间顺序拼起来形成一张二维时频图。STFT 的数学定义是[ G(f, t) \int_{-\infty}^{\infty} x(\tau) w(t - \tau) e^{-j2\pi f \tau} d\tau ]这里的 ( w(t-\tau) ) 是窗函数它决定了在时间 t 附近哪些采样点参与本次 FFT。窗函数沿着时间轴滑动每次覆盖一小段信号得到一个频谱切片。把所有的切片堆叠成矩阵用imagesc或者surf画出来就是常说的语谱图或时频图。这个过程在 MATLAB 里最常用的函数是spectrogram底层就是这个公式的离散实现。LFM 信号的瞬时频率变化快如果窗选太长窗内频率变化范围大频谱会被展宽时频图上的线条会变糊如果窗选太短频率分辨率差同样看不清频率细节。时频分辨率冲突海森堡不确定原理在 LFM 这种快变信号上体现得格外明显这也是摘要里提到 D 倍抽取技术的动机——后面第四章细说。3.2 基于 spectrogram 的 STFT 代码与窗长选择LFM.m中 STFT 部分常见写法如下%% STFT 参数 win_len 128; % 窗长 128 点对应 1.28 us hop 32; % 帧移 32 点时间轴平滑度 nfft 512; % FFT 点数不够补零 % 汉明窗 win hamming(win_len); % STFT 时频矩阵 [S, f, t_stft] spectrogram(s, win, win_len - hop, nfft, fs); % 绘制灰度时频图 figure; imagesc(t_stft*1e6, f/1e6, 20*log10(abs(S) eps)); axis xy; xlabel(时间 (us)); ylabel(频率 (MHz)); title(LFM 信号 STFT 时频图); colorbar;spectrogram的关键参数有三个窗长win_len、帧移hop、FFT 点数nfft。窗长决定频率分辨率和时间分辨率的平衡点——128 点窗在 100 MHz 采样率下对应 1.28 us频率分辨率约等于 ( fs/win_len 0.78 ) MHz。帧移 32 点意味着相邻两帧有 96 点重叠时间轴不会出现明显的块状感。nfft设为 512 大于窗长是补零做插值让频率轴的网格更细但补零不会提升真实频率分辨率只是让频谱曲线看起来更平滑。加上eps是为了防止log10(0)出现负无穷警告这是绘图时的小习惯。绘制时用axis xy把 y 轴翻转成常规坐标方向否则imagesc默认 y 轴从左上角开始往下增长时频图看起来是倒的。3.3 窗函数对比汉明窗、海明窗与矩形窗的取舍窗函数主瓣宽度旁瓣衰减适用场景矩形窗最窄-13 dB频率分辨率优先且信号本身无泄漏问题时汉明窗较宽-43 dB通用场景谱泄漏小本项目默认选择海明窗较宽-53 dB需要更低旁瓣且能接受主瓣略宽时LFM 信号的频谱是连续扫过的用矩形窗会造成明显的频谱泄漏——时频图上除了主斜线之外还会有平行的亮度条纹那是旁瓣泄漏出来的虚假能量。汉明窗旁瓣衰减足够压制大部分泄漏且主瓣宽度增加不多在 LFM 时频分析里是性价比最高的选择。海明窗旁瓣更低但主瓣更宽如果两个频率分量靠得近可能分辨不开。工程上我会先用汉明窗跑一版看到谱线干净了就不动了如果有强干扰信号需要抑制再换海明窗试。nfft值的设置也值得说一句。如果原始信号长度 N 只有 1000 点win_len取 128nfft取 512时频矩阵的维度是 257 行、约 28 列取决于帧移画出来时间分辨率会比较粗。如果希望时频图时间轴更细腻可以把hop调小到 16 甚至 8代价是计算量增加和相邻帧相关性变强视觉上更平滑但并非信息量变多。4. D 倍抽取改善时频分辨率原理、实现与参数权衡4.1 为什么 STFT 之后时频分辨率还不够——问题的根源前面提到STFT 的时频分辨率受窗长限制窗长越短时间分辨率越高窗长越长频率分辨率越高。但 LFM 信号的瞬时频率随时间线性变化在任何一个窗内信号都不是纯单频而是覆盖一段频率范围 ( \beta \cdot T_{win} )。这个频率范围就是窗内频率展宽的来源。设窗长为 ( T_{win} )则窗内 LFM 信号的频带宽度为[ \Delta f_{win} \beta \cdot T_{win} ]而 STFT 能分辨的最小频率间隔大约是 ( 1 / T_{win} )。两相对比( \Delta f_{win} ) 远大于 ( 1/T_{win} )时频图上那条线就变粗了。比如前面参数下128 点窗时长 1.28 us窗内频率变化范围约 2.56 MHz而频率分辨率只有 0.78 MHz频率被展宽了 3 倍多。D 倍抽取的思路是既然窗内频率展宽是由信号本身扫频过快导致的那我把信号的采样率降下来让同样的时间跨度内采样点变少等效于用更长时间窗覆盖同样的扫频范围吗不是这个逻辑。D 倍抽取是对信号本身做整数倍抽取采样率从 fs 降到 fs/D时间跨度拉长 D 倍再做 STFT 时同一窗长对应的实际时间变长频率分辨率变为 ( fs/(D \cdot win_len) )从而在频率维度上看得更清楚。4.2 抽取的完整实现先抽取再用新采样率做 STFT%% D 倍抽取处理 D 4; % 抽取因子根据信号最高频率与抗混叠需求选 % 先做抗混叠滤波低通截止频率 fs/(2*D) % 这里用最简单的 FIR 低通实际工程会换成等纹波/窗函数设计的滤波器 cutoff fs / (2*D); % 抽取后奈奎斯特频率 norm_cutoff cutoff / (fs/2); % 归一化截止频率 fir_len 64; % 滤波器阶数 b fir1(fir_len, norm_cutoff); % 设计低通滤波器 s_filt filter(b, 1, s); % 滤波 % 抽取每 D 个点取一个 s_dec s_filt(1 : D : end); % 新采样率 fs_dec fs / D; % 用新的采样率重新做 STFT win_len_dec 128; % 同样的窗长对应实际时间变为原来的 D 倍 hop_dec 32; nfft_dec 512; [S_dec, f_dec, t_dec] spectrogram(s_dec, hamming(win_len_dec), ... win_len_dec - hop_dec, nfft_dec, fs_dec); figure; imagesc(t_dec*1e6, f_dec/1e6, 20*log10(abs(S_dec) eps)); axis xy; xlabel(时间 (us)); ylabel(频率 (MHz)); title([D num2str(D) 抽取后 STFT 时频图]); colorbar;这段代码里filter是抗混叠低通滤波fir1设计一个截止频率为 ( fs/(2D) ) 的 FIR 滤波器。s_dec s_filt(1:D:end)是标准的 D 倍抽取操作每隔 D 个点取一个。抽取后采样率降低为原来的 1/D时间轴刻度需要乘以 D 才能对上真实物理时间。fir1的参数norm_cutoff是归一化截止频率取值范围 0 到 11 对应奈奎斯特频率 ( fs/2 )。这里cutoff fs/(2*D)归一化后是1/D。对 D4 来说就是保留 0 到 12.5 MHz 的频率范围大于信号最高频率 30 MHz不对这里有个陷阱——如果直接对原始 10~30 MHz 的带通信号做低通滤波加抽取信号会被滤掉大部分因为 LFM 信号的频率范围在 10~30 MHz低通截止 12.5 MHz 会砍掉 12.5 MHz 以上的所有分量。4.3 D 倍抽取的混叠陷阱与正确前提抽取前必须先搞清楚信号频谱的位置。LFM 信号载频 10 MHz扫频到 30 MHz占据 10~30 MHz 的频带。D 倍抽取后采样率变为 fs/D奈奎斯特频率降为 fs/(2D)。想要抽取后信号不混叠必须满足两个条件之一一是原始信号最高频率低于抽取后奈奎斯特频率即 ( f_{max} fs/(2D) )二是信号是带通信号且抽取后频谱搬移能避开混叠——这时候需要对信号做带通滤波而非低通滤波而且需要用到带通采样的频率搬移思想。对上面的参数fs100 MHz, 信号 10~30 MHzD4 时抽取后采样率 25 MHz奈奎斯特频率 12.5 MHz信号最高频率 30 MHz 严重超出。直接抽取必然混叠时频图上会看到折叠回来的虚假频率分量。所以抽取的对象通常不是原始的带通 LFM 信号而是先做下变频乘以 ( e^{-j2\pi f_c t} )把信号搬到零中频附近变成低通信号后再抽取。零中频后信号占据 -10 MHz 到 10 MHz 的频带D4 时奈奎斯特频率 12.5 MHz 勉强够用。%% 正确流程先下变频再抽取 % 下变频乘以复指数把频谱搬移到零中频 t_full (0 : N-1) / fs; s_baseband s .* exp(-1j * 2 * pi * fc * t_full); % 这时候 s_baseband 是复数信号占据 -10 MHz ~ 10 MHz % 再做抗混叠低通 D 倍抽取 s_baseband_filt filter(b, 1, s_baseband); s_baseband_dec s_baseband_filt(1 : D : end); fs_dec fs / D; % 对零中频信号直接做 STFT [S_bb, f_bb, t_bb] spectrogram(s_baseband_dec, hamming(win_len_dec), ... win_len_dec - hop_dec, nfft_dec, fs_dec); figure; imagesc(t_bb*1e6, f_bb/1e6, 20*log10(abs(S_bb) eps)); axis xy; xlabel(时间 (us)); ylabel(频率 (MHz)); title(零中频 LFM 抽取后 STFT 时频图); colorbar;下变频的本质是频谱搬移乘上 ( e^{-j2\pi f_c t} ) 后原来 10~30 MHz 的频谱被搬到了 -10 MHz~10 MHz。这时再做低通滤波抽取滤波器截止频率设为 fs/(2D) 就合理了。抽取后的时频图中心在 0 MHz扫频线从 -10 MHz 斜着上升到 10 MHz物理意义等价于原信号的瞬时频率减去载频。如果需要把频率轴还原成真实频率显示时加回 fc 即可f_display f_bb fc。D 的取值直接决定时频图质量和计算开销。D 越大抽取后采样率越低同样窗长下的频率分辨率越高但代价是频率观测范围变窄和可能的混叠风险。一般工程上 D 取 2 到 8 之间具体看信号带宽和中间频率的比值。判断 D 是否合适的简单办法抽取后最高频率不超过新奈奎斯特频率的 80%留出过渡带余量。5. 避坑手册LFM 仿真与 STFT 时频分析的五个常见问题5.1 扫频方向看反了相位系数写错导致时频图斜率方向错误现象时频图上瞬时频率线在减小从高频到低频但理论值明明是频率随时间增加。原因相位表达式中t.^2的系数写错。正确写法是pi*beta*t.^2如果写成2*pi*beta*t.^2等效扫频率变成 2β瞬时频率公式变成 ( f_c 2\beta t )。有些情况下甚至写成了-pi*beta*t.^2扫频方向就反了。还有一种可能是 beta 本身取了负值下扫频设计但时频图判读时没注意符号。解决先在代码里加一行验算——直接用相位对时间求导的数值解f_inst_theory fc beta * t再把理论线叠加到时频图上看斜率是否吻合。代码层面对照公式检查系数cos(2*pi*fc*t pi*beta*t.^2)系数 pi 不对就改回来。5.2 抽取后频谱混叠出一堆假信号现象D 倍抽取后 STFT 时频图上出现多余的不规则亮线真实扫频线反而看不清。原因没有先下变频对带通信号直接低通滤波 抽取。原始信号频率 10~30 MHz 超过抽取后奈奎斯特频率高频分量折叠回低频区域产生虚假分量。解决抽取前先做频谱搬移乘法操作s .* exp(-1j*2*pi*fc*t)把信号中心频率移到 0。搬移后确认信号的最高频率低于fs/(2D)。如果不方便下变频就需要把低通滤波器换成带通滤波器让感兴趣的频带在抽取后落在新奈奎斯特频率以内。5.3 时频图上下颠倒imagesc 的 y 轴方向陷阱现象频率轴显示为从大到小低频在图片上方高频在下方跟教科书习惯相反。原因MATLAB 的imagesc默认 y 轴方向是反的第一个数组行对应图像顶部。STFT 矩阵的第一行是低频分量于是低频被画到了图片上方。解决绘制后紧跟axis xy;命令翻转 y 轴方向。或者用set(gca, YDir, normal)效果相同。这条几乎每次用imagesc画频谱都会遇到属于高频复发型踩坑。5.4 时频线条太粗窗长与扫频速率的匹配问题现象时频图上的扫频线不是一条细线而是模糊的宽带带看不出瞬时频率的精确位置。原因窗长内频率变化范围 ( \beta \cdot T_{win} ) 过大。窗长 128 点对应 1.28 us扫频率 2e12 Hz/s 时窗内频率变化 2.56 MHz而频率分辨率只有 0.78 MHz线条被展宽到几个像素宽度。解决优先缩减窗长而不是依赖 D 倍抽取。把win_len从 128 降到 64窗内频率变化范围减半到 1.28 MHz。代价是频率分辨率降到 1.56 MHz所以需要配合 D 倍抽取来补频率分辨率先抽取再缩短窗长才能同时保时间和频率的分辨率。这一招是从工程调试里学来的两者联动调而不是单独调一个参数。5.5 频率轴和理论值对不上采样率参数传递错误现象时频图上扫频线起点和终点频率跟设定的 10 MHz/30 MHz 误差很大或者整体偏移。原因spectrogram的采样率参数fs与信号生成时的fs不一致。常见于抽取后没有把fs更新为fs/D还沿用原来的 100 MHz 传参导致频率标尺放大 D 倍。另一个可能原因是下变频后忘记加回 fc频率轴显示的是相对频率。解决每次改动采样率后在代码里用assert(fs_dec fs/D)做一次断言。绘制频率轴时明确注释是绝对频率还是相对频率。显示时如果需要绝对频率用f_abs f_bb fc生成新的频率轴向量再传入imagesc。6. 时频分析结果的验证技巧从肉眼确认到定量校验时频图画出来只是第一步怎么确认它是对的才是真正考验工程判断力的地方。我的习惯是同时做三件事第一把理论瞬时频率线叠加到时频图上。第二个方法能交叉验证 STFT 和 D 倍抽取的联合正确性——用已知结果的信号跑整个链路如果输出符合预期链路就没问题。第三个方法最实用也最容易被人忽略抽取后的时频图频率轴需要跟理论带宽对得上。%% 验证理论瞬时频率叠加到时频图上 % 理论瞬时频率绝对频率 f_inst_theory fc beta * t; % t 是原始信号时间轴 % 抽取后的时频图频率轴是相对频率加回 fc 才能比较 f_abs_display f_bb fc; figure; imagesc(t_bb*1e6, f_abs_display/1e6, 20*log10(abs(S_bb) eps)); axis xy; hold on; plot(t*1e6, f_inst_theory/1e6, r--, LineWidth, 1.5); xlabel(时间 (us)); ylabel(频率 (MHz)); title(理论瞬时频率与时频图叠加验证); legend(STFT 时频图, 理论瞬时频率);理论线f_inst_theory是从 10 MHz 到 30 MHz 的直线斜率为 β。如果抽取链路正确这条红线应该刚好落在时频图亮带的中心位置。有小偏差是正常的STFT 的时间-频率分辨率有限但偏差不应该超过一个窗长对应的频率范围 ( \beta \cdot T_{win} )。如果偏差超过这个范围说明滤波器的群延迟没补偿——FIR 滤波有固定的群延迟抽取后的时间轴整体偏移了半个滤波器阶数的采样点数。群延迟补偿是最后一个容易被忽略的坑。filter函数输出的群延迟约为(fir_len - 1) / 2个采样点在抽取后的时间轴上对应(fir_len - 1) / (2 * D)个点。画图时把时间轴减去这个偏移量t_corrected t_bb - (fir_len - 1) / (2 * D * fs_dec);叠加理论线时才不会出现固定时间偏移。从那以后我每次分析时频图都强制走一遍「下变频 → 滤波 → 抽取 → STFT → 叠加理论线」的完整链路确认所有部件都正常了才输出结论。希望这个方法能帮你在自己的 LFM 仿真里少走几步弯路。本文还有配套的精品资源点击获取
