简介:该资源围绕线性调频(LFM)信号生成及其短时傅里叶变换(STFT)时频分析,提供一份可直接运行的MATLAB脚本,适合学习雷达信号处理、通信系统信号检测以及非平稳信号时频分析的高校学生或工程初学者。代码涵盖LFM信号数学模型构建、窗函数选择、STFT时频图绘制,并引入D倍抽取策略以改善瞬时频率快速变化导致的时频分辨率不足,同时提示下采样因子选取不当可能引发的混叠风险。压缩包内仅含1个m文件,整体大小约1KB,结构精简,便于读者逐行阅读、修改参数并对比不同窗函数与抽取倍数对时频结果的影响,也可作为教学演示或课程设计的基础模板。目前该资源已有881人学习下载,对于希望快速掌握LFM信号特性及STFT工程实现细节的读者,是一份轻量实用的入门参考。
1. 把LFM信号看明白:STFT时频分析与D倍抽取仿真代码怎么用
做雷达或通信的人大概都有这种经历:拿到一段LFM信号,时域波形看半天只是幅度不变的震荡,换成时频图才能看到那条从低频扫到高频的斜线。这份LFM仿真代码包(LFM.rar 里的 LFM.m)做的就是这件事——生成LFM信号、经过STFT处理画时频图,再用D倍抽取把时频分辨率往上提一档。它解决三个核心问题:LFM怎么从数学表达式变成可运行的MATLAB代码;STFT的窗长、hop怎么选才不出废图;时频图不清楚时,D倍抽取怎么落进代码、边界在哪。做雷达信号课设、通信仿真,或刚接触chirp信号分析的人,可以直接拿这个脚本对照改参数。下面按信号生成、STFT、D倍抽取、避坑的顺序拆开讲。
2. LFM信号生成:数学表达式落地成MATLAB代码的关键参数
2.1 LFM数学表达:扫频率β才是核心
LFM信号的标准表达式是 s(t) = A·cos(2π(f_c·t + βt²/2)),式子本身不难,但很多人在写代码前没搞明白一个点:β才是这个信号的灵魂。β是扫频率,单位是Hz/s,它决定了信号频率每秒变化多少。对相位项 2π(f_c·t + βt²/2) 求导再除以2π,就得到瞬时频率 f_inst(t) = f_c + βt。正是这个一次函数关系,让LFM在时频图上呈现一条斜线。
为什么相位里要放一个二次项?因为频率是相位的导数,要让频率随时间线性变化,相位就必须是时间的二次函数。这个关系想通了,后面写代码和验证时频图就都有了理论基准。工程上LFM大量用于雷达脉冲压缩:发射一个宽脉冲,接收时通过匹配滤波压缩成窄脉冲,而脉冲压缩比恰好等于时宽带宽积BT。所以仿真时B和T不是随便填的,它们共同决定了信号的性能指标。
知道了这些再回头看参数:fc是起始频率,B是扫频带宽,T是脉宽,β = B/T。例如fc=200Hz、B=400Hz、T=0.5s时,信号从200Hz线性扫到600Hz,扫频率β=800Hz/s。代码里所有与频率相关的量都要围绕这套关系展开,下面是这三个参数的典型关系:
| 参数 | 符号 | 本例取值 | 单位 | 说明 |
|---|---|---|---|---|
| 采样率 | fs | 2000 | Hz | 必须大于最高瞬时频率的两倍 |
| 起始频率 | fc | 200 | Hz | 扫频起点 |
| 带宽 | B | 400 | Hz | 扫频范围,决定脉压比 |
| 脉宽 | T | 0.5 | s | 信号持续时间 |
| 扫频率 | β | 800 | Hz/s | β = B/T,频率变化快慢 |
2.2 LFM.m中的信号生成代码与逐行拆解
LFM.m的习惯写法是先设参数、再建时间轴、最后用一行核心表达式生成信号。下面这段就是脚本的开头:
% LFM信号生成 —— LFM.m 开头 clear; close all; clc; fs = 2000; % 采样率 2000 Hz fc = 200; % 起始频率 200 Hz B = 400; % 带宽 400 Hz T = 0.5; % 脉宽 0.5 s A = 1; % 幅度 beta = B / T; % 扫频率 = 800 Hz/s t = 0 : 1/fs : T - 1/fs; % 时间轴,末尾减一个采样间隔 N = length(t); % 总样本点数 % 核心表达式:把瞬时相位写进cos里 s = A * cos(2 * pi * (fc * t + 0.5 * beta * t.^2));逻辑说明:时间轴用冒号生成,从0开始,步长1/fs,到T-1/fs结束,总长度正好N = fs*T,不会多出第N+1个点。t.^2用点乘方逐元素平方,不能用t^2,矩阵乘方会直接报错或算出完全错误的结果。核心表达式把fc*t和0.5*beta*t.^2相加再乘2π,正好对应相位项2π(f_c·t + βt²/2),cos外面乘A就是幅度。
参数说明:fs=2000Hz是后面STFT能正常工作的前提。这个信号从200Hz扫到600Hz,最高频率600Hz,奈奎斯特要求采样率至少1200Hz,2000Hz留了余量。如果把fs改成800Hz,高频段会折返混叠,时频图上出现一条向下弯的假轨迹,这是初学最容易踩的坑之一。
生成之后顺手画一下时域波形,验证振幅和频率变化是否正常:
figure(1); plot(t, s); xlabel('时间 (s)'); ylabel('幅度'); title('LFM 信号时域波形'); grid on;这段没有技术含量,但值得看一眼波形对不对:LFM时域波形是等幅震荡,频率肉眼可见地变密再变疏(取决于扫频方向)。如果画出的是明显被调幅的包络,多半是相位项写错,比如漏了0.5倍,或者把t.^2写成了t*2。
这里顺带提一个工程细节:雷达系统里LFM经常用复信号(解析信号)表示,s(t) = A·exp(j·2π(f_c·t + βt²/2)),好处是瞬时频率可以直接从相位提取,不用考虑负频率。LFM.m里用实信号cos版本,适合做STFT和时频图展示;如果你后面要接脉冲压缩或匹配滤波,建议改成复信号版本,两者在STFT上的表现差别不大,但复信号方便直接算瞬时频率。
2.3 采样率、带宽与脉宽怎么配:三个常犯错误
三个参数不是独立变量,它们通过β = B/T和奈奎斯特条件相互约束。先说最常见的错误:采样率取得太随意。LFM的最高瞬时频率是fc+B,不是fc,采样率必须满足fs > 2(fc+B)。很多人只记得fs > 2fc,忽略了扫频带宽B,结果高频段混叠。判断方法很简单:画出时频图后,如果最高频率处出现一条"倒钩"状轨迹,基本就是采样率不够。
第二种错误是时宽带宽积配得不合理。BT过小(比如B=10、T=0.1,BT=1),信号几乎没有调频展宽,时频图上只是窄带的一条横线,完全看不出LFM特性。一般仿真希望看到明显的斜线轨迹,BT至少取几十以上。前面例子B=400、T=0.5,BT=200,时频图上扫频特征非常明显。
第三种错误是不检查β的量级。β=B/T,如果T很小而B很大,扫频率会高到每秒钟扫过极大带宽,STFT窗长稍长,窗内频率变化就超过一个频率分辨单元,时频图直接糊掉。参数配好后,建议先手算一下:窗长为win_len时,窗内频率变化量是β×win_len/fs。这个值最好小于两倍的频率分辨率fs/nfft,否则窗内频率变化太大,STFT结果必然模糊。比如β=800Hz/s、win_len=256、fs=2000时,窗内频率变化约102Hz,而频率分辨率约3.9Hz,差了二十多倍。这个矛盾正是后面引入D倍抽取的动机。
3. STFT时频分析:窗函数、hop与时频分辨率的取舍
3.1 从FFT到STFT:为什么LFM必须看时频
对整段LFM信号直接做FFT,得到的是幅度谱:在fc到fc+B的带宽内近似一个矩形平台,看不到"频率随时间变化"的信息,因为FFT把时间维度完全积分掉了。LFM是非平稳信号,频率是时间的函数,要观察它就必须用能同时给出时间和频率两个维度的工具。STFT做的就是这件事:把信号切成一段一段,对每一段分别做FFT,按时间顺序拼成一张二维图。
STFT公式 G(f,t) = ∫x(τ)w(t-τ)e^{-j2πfτ}dτ 里,w(t)是窗函数。窗的作用是只保留t附近的一小段信号,屏蔽远处。窗越长,参与FFT的数据越多,频率分辨率越高(Δf ≈ fs/win_len);但窗覆盖的时间范围也越大,时间分辨率越差。这是STFT与生俱来的矛盾:不可能同时获得高的时间和频率分辨率。对LFM这种频率快速变化的信号,这个矛盾尤其尖锐。
一个具体的数感:假设fs=2000Hz,窗长256点,频率分辨率约2000/256 ≈ 7.8Hz;窗的时间跨度是256/2000 = 128ms。在这128ms内,如果β=800Hz/s,频率已经变化了102Hz,相当于十几个频率分辨单元。也就是说,这个窗内根本不是"一段近似单频",而是一段扫频信号,STFT结果会拉宽成一块,而不是一条清晰的线。这就是为什么有时候时频图糊成一片,不是代码错了,是窗长和信号的调频速率根本不匹配。
3.2 窗函数与hop:时频分辨率的跷跷板
STFT里有两个参数直接决定时频图质量:窗函数和hop(帧移)。窗函数常见有矩形窗、汉宁窗、海明窗。矩形窗主瓣最窄、频率分辨率最高,但旁瓣很高,会在时频图上产生频谱泄漏造成的横向条纹;汉宁窗和海明窗旁瓣低很多,代价是主瓣略宽。处理LFM这种宽带信号,默认用汉宁窗,主瓣展宽的损失可以接受,漏出来的假条纹更烦人。
hop决定相邻两帧起点之间的间隔。hop=128(50%重叠)是常见起点,重叠越多时间轴越密,但计算量增大,而且相邻帧高度相关,并不会带来新的信息。对LFM时频图来说,hop取窗长的1/4到1/2通常就够。如果希望时频图上的轨迹更平滑,先加帧间插值,而不是一味减小hop。窗长和hop的选择顺序应该是:先定窗长满足频率分辨率需求,再定hop满足时间轴密度需求,最后检查窗内频率变化量。
三种窗的取舍大致如下:
| 窗函数 | 主瓣宽度 | 旁瓣抑制 | 适用场景 |
|---|---|---|---|
| 矩形窗 | 最窄 | 差,约-13dB | 不推荐用于LFM时频图 |
| 汉宁窗 | 较宽 | 好,约-31dB | LFM时频分析默认选择 |
| 海明窗 | 较宽 | 较好,约-41dB | 幅度精度优先时使用 |
3.3 STFT代码实现与LFM时频图的读图方法
手动实现STFT的好处是每个参数都可见,比直接调spectrogram函数更容易排查问题。LFM.m里的思路基本是下面这段:
% STFT参数 win_len = 256; % 窗长,决定频率分辨率 hop = 64; % 帧移,相邻帧时间间隔 nfft = 512; % FFT点数,不足时自动补零 win = hann(win_len, 'periodic'); % 汉宁窗,周期型更适合谱分析 num_frames = floor((N - win_len) / hop) + 1; spec = zeros(nfft/2, num_frames); for k = 1 : num_frames idx = (k-1)*hop + 1 : (k-1)*hop + win_len; seg = s(idx) .* win; % 加窗 fft_out = fft(seg, nfft); spec(:, k) = abs(fft_out(1:nfft/2)); % 只取正频率一半 end % 画时频图 f_axis = (0 : nfft/2 - 1) * fs / nfft; t_axis = (0 : num_frames - 1) * hop / fs; imagesc(t_axis, f_axis, 20*log10(spec / max(spec(:)) + eps)); axis xy; xlabel('时间 (s)'); ylabel('频率 (Hz)'); colorbar;逻辑说明:num_frames的公式保证最后一段窗口不会越过信号末尾;每一帧先取出长度为win_len的片段,与窗函数逐点相乘,再做nfft点FFT;spec只保留前nfft/2行,因为实信号的频谱关于奈奎斯特频率对称,后半段冗余。imagesc把spec显示成二维图,用20*log10转成dB,避免幅度差异大的时候弱信号看不见。
参数说明:nfft=512比win_len大,FFT会自动补零,补零不增加真实分辨率,但让频谱显示更平滑。频轴f_axis的第k个点对应频率k·fs/nfft;时间轴t_axis第k帧对应时间(k-1)·hop/fs。这里最容易错的是漏掉axis xy,不写的话图像纵轴会上下颠倒——imagesc默认把第一行画在顶部,而spec的第一行对应最低频率,所以必须翻转。
用这组参数跑出来的时频图,应该看到一条从200Hz斜向上扫到600Hz的亮线。如果看到的是断断续续的亮斑,十有八九是hop太大或者窗长太短。另外,20*log10里加eps防止log(0),不然动态范围压到-Inf会显示成一片黑。
4. D倍抽取:时频分辨率改进的具体做法与边界
4.1 D倍抽取的思路:两种做法与选择理由
上一章末尾提到,窗内频率变化量可能超过频率分辨率一个数量级,时频轨迹被拉宽成块。一个直接的改进思路:降低信号的有效采样率,让同样的窗长覆盖更长的时间,频率变化在窗内更平缓。D倍抽取(降采样)就是这个思路的工程实现——对信号每隔D个点取一个,等效采样率从fs降到fs/D,信号总时长不变,样本点数变成N/D。再对抽取后的信号做STFT,同样的窗长覆盖的实际时间跨度扩大D倍,窗内频率变化量相应减小,频率轨迹更清晰。
但抽取不是免费的:降采样后奈奎斯特频率变成fs/(2D),原信号里超过这个上限的高频成分会折叠回低频区,产生混叠。所以抽取前必须先做低通滤波,把fs/(2D)以上的分量滤掉,否则时频图上会出现根本不属于原始信号的假轨迹。这是D倍抽取最容易翻车的点。
另一种做法是在STFT输出端做抽取:对时频矩阵按时间轴每D帧取一帧,等效把hop放大D倍。这种做法不改变原始数据,计算量小、实现简单,LFM.m里采用的就是这个思路。它提升的不是物理分辨率,而是让显示出来的时频网格更干净,同时去掉高度冗余的相邻帧。两种做法本质上是两套东西:前者真的在改变时频分析的分辨率特性,后者是显示层面的整理。搞清楚代码里用的是哪一种,才不会改错地方。
4.2 LFM.m里D倍抽取的实现与参数说明
按摘要描述,LFM.m把D倍抽取放在STFT之后、可视化之前,处理对象是时频矩阵spec。常见实现如下:
% D倍抽取:对STFT输出时频矩阵沿时间轴抽取 D = 4; % 抽取因子 spec_d = spec(:, 1:D:end); % 每D帧保留1帧 t_axis_d = t_axis(1:D:end); % 时间轴同步抽取 figure(2); imagesc(t_axis_d, f_axis, 20*log10(spec_d / max(spec_d(:)) + eps)); axis xy; xlabel('时间 (s)'); ylabel('频率 (Hz)'); title(sprintf('LFM 时频图(D=%d)', D)); colorbar;逻辑说明:spec(:, 1:D:end)是对矩阵做列方向索引,保留第1、1+D、1+2D…列,也就是每隔D帧保留一帧,行数(频率维)不动,因此频轴f_axis保持不变。时间轴按同样的索引方式抽取,保证横坐标和保留的帧一一对应。抽取后帧数从num_frames变成约num_frames/D,显示密度降低,但每一帧仍是完整的加窗FFT结果,幅度信息没有损失。
参数说明:D=4意味着显示帧数是原来的四分之一,相邻保留帧的时间间隔从hop/fs变成D·hop/fs。以fs=2000、hop=64为例,原始帧间隔32ms,抽取后变成128ms,肉眼看起来轨迹更稀疏但更清晰。D取2到8是比较常用的范围。D=1等于不抽取,D超过16之后,稀疏的帧会让LFM的斜线轨迹变成离散点列,反而看不出扫频特征。
如果要在信号路径上做抽取(前面说的第一种做法),MATLAB里用resample一步完成滤波和抽取:
% 可选:先对信号做D倍抽取再STFT(另一种实现) D_sig = 4; s_dec = resample(s, 1, D_sig); % 内部自带抗混叠低通滤波 fs_dec = fs / D_sig; % 等效采样率要同步修改这段在LFM.m里一般作为注释块出现,因为它改变的是整个STFT的输入,使用时要同步修改fs、fc、B等参数的对应关系,改动面比较大。初学建议先读懂主路径的矩阵抽取,再研究信号抽取。
4.3 D值怎么选:混叠风险与轨迹完整性
选择D没有解析公式,更多是看图说话。判断标准有两条:一是保留帧的时频轨迹是否连续,LFM的瞬时频率是线性函数,时频图上的亮轨迹应该是连续斜线,如果D太大,帧与帧之间间隔过大,轨迹会呈现"楼梯"甚至"虚线"状;二是抽取后是否出现本不该有的频率分量,这主要发生在信号路径抽取的场景。
实际调试套路:先D=1跑一遍,记住时频图的样子;然后D=2、D=4、D=8各跑一遍,对比轨迹清晰度和连续性。D=4通常是一个甜点值——既能把STFT帧间冗余消掉,又不会让轨迹稀疏到丢失扫频趋势。如果D=4仍然模糊,问题大概率不在抽取,而在窗长和hop的搭配,回去重新调STFT参数更有效。
还有一点需要提醒:矩阵抽取和信号抽取是两回事。如果代码里信号抽取和矩阵抽取同时存在,等效的显示帧间隔是两者的乘积(D_sig × D × hop/fs),这个值一旦超过窗长时间,时频图的轨迹就会明显断裂。所以要么只用一种抽取,要么把两个D乘起来作为总的D来评估。
5. 避坑指南:LFM时频仿真最容易翻车的五个点
时频图这玩意,跑通容易跑准难。下面五个坑我都在调试LFM时频仿真时踩过,每一条都按"现象→原因→解决"写清楚,照着排查能省不少时间。
5.1 频率轴标定错误:时频图的频率坐标对不上理论值
现象:时频图跑出来了,斜线形状也对,但起点不是200Hz、终点也不是600Hz,整体往上或往下偏。
原因:f_axis生成公式与spec的行数不匹配。spec取了fft结果的前nfft/2行,对应的频率轴应该是 (0:nfft/2-1)·fs/nfft,有人直接用 (0:nfft-1)·fs/nfft,矩阵只有256行,频率轴却给了512个点,imagesc把两者按位置硬对齐,频率坐标整体错位。
解决:养成一个习惯,画图前先断言length(f_axis) == size(spec, 1)。调试时在时频图上叠一条理论瞬时频率线f_theory = fc + beta·t_axis,如果斜线和理论线不重合,说明不是画图问题就是参数本身设错了。
5.2 窗长过短:时频图糊成一片
现象:清晰斜线看不到,只有一块模糊亮带,扫频起点和终点都分不清。
原因:win_len太短,频率分辨率Δf = fs/win_len远大于LFM瞬时频率的变化速率,窗口内信号已经不是单频分量,FFT把一小段扫频谱拉平了。
解决:win_len往大调,nfft跟着调。一个可用的经验值:Δf至少小于B/20。以B=400为例,Δf应小于20Hz,fs=2000时win_len至少要100点,实际用256点更稳。注意win_len加大后,窗内频率变化量β×win_len/fs也会增大,需要同步评估,必要时配合D倍抽取使用。
5.3 D倍抽取过头:轨迹断成虚线
现象:D=16或更大时,时频图上的斜线变成一串孤立的亮点,扫频方向勉强能看出,但连续性没了。
原因:抽取后相邻保留帧的时间间隔D·hop/fs过大,LFM瞬时频率在相邻两帧之间的变化量超过了频率分辨率,时频轨迹接不上。等效来看,是时间采样太稀疏,信号本身的连续性在显示层面丢失了。
解决:把D降到4左右,或者先把hop减小再抽取。判断标准是β × D × hop / fs这个频率变化量最好小于Δf,这样轨迹才连续。如果D=4仍然断线,优先调hop而不是继续加D。
5.4 首尾边界效应:时频图两端出现异常亮线
现象:时频图最左边和最右边各有一道竖直亮纹,幅度比主体还高,把动态范围都占满了。
原因:信号起点和终点是硬截断,窗函数在边界处只覆盖了部分信号,幅度跳变产生宽频谱泄漏。STFT第一帧和最后一帧的窗内数据不完整,谱线被展宽成竖直亮纹。
解决:显示时把首尾各几帧裁掉;更稳的做法是生成信号时加一点淡入淡出,或者对信号做边缘延拓再STFT。LFM本身是等幅信号,边界跳变主要影响显示,不是信号本质问题。
5.5 dB动态范围没压住:弱信号被淹没
现象:图整体很亮,斜线旁边的小结构完全看不见,色彩条范围被最大幅度占满,背景噪声一片白。
原因:直接用imagesc(abs(spec))画线性幅度,动态范围跨几个数量级,低幅度细节全部被压缩到同一灰度。
解决:至少用20*log10(spec/max(spec(:)) + eps),再配合clims = [-40 0]之类限制显示范围。设置色彩范围后,时频图的背景噪声和弱信号都能看到。这个习惯对所有谱图类可视化都适用。
6. 验证时频图可信度:瞬时频率提取与理论曲线对比的最后一公里
6.1 峰值提取法:把时频图变回一条频率曲线
时频图是给人看的,机器验证还需要把它变成可对比的数据。对LFM来说,最好的验证方式就是从时频矩阵里提取每帧的峰值频率,和理论瞬时频率f_theory = fc + βt直接对比:
% 从D倍抽取后的时频矩阵提取瞬时频率 [~, idx_max] = max(spec_d, [], 1); % 每帧最大幅度的行索引 f_est = f_axis(idx_max); % 映射为频率 % 理论瞬时频率 f_theory = fc + beta * t_axis_d; % 对比 figure(3); plot(t_axis_d, f_est, 'b.', t_axis_d, f_theory, 'r-', 'LineWidth', 1); xlabel('时间 (s)'); ylabel('频率 (Hz)'); legend('STFT峰值提取', '理论瞬时频率'); grid on; % 误差统计 err = f_est - f_theory; fprintf('峰值提取误差: mean=%.2f Hz, std=%.2f Hz\n', mean(err), std(err));逻辑说明:max函数沿频率维找每帧幅度最大的位置,idx_max是行索引,再通过f_axis映射成实际频率。LFM时频图的主体就是瞬时频率轨迹,峰值位置应当接近理论值。误差的均值反映系统偏差(比如频率轴标定错误),标准差反映随机抖动(主要是窗长和噪声带来的)。
参数说明:这段验证代码必须在D倍抽取之后、显示参数确定之后运行,因为f_est的横坐标要和t_axis_d对齐。如果误差均值的绝对值超过频率分辨率Δf,基本可以断定标定或参数有问题,回头查第5章的坑。
6.2 一个可以带走的调参习惯
我个人养成的习惯是:任何LFM时频图,不直接信眼睛,先叠加理论瞬时频率曲线,再看误差统计。两张图叠在一起,误差均值能压到Δf的1/3以下,我才认为这组参数可信。从那以后,每次跑LFM仿真我都强制走一遍这个验证流程——先生成信号、再STFT、再D倍抽取、最后峰值提取对比,一条路走完参数对不对就清楚了。这套流程看着多几步,但能避免"图看着差不多、实际参数全错"的尴尬。希望帮到你。
本文还有配套的精品资源,点击获取