1. 从一张平均过的频谱说起雷达信号为什么非得换种看法接手雷达信号时频分析的头一个月我曾在频谱分析仪前坐了整整半天示波器上明明是两段完全不同的脉冲——一段频率从低频往高频扫另一段干脆在高频和低频之间来回跳FFT却给出了几乎一模一样的宽馒头谱线。那一刻我意识到傅里叶变换只回答了有哪些频率成分却把这些频率什么时候出现的信息彻底抹掉了。对雷达信号这种典型的非平稳信号光看频谱等于把所有时间切片叠在一起洗牌。后来我把小波变换接进了 MATLAB 工作流才真正看清每个脉冲内部的频率细节。这篇文章就围绕小波变换 雷达信号 MATLAB展开原理、程序、参数坑一次说透适合正在做信号处理、做电子侦察或雷达仿真、以及准备把时频分析工具应用在工程里的同学参考。1.1 FFT 把时间信息藏起来了平稳信号还好说频率成分在整个观察时间内一成不变FFT 的结果就是一段历史的忠实记录。可雷达信号几乎都是非平稳的线性调频LFM脉冲的频率在脉宽内线性扫过几十到几百兆赫频率捷变雷达的载波在不同脉冲间跳变脉压雷达的一串脉冲内部还有幅度调制和相位编码。你拿整段信号做 FFT相当于要求现在正在发生的频率和一秒前已经结束的频率同时存在于同一个频谱里这本身就违背了物理直觉。一个很直观的类比把运动物体的频闪照片叠印在同一张底片上剑桥画出来就是一团模糊的轮廓你看不出运动员是先快后慢还是先慢后快。FFT 干的就是这件事——它把每一个时刻的频率证据全部叠加平均时间维度彻底丢失。我在实测里遇到最典型的情况是两种不同斜率的 LFM 脉冲FFT 长几乎一样都是一片宽带宽的隆起只是中心频率和 3dB 带宽略有差异。这时候你想通过频谱区分它们是哪类调制、有没有跳频规律基本是瞎猜。1.2 短时傅里叶变换能救场但有一道卡死的墙既然全局 FFT 不行最自然的想法是把信号切成小段逐段做 FFT这就是短时傅里叶变换STFTX(τ, f) ∫ x(t) w(t − τ) e^{−j2πf t} dt问题在于窗长 w 是固定不变的。窗短时间定位精细但频率分辨率差窗长频率分辨率高但时间分辨率惨不忍睹。这是海森堡不确定性原理在时频分析里的具体表现——时频联合分辨率存在下限你没法让时间分辨率和频率分辨率同时无限提高。对雷达 LFM 脉冲来说这道坎非常致命。我常用的一组仿真参数是脉宽 0.2 秒、带宽 150 Hz为演示我把频率降到了基带可处理的量级工程上同理换算。要分辨 2 Hz 的频率差STFT 窗长至少需要 0.5 秒这已经比整个脉冲还长了可脉冲内部的频率变化要求时间分辨率达到毫秒级。二者不可兼得所以 STFT 图上的 LFM 脊线要么糊成一条宽带要么被切得断断续续。这时候就需要一把会自己伸缩的尺子——小波变换。2. 小波变换的思路一把会自己伸缩的尺子2.1 CWT 的数学骨架信号与小波做内积试穿连续小波变换CWT的定义长这样CWT_x(a, τ) (1/√a) ∫ x(t) ψ*((t − τ)/a) dt其中 ψ(t) 是母小波a 是尺度因子τ 是时间平移量星号代表共轭。直观理解就是拿一个波形模板把它压扁或者拉伸然后沿着时间轴滑动看它和信号的贴合程度——贴合得越狠系数越大说明这一时刻存在与该模板形状和频率都匹配的成分。尺度 a 小小波被压缩持续时间短、震荡频率高擅长捕捉高频细节和瞬态尺度 a 大小波被拉伸持续时间长、震荡频率低擅长把握低频轮廓。这正是雷达信号最需要的频率高时自动用短窗保证时间精度频率低时自动用长窗保证频率精度。对比 STFT 的区别一句话就能说明白STFT 是拿固定焦距的镜头拍整条胶片小波是变焦镜头——高频区域自动对焦到局部低频区域自动拉远看全局。2.2 尺度到频率的换算伪频率到底伪在哪里小波的横轴是尺度 a要把它变成频率轴工程上用的是伪频率公式f Fc × fs / aFc 是小波基函数的名义中心频率。对常用的复数 Morlet 小波MATLAB 里的cmor1-1Fc 1对默认的 Morse 小波中心频率在某个固定值附近具体以程序输出为准。代入示例fs 1000 HzFc 1a 100 时对应的伪频率就是 10 Hz。尺度 a伪频率 HzFc1, fs100025005200101002050502010010为什么叫伪频率因为小波不是纯正弦波它是一段带着包络的震荡所谓频率其实是这一尺度下最主导的局部震荡频率只有对窄带信号才严格成立。但工程上画时频图、提取脊线、看扫频规律这个精度完全够用。新版 MATLAB 的cwt(x, fs)会直接返回物理频率向量 f不再逼你手动换算这个我在第 3 节展开。2.3 复数解析小波 vs 实数正交小波雷达信号该选哪边写代码之前得先分清小波的两大流派。雷达时频分析几乎都用复数解析小波比如 Morse 小波、复数 Morlet 小波cmor、解析 Morlet 小波amor。这类小波给出的是复系数幅值反映信号的局部强度相位则携带瞬时频率的细节非常适合提取 LFM 脊线、估计瞬时频率、判断相位调制。另一类是在图像处理、数据压缩和降噪里非常火的实数正交小波db4、sym8 等配合离散小波变换 DWT 使用。热词里经常看到小波变换图像去噪小波变换图像增强 python那些场景用的基本都是这套。但做雷达信号的窄带时频分析时拿 db4 直接看时频图会很难受——实数小波系数正负震荡不直接对应窄带频谱幅度判读难度大。我的习惯是两者搭配先用 db4 / sym8 这类实数小波做一遍去噪预处理压掉底噪和毛刺再用复数小波做时频分析效果比直接拿原始带噪信号出图干净得多。3. MATLAB 实战从造信号到看清时频图3.1 先造一管雷达样式的 LFM 脉冲串没有实测数据时仿真是检验算法的第一步。下面这段代码生成一串线性调频脉冲每个脉冲内部频率从 100 Hz 扫到 250 Hz脉冲宽度 0.2 秒重复间隔 0.5 秒并把复基带接收噪声叠加上去clear; clc; close all; fs 1000; % 采样率 Hz T 2; % 总时长 s t (0:1/fs:T-1/fs).; % 时间列向量 f0 100; % LFM 起始频率 Hz B 150; % 带宽 Hz pw 0.2; % 脉冲宽度 s pri 0.5; % 脉冲重复间隔 s x zeros(size(t)); for k 0:floor((T-pw)/pri) idx (t k*pri) (t k*pri pw); tau t(idx) - k*pri; x(idx) exp(1j*2*pi*(f0*tau (B/(2*pw))*tau.^2)); end x x 0.05*(randn(size(t)) 1j*randn(size(t)));为什么用复信号而不是实信号雷达接收机解调后基本都是 I/Q 复基带保留复数意味着同时保留幅度和相位后面做脊线提取、估计瞬时频率都要用相位信息。用实信号也能算但负频率会多出半个谱判读时看着别扭。3.2 路线 A内置 cwt 三行出图先看个大概调试时我的第一个动作永远是直接画figure; cwt(x, fs);一行代码默认用 Morse 解析小波自动选尺度自动画影响锥COIY 轴直接就是物理频率Hz。跑完之后你会立刻看到四条清晰的斜线——那就是四个 LFM 脉冲在时频平面上的轨迹斜率就是调频斜率 B/pw。如果想自己控制后续的定制绘图就取系数[wt, f] cwt(x, fs); figure; P 10*log10(abs(wt).^2 eps); % 转 dB方便看动态范围 P P - max(P(:)); % 归一化到 0 dB imagesc(t, f, P); axis xy; % 让 Y 轴从小到大 xlabel(时间 / s); ylabel(频率 / Hz); colorbar;cwt在不同 MATLAB 版本里语法略有差异老版本R2016a 以前是cwt(x, scales, wname)返回系数矩阵新版本是cwt(x, fs)直接给你物理频率。我建议有新版就用新版少走换算的弯路。3.3 路线 B手写一个最小 CWT把公式落成代码只看内置函数容易把原理忘掉。我面试新人时经常让他们手写一个 CWT这里给一个教学版实现——通过卷积完成内积运算逻辑和公式一一对应function cfs manual_cwt(x, fs, scales, fc, fb) % 教学版连续小波变换用卷积实现 % x - 复数信号 % fs - 采样率 Hz % scales - 尺度序列 % fc - 复数Morlet中心频率默认1 % fb - 带宽参数默认1 x x(:).; n length(x); cfs zeros(length(scales), n); halfwidth 8; % 小波支撑半宽u∈[-8,8] for k 1:length(scales) a scales(k); m -round(halfwidth*a):round(halfwidth*a); u m / a; % 归一化时间 psi (1/sqrt(pi*fb)) * exp(1j*2*pi*fc*u) .* exp(-u.^2/fb); h conj(fliplr(psi)) / sqrt(a); % 卷积核翻转共轭能量归一 cfs(k,:) conv(x, h, same); end end调用targetF 20:5:450; % 想看的频率范围 Hz scales fc * fs ./ targetF; % 反推尺度这里 fc1 cfs manual_cwt(x, fs, scales, 1, 1); figure; imagesc(t, targetF, abs(cfs)); axis xy; xlabel(时间 / s); ylabel(频率 / Hz);注意这个教学版没有做完整的边界处理和精确归一化常数因子不会改变时频图的形态但工程级的参数测量请以内置cwt为准。我坚持手写一遍的原因很简单只有亲手落过一次卷积你才会理解为什么尺度小的区域时间分辨率高、尺度大的区域频率分辨率高也才能在后续遇到频率轴偏移时迅速定位问题。3.4 时频图显示细节dB、色标裁剪与脊线提取雷达时频分析最终是给人看的显示效果直接影响判读。三项基本功我每次都做。第一转 dB 并做百分位裁剪。原始系数动态范围太广最强脉冲会把色标整个拉满弱信号全变底色。用百分位锁定色标区间PdB 10*log10(abs(wt).^2 eps); PdB PdB - max(PdB(:)); lo prctile(PdB(:), 3); hi prctile(PdB(:), 97); imagesc(t, f, PdB, [lo hi]); axis xy;第二提脊线估计瞬时频率。雷达 LFM 信号的脊线就是时频图上能量最大的那条线逐列找abs(wt)的最大值位置再映射回频率[~, ridx] max(abs(wt), [], 1); ifreq f(ridx); hold on; plot(t, ifreq, r-, LineWidth, 1.5);对一条干净的 LFM 脉冲提出来的ifreq应该是严格线性上升的用polyfit(t, ifreq, 1)就能拟合出调频斜率直接和仿真参数核对。第三时刻记住时间和示波器波形对照。不要把时频图当成孤立图片看同一时间段里信号时域波形的有无、包络强弱、时频图上的能量分布应该相互印证。我在第 5 节会给出一个具体的对照检查方法。4. 参数选不好时频图就是一团墨小波基、频率范围与 COI4.1 小波基怎么挑先确认你要提取什么特征内置cwt默认的 Morse 小波在时间定位和频率定位之间平衡很好是雷达信号的默认首选。我一般先用它出图如果发现 LFM 脊线有锯齿状阶梯就把voicesperoctave调高让尺度轴更细[wt, f] cwt(x, fs, voicesperoctave, 16);voicesperoctave表示每个倍频程里划分多少条尺度网格。默认是 10调成 16 之后时频图上的脊线明显更顺滑代价是计算量线性上升。对 2 秒、1000 Hz 采样的信号16 几乎没有体感但对长记录、高采样率的实测数据多试几个值再定。复数 Morlet 小波cmor的优点是公式直观老代码里常配scal2frq手动换算频率。它的带宽参数 fb 和中心频率 fc 可以单独调fb 越小频率选择性越好但时间分辨率变差。我的经验雷达 LFM 分析用cmor1-1.5这类中档带宽足够调试起来比 Morse 更容易看到参数变化的直接效果。选基的判断标准不是哪个更高级而是结论是否稳定。同一个信号换三种小波Morse、cmor、amor都出图如果脊线位置基本一致说明结果可信如果换一个小波图就面目全非那大概率是参数设置或边界问题而不是小波本身的问题。4.2 频率范围怎么定从雷达脉冲反推而不是让程序默认内置cwt的默认频率范围会横跨能算的所有频率但这不一定匹配你的信号。LFM 脉冲中心 100 Hz、带宽 150 Hz 时能量集中在 100~250 Hz如果默认图把 0~500 Hz 全铺开脊线会细得像发丝调频斜率肉眼根本看不清。正确做法是主动限定频率范围。新版 MATLAB 可以直接传[wt, f] cwt(x, fs, FrequencyLimits, [30 450]);这里 30 Hz 留出低频余量450 Hz 略低于奈奎斯特频率 500 Hz。用手写 CWT 时反向反推尺度targetF 30:5:450; % 想要的频率网格 scales 1 * fs ./ targetF; % fc1频率上限不能超过 fs/2否则混叠下限也别太低因为一个 2 秒的信号里 0.5 Hz 只有 1 个完整周期属于不可信的极低频。雷达信号一般也不用看直流附近先做一步去直流更稳x x - mean(x);不去直流的话0 Hz 附近一条横贯全图的高能量带会把色标拉满弱目标回波全部被压成背景色。4.3 边界效应与 COI雷达短脉冲的信任边界小波有支撑长度当小波窗滑动到信号边界附近时有一部分伸出信号外面系数就不可靠。影响锥Cone of Influence, COI标出的正是这种不可信区域。内置cwt绘图的默认界面里COI 以半透明阴影标出老版本还会画两条对称的曲线。雷达信号的边界效应比普通音频信号更麻烦——脉冲可能正好落在记录开头或结尾。我遇到过把强信号在 COI 内形成的亮斑当成新目标的教训。排查方法很简单把时间轴挪到脉冲完全位于信号中段的位置重算一次如果亮斑位置和形状都变了那基本就是边界效应。工程上我的处理习惯是给信号两端补一段零值pad做完 CWT 后再裁掉对应区域相当于把不可信边缘推出显示区域。实测下来在保留原始信息的前提下抬高了边缘结果的可信度。注意 padding 本身会引入人为的截断所以 pad 长度不要超过 10% 的总时长而且显示时一定得裁掉。5. 三次翻车复盘雷达信号小波分析里的典型坑5.1 翻车一默认一把梭频率范围太宽LFM 脊线被压成细线第一次跑通程序时我兴冲冲地直接cwt(x, fs)出来的图确实有时频结构但四条 LFM 脊线细得像头发丝调频斜率完全看不出换成不同仿真的脉冲对比时也抓不到差异。原因就是默认频率范围 0~500 Hz 全铺满了信号能量集中在 100~250 Hz 这一段显示时被压缩到很小的高度范围。修复就是第 4.2 节的频率限定把FrequencyLimits设成 [30 450]再把voicesperoctave提到 16。重画之后四条平行的斜线从 100 Hz 爬到 250 Hz端点清清楚楚四个脉冲的时序、脉宽、斜率一次看明白。从那以后我养成了先查信号频谱、再定时频图频率范围的习惯绝不让程序替我做这个决定。5.2 翻车二手动换算尺度时频率轴差了一个数量级有一次我用老版cwt语法配合scal2frq手动换算频率轴画出来的 LFM 脊线从 10 Hz 一直扫到 25 Hz可仿真参数明明写的是 100~250 Hz差了整整十倍。排查半天才发现问题出在中心频率上——我当时用了一个默认 Morse 小波的中心频率值当 1 算实际不是 1换算公式f Fc*fs/a里的 Fc 没对齐。从那以后我总结了两个铁律。第一凡是怀疑频率轴有问题先扔一个单频探针进去验证ftest 150; probe exp(1j*2*pi*ftest*t); [wtp, fp] cwt(probe, fs); [~, mi] max(sum(abs(wtp), 2)); % 全行能量最大处 fprintf(检测到的频率峰%.2f Hz\n, fp(mi));如果输出不是 150说明频率标定有问题。第二做手写 CWT 时只用cmor1-1这类 Fc 确实为 1 的小波堵死换算错误。单频探针这个小技巧成本极低应该成为每次时频分析首跑之前的固定动作。5.3 翻车三强脉冲加色标饱和弱目标回波被涂成一片红雷达场景经常是几个强脉冲旁边藏一个弱脉冲或者强直达波盖住弱反射回波。直接画abs(wt)的线性色标时强脉冲亮黄色一片弱回波在色标映射下全部落在最低端和底色完全分不出来我一度以为仿真没有产出第二个脉冲。转 dB 是不够的线性转对数只改变动态范围色标上限依然被最强点拽死。必须做第 3.4 节的百分位裁剪把色标下限定在第 3 百分位、上限定在第 97 百分位弱回波对应的能量才能映射到可见颜色区间。裁剪之后的图里强脉冲不是亮到刺眼弱脉冲也不是黑成一片整体结构清楚得多。另外别忘了去直流0 Hz 高能量带对色标的破坏比强脉冲还狠。5.4 快速检查清单我贴在显示器边上的小抄检查项方法常见后果频率轴是否正确单频探针验证最高峰位置脊线整体偏移可能差一个中心频率倍数色标是否饱和观察弱回波是否可见弱脉冲、弱分量消失目标是否进入 COI对照内置绘图阴影区把边缘亮斑当成真实信号采样率是否撑得住确认最高关注频率小于 fs/2高频段混叠、脊线回折直流是否去除看 0 Hz 附近是否有横贯亮带0 Hz 强带压掉其他细节时间轴是否对齐与示波器包络上下对照脉冲位置错位、脉宽判断错误这六项我基本每次算新数据都要过一遍其中单频探针和百分位裁剪几乎可以说是雷达信号时频分析的保命动作。最后说一句体会。小波时频分析不是万金油——单频连续波信号FFT 几秒钟就给你答案没必要上小波多分量信号交叉严重的时候Wigner-Ville 分布其实有交叉项的麻烦而小波是线性变换天然没有交叉项干扰这才是它在我雷达脉冲分析工作流里长期占位置的原因。我的做法是 STFT 先粗看全局小波再细看脊线和调制细节两者配合比你单独押注任何一套都稳。如果你也正在做雷达信号的时频分析建议从今天这串 LFM 仿真代码开始先跑通图像再动手提脊线、估斜率——整套流程走一遍你会比看十篇原理文章都明白得快。
