简介频谱分析与时频分析是通信与信号处理领域的基础工具本资源面向电子工程、通信技术方向的学习者与研发人员聚焦调制解调信号观测、时频谱图绘制与频谱特征提取等常见需求适合在课程实验或项目仿真中参考使用。资源共3个文件均为MATLAB脚本.m压缩包仅2KB体积小巧但代码结构集中覆盖信号产生、低通等效变换与包络相位计算等关键环节可直接用于快速验证频谱和时频分析方法。目前已有717人学习下载。通过研读这些源码读者可以掌握如何构造分析信号、调用时频变换函数、绘制时频谱图以及解读频谱分布同时对照描述中的调制解调、STFT、功率谱密度等概念加深对非平稳信号处理流程的理解并为后续通信系统设计与故障检测分析提供可复用的脚本基础。1. 频谱分析不是只有一个 FFT从三个 m 文件看一条完整信号链拿到一套只有main.m、loweq.m、env_phas.m三个脚本的 MATLAB 工程时我一开始也以为这只是个画频谱图的演示程序。跑完才发现这三个文件合起来实际上是一台简化版的软件无线电分析仪loweq.m负责把基带信息“搬”到等效低通信号上main.m完成波形生成和频谱测量env_phas.m再从解析信号里把包络、瞬时相位和瞬时频率一一抽出来。换句话说它覆盖了调制、传输、解调、频域测量和时频图绘制的完整闭环。对做通信物理层、信号处理和嵌入式 FFT 分析的人来说这套代码的价值不在于某一个算法多高级而在于它把“频谱分析”从单个 FFT 函数扩展成了一条可验证的链路。想搞清楚非平稳信号怎么观察、调制参数怎么定、解调误差从哪里来顺着这三个脚本走一遍比读十页公式更直观。2. 调制端建模用 loweq.m 把等效低通和真实载频解耦2.1 为什么要用“等效低通”而不是直接生成载波很多人在 MATLAB 里仿真调制信号时习惯直接写y cos(2*pi*fc*t phi)把载频fc取到几百 MHz采样率跟着拉到 GHz 量级。这样做的直接后果是一个 1 秒的仿真就要存几十 GB 数据FFT 点数大得离谱而真正有用的信息带宽往往只有几 kHz。纯属浪费算力。通信系统仿真里更常见的做法是**等效低通complex baseband / equivalent lowpass**模型。原理很简单实带通信号x(t) Re{s(t) * exp(j*2*pi*fc*t)}其中s(t)是复包络。只要把复包络算出来载频只是一个乘法因子甚至可以在分析阶段完全忽略。这样一来采样率只需要覆盖信号带宽而不是覆盖载频仿真数据量直接下降两三个数量级。这也是 GMSK 这类恒包络调制能高效仿真的原因。GMSK 调制解调原理看着复杂本质就是在 IQ 平面上做相位累积再让相位路径经过高斯滤波变得平滑从而压缩频谱旁瓣。用等效低通模型来写核心代码不会超过二十行。2.2 loweq.m 的典型实现与参数含义下面给出一份和资源中loweq.m功能等价的实现。它把 AM 和 FM 两种调制方式统一在同一个函数里输出复基带信号。function [s, t] loweq(mode, info, Fs, fm, fc, mf) % [s, t] loweq(mode, info, Fs, fm, fc, mf) % mode : am 幅度调制 / fm 频率调制 % info : 基带信息序列单声道列向量 % Fs : 采样率必须大于 2*fm这里 fm 是信息最高频率 % fm : 信息信号最高频率用于定时间轴和成型 % fc : 等效载频只影响相位旋转不参与采样率设定 % mf : 调制度AM 用调幅深度FM 用频偏系数 Ts 1/Fs; N length(info); t (0:N-1). * Ts; switch lower(mode) case am % 直流分量 1 保证包络不会过零mf 取 0~1 之间 s (1 mf * info) .* exp(1j*2*pi*fc*t); case fm % 相位是信息的积分mf 等效为频偏/信息频率 phase 2*pi*mf*fm * cumsum(info) * Ts; s exp(1j * (phase 2*pi*fc*t)); otherwise error(loweq: unknown mode %s, mode); end end这段代码的关键点有三个。第一AM 模式里exp(1j*2*pi*fc*t)是未被显式调用的复载波实际信号是real(s)相当于把载波搬到了fc处。第二FM 模式用cumsum(info) * Ts做离散积分这正是角度调制里“相位是调制信号的积分”这一结论的直接翻译。第三Fs只由info的带宽决定和fc没有关系这是等效低通建模最反直觉也最容易被忽略的地方。参数调优时mf的选择直接影响解调难度。AM 的mf超过 1 会导致过调幅包络出现截断包络检波后谐波失真明显FM 的mf决定最大频偏频偏过大则瞬时频率可能超出分析带宽时频谱图上会出现折叠。我在调试时一般先把mf设为 0.5 做 AM、1.0 做 FM在确认整条链路无误后再加大。2.3 验证调制器输出是否正确写完调制函数后先做两个快速检查一是看频谱形状二是看瞬时频率是否符合预期。Fs 100e3; fm 1e3; N 8192; t (0:N-1)./Fs; info sin(2*pi*fm*t); % 1 kHz 单音信息 s_am loweq(am, info, Fs, fm, 10e3, 0.5); s_fm loweq(fm, info, Fs, fm, 10e3, 1.0); % 频谱检查AM 应出现载波 ± 1kHz 两根边带FM 的谱线间隔等于 fm f_axis (0:N-1)/N*Fs - Fs/2; figure; subplot(2,1,1); plot(f_axis/1e3, fftshift(abs(fft(s_am)))); title(AM spectrum); xlim([-5 20]); subplot(2,1,2); plot(f_axis/1e3, fftshift(abs(fft(s_fm)))); title(FM spectrum); xlim([-5 20]);看 AM 频谱时理想情况下载波谱线在 10 kHz 处两侧 9 kHz 和 11 kHz 各有一根对称边带幅度是载波的一半。FM 谱相对复杂但谱线间隔严格等于fm且高次边带随阶数衰减。如果谱线间隔不对先检查Fs和fm的量纲是否匹配如果出现非预期的大块噪声底则说明cumsum的累积误差在起作用可以考虑改用filter(1, [1 -1], info)做抗饱和积分。3. main.m 里的频谱与极坐标图FFT 参数、频率轴标定和时频谱图的边界效应3.1 main.m 的主流程拆解main.m在整个工程里承担的角色是“总装车间”一般包含五个步骤设置仿真参数、调用loweq生成信号、叠加噪声、进行频谱分析、绘制结果。其中频谱分析部分写得是否严谨直接影响后续对调制方式的判断。很多人在这一步把fft()的结果直接plot然后对着横轴发呆因为横轴是“点数”而不是“频率”。正确的频率轴标定公式是f (0:N-1) * Fs / N对应到单边谱则是f (0:N/2-1) * Fs / N。采样率是频率轴唯一的“定标尺”不知道采样率FFT 只能给出归一化频率0 到 1 对应 0 到 Fs这也是“不知道采样率怎么求频率频谱”这个问题最常见的场景——工程数据一旦丢了采样率频域分析基本只能看形状无法读数。3.2 频率分辨率、补零和窗函数频谱分析的核心指标是频率分辨率Δf Fs / N。注意补零能改变N从而细化频点间隔但它不能提高真正的分辨率因为有效信息长度没有变。想分辨两根间隔 100 Hz 的谱线FFT 点数至少要大于Fs / 100补零只是让谱线看起来更圆润。Fs 100e3; N 10000; t (0:N-1)./Fs; x sin(2*pi*1000*t) 0.5*sin(2*pi*1020*t); % 间隔 20 Hz 的两根谱线 % 直接 FFT矩形窗 X1 fft(x); % 加 hamming 窗后再 FFT w hamming(N); X2 fft(x .* w); % 频率轴单边谱 f_axis (0:N/2-1) * Fs / N; figure; subplot(2,1,1); plot(f_axis, 20*log10(abs(X1(1:N/2))/N)); title(Rectangular window); xlim([900 1100]); subplot(2,1,2); plot(f_axis, 20*log10(abs(X2(1:N/2))/N)); title(Hamming window); xlim([900 1100]);这段代码直接体现了窗函数的作用矩形窗下 1000 Hz 和 1020 Hz 两根谱线大概率粘连成一个峰因为 20 Hz 的间隔小于Δf Fs/N 10 Hz。加 Hamming 窗后主瓣变宽旁瓣抑制增强但 20 Hz 间隔在 10 Hz 分辨率下依然不一定能分开只是不会再被旁瓣淹没。参数选择上hamming(N)的等效噪声带宽比矩形窗宽约 1.36 倍所以幅度读数要按窗的相干增益修正否则测出来的功率偏低。3.3 从全局 FFT 到时频谱图分段 FFT 的正确姿势全局 FFT 适合平稳信号换成跳频、扫频或调频信号就失效了。观察这类非平稳信号标准做法是短时傅里叶变换STFT把信号切成一小段一小段每段做 FFT再把结果按时间顺序拼成二维矩阵。这就是“分段频谱图”的实质也就是时频分析中时频谱图的来源。分段 FFT 有两个参数需要手工调窗口长度wlen和重叠率overlap。窗口长度决定时间分辨率和频率分辨率的折中窗越长频率分辨率越高但时间上越模糊重叠率越高时频谱图越平滑但计算量越大。MATLAB 自带的spectrogram函数可以直接用但用矩阵切片手写一段会更清楚每一步在做什么。function [S, f, t] my_stft(x, Fs, wlen, overlap) % 手动实现短时傅里叶变换 % wlen : 窗长度点数 % overlap: 重叠率0~0.99 hop round(wlen * (1 - overlap)); nfft max(256, 2^nextpow2(wlen)); win hamming(wlen, periodic); % 计算时间轴和帧数 N length(x); nframes floor((N - wlen) / hop) 1; S zeros(nfft/2, nframes); t zeros(nframes, 1); for i 1:nframes idx_start (i-1)*hop 1; seg x(idx_start : idx_startwlen-1) .* win; X fft(seg, nfft); S(:, i) abs(X(1:nfft/2)); t(i) (idx_start wlen/2) / Fs; end f (0:nfft/2-1) * Fs / nfft; end时序循环里做了三件容易被忽略的事第一hop wlen * (1-overlap)定义了帧移也就是相邻两帧起点的时间差第二nfft通过nextpow2取到不小于窗长的 2 的幂方便 FFT 计算但引入了补零第三时间轴t(i)取窗的中心位置对应那一帧频谱的“时刻”而不是窗口起点这能避免时频谱图在视觉上整体右移半个窗长。用这个函数观察一个 1 kHz 到 5 kHz 的线性扫频信号时频谱图上能看到一条从左下到右上的连续亮线。如果亮线明显断裂说明overlap太小或hop大于扫频速率可接受的粒度如果亮线很粗但位置模糊说明wlen太长时间分辨率不足。这是调 STFT 参数最直接的物理反馈。参数作用调优方向wlen频率分辨率与时间分辨率的折中增大则频率变清晰、时间变模糊overlap控制帧间连续性越大亮线越连续开销越高nfftFFT 长度取2^nextpow2(wlen)即可增大只做补零hop帧移期望每帧时间步长直接影响时间轴精度实践里我一般把wlen设为目标信号最低频率对应周期的 5 到 10 倍。比如最低关注频率是 50 Hz周期 20 ms那么wlen取 100 到 200 ms对应采样率 100 kHz 下就是 10000 到 20000 点。比这个短会看不到低频细节比这个长则高频事件被平均掉。4. env_phas.m 解包络与相位从解析信号到瞬时频率的调制解调验证4.1 为什么解调不用直接取 abs时频分析给出了“信号里有什么”但要还原发送的信息还得把调制过程逆回来。包络检波是 AM 最直观的解调方式取abs(x)再低通滤波信息就回来了。这个做法在纯单音信号上没问题一旦信号混入噪声、多径或经过频带限制直接取绝对值会引入大量高频纹波且对直流偏置敏感。更稳的做法是先构造解析信号再从中分离包络和瞬时相位。Hilbert 变换是构造解析信号的标准工具z(t) x(t) j*hilbert(x(t))。解析信号的实部是原信号虚部是原信号的 90 度相移版本称为正交分量。有了z(t)瞬时幅度就是abs(z)瞬时相位是angle(z)瞬时频率是相位的导数除以2*pi。这个过程在 MATLAB 里只需要一行hilbert但背后的“单边谱”假设值得留意解析信号只保留正频率分量所以 Hilbert 变换天然丢弃负频率如果原始信号在Fs/2附近有能量构造结果会出现混叠。4.2 env_phas.m 的核心逻辑env_phas.m名字里的 env 和 phas 分别对应包络envelope和相位phase实现的就是上面这套流程。下面给出一份可直接使用的等价实现并加上瞬时频率输出。function [env, inst_phase, inst_freq] env_phas(s, Fs) % [env, inst_phase, inst_freq] env_phas(s, Fs) % s : 输入信号实信号或复基带信号 % Fs : 采样率 % env : 瞬时包络长度与 s 相同 % inst_phase : 瞬时相位rad已做 unwrap % inst_freq : 瞬时频率Hz用中心差分计算 % 解析信号对实信号直接 hilbert对复信号只需取原值 if isreal(s) z hilbert(s); else z s; end env abs(z); inst_phase unwrap(angle(z)); % 中心差分求瞬时频率边界用单侧差分 inst_freq zeros(size(inst_phase)); inst_freq(2:end-1) (inst_phase(3:end) - inst_phase(1:end-2)) * Fs / (4*pi); inst_freq(1) (inst_phase(2) - inst_phase(1)) * Fs / (2*pi); inst_freq(end) (inst_phase(end) - inst_phase(end-1)) * Fs / (2*pi); % 去除因差分造成的边缘跳变 inst_freq(abs(inst_freq) Fs/2) NaN; end代码里有三个地方是工程上容易踩坑的。第一unwrap把相位从[-pi, pi]展开成连续曲线但它是逐点判断跳变的如果信噪比过低导致相邻两点相位差超过piunwrap会把真实跳变误判为翻折产生虚假的瞬时频率尖峰。第二瞬时频率用中心差分计算频率分辨率受限于采样率与相位噪声的比值数据里如果还有直流偏置angle(z)会被直流分量抬高导致inst_freq整体偏移。第三对复基带信号直接解析而不再过 Hilbert因为s本身已经是单边带表示。4.3 调制解调验证实验端到端跑通把三个文件串起来做一次完整的调制-解调验证用 1 kHz 正弦作为信息分别调 AM 和 FM再通过env_phas解调对比解调输出与原信息。% 端到端验证脚本 Fs 200e3; fm 1e3; N 20000; t (0:N-1)./Fs; info 0.8 * sin(2*pi*fm*t); % 幅度 0.8避免过调幅 % AM 调制与解调 s_am loweq(am, info, Fs, fm, 20e3, 0.5); rx_am real(s_am); % 模拟接收端取实部 [env_am, ~, ~] env_phas(rx_am, Fs); demod_am env_am - mean(env_am); % 去除直流分量 % FM 调制与解调 s_fm loweq(fm, info, Fs, fm, 20e3, 1.0); [~, ~, inst_freq] env_phas(real(s_fm), Fs); demod_fm inst_freq - mean(inst_freq); % 指标对比 mse_am mean((demod_am - info).^2); mse_fm mean((demod_fm - info).^2); fprintf(AM 解调 MSE: %.2e\n, mse_am); fprintf(FM 解调 MSE: %.2e\n, mse_fm);AM 解调输出与原始info之间只差一个直流偏置和缩放因子因此减均值再归一化后能高度重合。FM 解调的关键是把“相位导数”这个中间量换算回幅度域理想情况下inst_freq的交流部分就是信息信号本身但要先减均值以去掉固定频偏。实际运行时如果mse_am远大于mse_fm优先怀疑 AM 调制时mf设置过低导致包络变化太微弱如果mse_fm异常检查unwrap后是否还有残留跳变。这套包络-相位提取方法也可以直接用在振动信号分析上。比如滚动轴承故障时产生的冲击会调制在高频固有振动上对包络信号再做一次 FFT 得到“包络谱”频谱上的故障特征频率就是轴承故障的直接证据。这也是供水管网噪声记录仪做频谱分析、频带划分时常用的前端思路先提取包络再对包络做窄带分析。5. 参数自查与诊断脚本让这套时频分析工具在新数据上跑通5.1 拿到陌生数据时先查这四个参数把main.m里的仿真信号换成现场采集的真实数据第一轮跑出来的结果通常比较难看。过去踩过的坑集中在这四个地方症状可能原因修法频谱横轴数值离谱忘记把 FFT 点数换算成f (0:N-1)*Fs/N核对采样率缺失就先归一化看形状谱线全部挤在零频附近数据含强直流或趋势项先detrend再去均值必要时高通滤波时频谱图有斜向条纹分段帧移与分析带宽不匹配检查hop与窗长关系增大重叠率到 0.75 以上解调输出噪声起伏剧烈采样率低于信号最高频率的 2 倍确认前端抗混叠滤波是否开启5.2 一个可嵌入 main.m 的自动检查块与其每次手动猜不如在频谱分析之前加一段自动诊断降低排查成本。% --- 自动诊断块嵌入 main.m 的数据加载之后 --- assert(~isempty(x), 输入信号为空); N length(x); Fs_est max(2 * abs(diff(angle(hilbert(x)))) * Fs / (2*pi)); % 粗估最大瞬时频率 fprintf(信号长度: %d\n, N); fprintf(实测最高瞬时频率: %.2f Hz\n, Fs_est); % 直流检查直流占比超过 30% 时提示去趋势 dc_ratio abs(mean(x)) / (std(x) eps); if dc_ratio 0.3 warning(直流分量偏大建议先 detrend 再分析); end % 峭度检查峭度过高说明存在明显脉冲优先看包络谱 k kurtosis(x); if k 5 fprintf(峭度 %.2f 5建议先做包络分析再观察频谱\n, k); end这段代码里的assert防止空数据进入后续流程Fs_est用 Hilbert 相位差粗估信号实际占用带宽一旦发现接近Fs/2就说明混叠风险很高。峭度是振动分析里常用指标大于 5 往往意味着信号里有冲击成分这时直接看功率谱会把冲击的宽频能量和连续谱混在一起先解包络再做窄带频谱更有说服力。把这个诊断块放在main.m的头部每次跑新数据都能先得到一组“健康指标”再决定是用全局 FFT、时频谱图还是包络谱。这套工作流并不限于 MATLABPython 里scipy.signal.hilbert和scipy.signal.spectrogram可以做完全相同的操作只是参数语义略有差别。分析逻辑不变换语言只是换函数名而已。本文还有配套的精品资源点击获取 基于MATLAB的等效低通调制与解调：从FFT频谱分析到瞬时频率提取
简介频谱分析与时频分析是通信与信号处理领域的基础工具本资源面向电子工程、通信技术方向的学习者与研发人员聚焦调制解调信号观测、时频谱图绘制与频谱特征提取等常见需求适合在课程实验或项目仿真中参考使用。资源共3个文件均为MATLAB脚本.m压缩包仅2KB体积小巧但代码结构集中覆盖信号产生、低通等效变换与包络相位计算等关键环节可直接用于快速验证频谱和时频分析方法。目前已有717人学习下载。通过研读这些源码读者可以掌握如何构造分析信号、调用时频变换函数、绘制时频谱图以及解读频谱分布同时对照描述中的调制解调、STFT、功率谱密度等概念加深对非平稳信号处理流程的理解并为后续通信系统设计与故障检测分析提供可复用的脚本基础。1. 频谱分析不是只有一个 FFT从三个 m 文件看一条完整信号链拿到一套只有main.m、loweq.m、env_phas.m三个脚本的 MATLAB 工程时我一开始也以为这只是个画频谱图的演示程序。跑完才发现这三个文件合起来实际上是一台简化版的软件无线电分析仪loweq.m负责把基带信息“搬”到等效低通信号上main.m完成波形生成和频谱测量env_phas.m再从解析信号里把包络、瞬时相位和瞬时频率一一抽出来。换句话说它覆盖了调制、传输、解调、频域测量和时频图绘制的完整闭环。对做通信物理层、信号处理和嵌入式 FFT 分析的人来说这套代码的价值不在于某一个算法多高级而在于它把“频谱分析”从单个 FFT 函数扩展成了一条可验证的链路。想搞清楚非平稳信号怎么观察、调制参数怎么定、解调误差从哪里来顺着这三个脚本走一遍比读十页公式更直观。2. 调制端建模用 loweq.m 把等效低通和真实载频解耦2.1 为什么要用“等效低通”而不是直接生成载波很多人在 MATLAB 里仿真调制信号时习惯直接写y cos(2*pi*fc*t phi)把载频fc取到几百 MHz采样率跟着拉到 GHz 量级。这样做的直接后果是一个 1 秒的仿真就要存几十 GB 数据FFT 点数大得离谱而真正有用的信息带宽往往只有几 kHz。纯属浪费算力。通信系统仿真里更常见的做法是**等效低通complex baseband / equivalent lowpass**模型。原理很简单实带通信号x(t) Re{s(t) * exp(j*2*pi*fc*t)}其中s(t)是复包络。只要把复包络算出来载频只是一个乘法因子甚至可以在分析阶段完全忽略。这样一来采样率只需要覆盖信号带宽而不是覆盖载频仿真数据量直接下降两三个数量级。这也是 GMSK 这类恒包络调制能高效仿真的原因。GMSK 调制解调原理看着复杂本质就是在 IQ 平面上做相位累积再让相位路径经过高斯滤波变得平滑从而压缩频谱旁瓣。用等效低通模型来写核心代码不会超过二十行。2.2 loweq.m 的典型实现与参数含义下面给出一份和资源中loweq.m功能等价的实现。它把 AM 和 FM 两种调制方式统一在同一个函数里输出复基带信号。function [s, t] loweq(mode, info, Fs, fm, fc, mf) % [s, t] loweq(mode, info, Fs, fm, fc, mf) % mode : am 幅度调制 / fm 频率调制 % info : 基带信息序列单声道列向量 % Fs : 采样率必须大于 2*fm这里 fm 是信息最高频率 % fm : 信息信号最高频率用于定时间轴和成型 % fc : 等效载频只影响相位旋转不参与采样率设定 % mf : 调制度AM 用调幅深度FM 用频偏系数 Ts 1/Fs; N length(info); t (0:N-1). * Ts; switch lower(mode) case am % 直流分量 1 保证包络不会过零mf 取 0~1 之间 s (1 mf * info) .* exp(1j*2*pi*fc*t); case fm % 相位是信息的积分mf 等效为频偏/信息频率 phase 2*pi*mf*fm * cumsum(info) * Ts; s exp(1j * (phase 2*pi*fc*t)); otherwise error(loweq: unknown mode %s, mode); end end这段代码的关键点有三个。第一AM 模式里exp(1j*2*pi*fc*t)是未被显式调用的复载波实际信号是real(s)相当于把载波搬到了fc处。第二FM 模式用cumsum(info) * Ts做离散积分这正是角度调制里“相位是调制信号的积分”这一结论的直接翻译。第三Fs只由info的带宽决定和fc没有关系这是等效低通建模最反直觉也最容易被忽略的地方。参数调优时mf的选择直接影响解调难度。AM 的mf超过 1 会导致过调幅包络出现截断包络检波后谐波失真明显FM 的mf决定最大频偏频偏过大则瞬时频率可能超出分析带宽时频谱图上会出现折叠。我在调试时一般先把mf设为 0.5 做 AM、1.0 做 FM在确认整条链路无误后再加大。2.3 验证调制器输出是否正确写完调制函数后先做两个快速检查一是看频谱形状二是看瞬时频率是否符合预期。Fs 100e3; fm 1e3; N 8192; t (0:N-1)./Fs; info sin(2*pi*fm*t); % 1 kHz 单音信息 s_am loweq(am, info, Fs, fm, 10e3, 0.5); s_fm loweq(fm, info, Fs, fm, 10e3, 1.0); % 频谱检查AM 应出现载波 ± 1kHz 两根边带FM 的谱线间隔等于 fm f_axis (0:N-1)/N*Fs - Fs/2; figure; subplot(2,1,1); plot(f_axis/1e3, fftshift(abs(fft(s_am)))); title(AM spectrum); xlim([-5 20]); subplot(2,1,2); plot(f_axis/1e3, fftshift(abs(fft(s_fm)))); title(FM spectrum); xlim([-5 20]);看 AM 频谱时理想情况下载波谱线在 10 kHz 处两侧 9 kHz 和 11 kHz 各有一根对称边带幅度是载波的一半。FM 谱相对复杂但谱线间隔严格等于fm且高次边带随阶数衰减。如果谱线间隔不对先检查Fs和fm的量纲是否匹配如果出现非预期的大块噪声底则说明cumsum的累积误差在起作用可以考虑改用filter(1, [1 -1], info)做抗饱和积分。3. main.m 里的频谱与极坐标图FFT 参数、频率轴标定和时频谱图的边界效应3.1 main.m 的主流程拆解main.m在整个工程里承担的角色是“总装车间”一般包含五个步骤设置仿真参数、调用loweq生成信号、叠加噪声、进行频谱分析、绘制结果。其中频谱分析部分写得是否严谨直接影响后续对调制方式的判断。很多人在这一步把fft()的结果直接plot然后对着横轴发呆因为横轴是“点数”而不是“频率”。正确的频率轴标定公式是f (0:N-1) * Fs / N对应到单边谱则是f (0:N/2-1) * Fs / N。采样率是频率轴唯一的“定标尺”不知道采样率FFT 只能给出归一化频率0 到 1 对应 0 到 Fs这也是“不知道采样率怎么求频率频谱”这个问题最常见的场景——工程数据一旦丢了采样率频域分析基本只能看形状无法读数。3.2 频率分辨率、补零和窗函数频谱分析的核心指标是频率分辨率Δf Fs / N。注意补零能改变N从而细化频点间隔但它不能提高真正的分辨率因为有效信息长度没有变。想分辨两根间隔 100 Hz 的谱线FFT 点数至少要大于Fs / 100补零只是让谱线看起来更圆润。Fs 100e3; N 10000; t (0:N-1)./Fs; x sin(2*pi*1000*t) 0.5*sin(2*pi*1020*t); % 间隔 20 Hz 的两根谱线 % 直接 FFT矩形窗 X1 fft(x); % 加 hamming 窗后再 FFT w hamming(N); X2 fft(x .* w); % 频率轴单边谱 f_axis (0:N/2-1) * Fs / N; figure; subplot(2,1,1); plot(f_axis, 20*log10(abs(X1(1:N/2))/N)); title(Rectangular window); xlim([900 1100]); subplot(2,1,2); plot(f_axis, 20*log10(abs(X2(1:N/2))/N)); title(Hamming window); xlim([900 1100]);这段代码直接体现了窗函数的作用矩形窗下 1000 Hz 和 1020 Hz 两根谱线大概率粘连成一个峰因为 20 Hz 的间隔小于Δf Fs/N 10 Hz。加 Hamming 窗后主瓣变宽旁瓣抑制增强但 20 Hz 间隔在 10 Hz 分辨率下依然不一定能分开只是不会再被旁瓣淹没。参数选择上hamming(N)的等效噪声带宽比矩形窗宽约 1.36 倍所以幅度读数要按窗的相干增益修正否则测出来的功率偏低。3.3 从全局 FFT 到时频谱图分段 FFT 的正确姿势全局 FFT 适合平稳信号换成跳频、扫频或调频信号就失效了。观察这类非平稳信号标准做法是短时傅里叶变换STFT把信号切成一小段一小段每段做 FFT再把结果按时间顺序拼成二维矩阵。这就是“分段频谱图”的实质也就是时频分析中时频谱图的来源。分段 FFT 有两个参数需要手工调窗口长度wlen和重叠率overlap。窗口长度决定时间分辨率和频率分辨率的折中窗越长频率分辨率越高但时间上越模糊重叠率越高时频谱图越平滑但计算量越大。MATLAB 自带的spectrogram函数可以直接用但用矩阵切片手写一段会更清楚每一步在做什么。function [S, f, t] my_stft(x, Fs, wlen, overlap) % 手动实现短时傅里叶变换 % wlen : 窗长度点数 % overlap: 重叠率0~0.99 hop round(wlen * (1 - overlap)); nfft max(256, 2^nextpow2(wlen)); win hamming(wlen, periodic); % 计算时间轴和帧数 N length(x); nframes floor((N - wlen) / hop) 1; S zeros(nfft/2, nframes); t zeros(nframes, 1); for i 1:nframes idx_start (i-1)*hop 1; seg x(idx_start : idx_startwlen-1) .* win; X fft(seg, nfft); S(:, i) abs(X(1:nfft/2)); t(i) (idx_start wlen/2) / Fs; end f (0:nfft/2-1) * Fs / nfft; end时序循环里做了三件容易被忽略的事第一hop wlen * (1-overlap)定义了帧移也就是相邻两帧起点的时间差第二nfft通过nextpow2取到不小于窗长的 2 的幂方便 FFT 计算但引入了补零第三时间轴t(i)取窗的中心位置对应那一帧频谱的“时刻”而不是窗口起点这能避免时频谱图在视觉上整体右移半个窗长。用这个函数观察一个 1 kHz 到 5 kHz 的线性扫频信号时频谱图上能看到一条从左下到右上的连续亮线。如果亮线明显断裂说明overlap太小或hop大于扫频速率可接受的粒度如果亮线很粗但位置模糊说明wlen太长时间分辨率不足。这是调 STFT 参数最直接的物理反馈。参数作用调优方向wlen频率分辨率与时间分辨率的折中增大则频率变清晰、时间变模糊overlap控制帧间连续性越大亮线越连续开销越高nfftFFT 长度取2^nextpow2(wlen)即可增大只做补零hop帧移期望每帧时间步长直接影响时间轴精度实践里我一般把wlen设为目标信号最低频率对应周期的 5 到 10 倍。比如最低关注频率是 50 Hz周期 20 ms那么wlen取 100 到 200 ms对应采样率 100 kHz 下就是 10000 到 20000 点。比这个短会看不到低频细节比这个长则高频事件被平均掉。4. env_phas.m 解包络与相位从解析信号到瞬时频率的调制解调验证4.1 为什么解调不用直接取 abs时频分析给出了“信号里有什么”但要还原发送的信息还得把调制过程逆回来。包络检波是 AM 最直观的解调方式取abs(x)再低通滤波信息就回来了。这个做法在纯单音信号上没问题一旦信号混入噪声、多径或经过频带限制直接取绝对值会引入大量高频纹波且对直流偏置敏感。更稳的做法是先构造解析信号再从中分离包络和瞬时相位。Hilbert 变换是构造解析信号的标准工具z(t) x(t) j*hilbert(x(t))。解析信号的实部是原信号虚部是原信号的 90 度相移版本称为正交分量。有了z(t)瞬时幅度就是abs(z)瞬时相位是angle(z)瞬时频率是相位的导数除以2*pi。这个过程在 MATLAB 里只需要一行hilbert但背后的“单边谱”假设值得留意解析信号只保留正频率分量所以 Hilbert 变换天然丢弃负频率如果原始信号在Fs/2附近有能量构造结果会出现混叠。4.2 env_phas.m 的核心逻辑env_phas.m名字里的 env 和 phas 分别对应包络envelope和相位phase实现的就是上面这套流程。下面给出一份可直接使用的等价实现并加上瞬时频率输出。function [env, inst_phase, inst_freq] env_phas(s, Fs) % [env, inst_phase, inst_freq] env_phas(s, Fs) % s : 输入信号实信号或复基带信号 % Fs : 采样率 % env : 瞬时包络长度与 s 相同 % inst_phase : 瞬时相位rad已做 unwrap % inst_freq : 瞬时频率Hz用中心差分计算 % 解析信号对实信号直接 hilbert对复信号只需取原值 if isreal(s) z hilbert(s); else z s; end env abs(z); inst_phase unwrap(angle(z)); % 中心差分求瞬时频率边界用单侧差分 inst_freq zeros(size(inst_phase)); inst_freq(2:end-1) (inst_phase(3:end) - inst_phase(1:end-2)) * Fs / (4*pi); inst_freq(1) (inst_phase(2) - inst_phase(1)) * Fs / (2*pi); inst_freq(end) (inst_phase(end) - inst_phase(end-1)) * Fs / (2*pi); % 去除因差分造成的边缘跳变 inst_freq(abs(inst_freq) Fs/2) NaN; end代码里有三个地方是工程上容易踩坑的。第一unwrap把相位从[-pi, pi]展开成连续曲线但它是逐点判断跳变的如果信噪比过低导致相邻两点相位差超过piunwrap会把真实跳变误判为翻折产生虚假的瞬时频率尖峰。第二瞬时频率用中心差分计算频率分辨率受限于采样率与相位噪声的比值数据里如果还有直流偏置angle(z)会被直流分量抬高导致inst_freq整体偏移。第三对复基带信号直接解析而不再过 Hilbert因为s本身已经是单边带表示。4.3 调制解调验证实验端到端跑通把三个文件串起来做一次完整的调制-解调验证用 1 kHz 正弦作为信息分别调 AM 和 FM再通过env_phas解调对比解调输出与原信息。% 端到端验证脚本 Fs 200e3; fm 1e3; N 20000; t (0:N-1)./Fs; info 0.8 * sin(2*pi*fm*t); % 幅度 0.8避免过调幅 % AM 调制与解调 s_am loweq(am, info, Fs, fm, 20e3, 0.5); rx_am real(s_am); % 模拟接收端取实部 [env_am, ~, ~] env_phas(rx_am, Fs); demod_am env_am - mean(env_am); % 去除直流分量 % FM 调制与解调 s_fm loweq(fm, info, Fs, fm, 20e3, 1.0); [~, ~, inst_freq] env_phas(real(s_fm), Fs); demod_fm inst_freq - mean(inst_freq); % 指标对比 mse_am mean((demod_am - info).^2); mse_fm mean((demod_fm - info).^2); fprintf(AM 解调 MSE: %.2e\n, mse_am); fprintf(FM 解调 MSE: %.2e\n, mse_fm);AM 解调输出与原始info之间只差一个直流偏置和缩放因子因此减均值再归一化后能高度重合。FM 解调的关键是把“相位导数”这个中间量换算回幅度域理想情况下inst_freq的交流部分就是信息信号本身但要先减均值以去掉固定频偏。实际运行时如果mse_am远大于mse_fm优先怀疑 AM 调制时mf设置过低导致包络变化太微弱如果mse_fm异常检查unwrap后是否还有残留跳变。这套包络-相位提取方法也可以直接用在振动信号分析上。比如滚动轴承故障时产生的冲击会调制在高频固有振动上对包络信号再做一次 FFT 得到“包络谱”频谱上的故障特征频率就是轴承故障的直接证据。这也是供水管网噪声记录仪做频谱分析、频带划分时常用的前端思路先提取包络再对包络做窄带分析。5. 参数自查与诊断脚本让这套时频分析工具在新数据上跑通5.1 拿到陌生数据时先查这四个参数把main.m里的仿真信号换成现场采集的真实数据第一轮跑出来的结果通常比较难看。过去踩过的坑集中在这四个地方症状可能原因修法频谱横轴数值离谱忘记把 FFT 点数换算成f (0:N-1)*Fs/N核对采样率缺失就先归一化看形状谱线全部挤在零频附近数据含强直流或趋势项先detrend再去均值必要时高通滤波时频谱图有斜向条纹分段帧移与分析带宽不匹配检查hop与窗长关系增大重叠率到 0.75 以上解调输出噪声起伏剧烈采样率低于信号最高频率的 2 倍确认前端抗混叠滤波是否开启5.2 一个可嵌入 main.m 的自动检查块与其每次手动猜不如在频谱分析之前加一段自动诊断降低排查成本。% --- 自动诊断块嵌入 main.m 的数据加载之后 --- assert(~isempty(x), 输入信号为空); N length(x); Fs_est max(2 * abs(diff(angle(hilbert(x)))) * Fs / (2*pi)); % 粗估最大瞬时频率 fprintf(信号长度: %d\n, N); fprintf(实测最高瞬时频率: %.2f Hz\n, Fs_est); % 直流检查直流占比超过 30% 时提示去趋势 dc_ratio abs(mean(x)) / (std(x) eps); if dc_ratio 0.3 warning(直流分量偏大建议先 detrend 再分析); end % 峭度检查峭度过高说明存在明显脉冲优先看包络谱 k kurtosis(x); if k 5 fprintf(峭度 %.2f 5建议先做包络分析再观察频谱\n, k); end这段代码里的assert防止空数据进入后续流程Fs_est用 Hilbert 相位差粗估信号实际占用带宽一旦发现接近Fs/2就说明混叠风险很高。峭度是振动分析里常用指标大于 5 往往意味着信号里有冲击成分这时直接看功率谱会把冲击的宽频能量和连续谱混在一起先解包络再做窄带频谱更有说服力。把这个诊断块放在main.m的头部每次跑新数据都能先得到一组“健康指标”再决定是用全局 FFT、时频谱图还是包络谱。这套工作流并不限于 MATLABPython 里scipy.signal.hilbert和scipy.signal.spectrogram可以做完全相同的操作只是参数语义略有差别。分析逻辑不变换语言只是换函数名而已。本文还有配套的精品资源点击获取