希尔伯特变换解调与DEMON谱：MATLAB实现水声轴频特征提取

简介面向水声目标识别与辐射噪声分析研究者的 MATLAB 源码小包聚焦希尔伯特变换解调在目标信号特征提取中的应用通过仿真信号和实际船舶辐射噪声信号的对比帮助理解不同解调方法在真实水声环境下的性能差异。资源共 2 个文件约 2KB包含一个 MATLAB 脚本和一个文本数据/记录文件脚本用于实现三种解调方法的对比文本文件对应仿真与实测信号处理结果便于复现与验证。压缩包体积小、内容聚焦适合信号处理或水声工程方向学生快速上手。已有 682 人学习/下载。从价值来看这份代码可作为希尔伯特变换解调与 DEMON 谱分析的入门参考既能看到算法实现的代码框架也能对比不同方法的解调谱效果进一步理解辐射噪声如何影响目标识别为后续使用 SVM 或神经网络等识别模型提供有效的特征输入。1. 辐射噪声里的“调制指纹”为什么值得用希尔伯特变换解调海试记录里目标船的螺旋桨轴频会周期性地调制宽带辐射噪声。直接对水听器信号做FFT轴频线谱被埋在连续谱里看不出区别但先把信号包络取出来再做一次FFT轴频和叶频就会像“指纹”一样出现在低频调制谱上。这就是DEMON谱的基本思想也是很多水声目标识别系统里成本最低、稳定性最高的特征来源。这次拆解的是两个MATLAB脚本P340217.m.txt是带数据读取动作的主流程sanzhongduibiP34.m则对同一条信号做了三种解调方法的横向比较。我会把希尔伯特变换解调的原理、代码实现和参数调整顺序过一遍适合刚跑通水声数据、想从辐射噪声里稳定提取轴频特征的人。2. 希尔伯特变换解调与DEMON谱的数学基础2.1 从实信号到解析信号包络提取为什么非它不可水听器接收到的辐射噪声x(t)可以看作一个宽带载波被螺旋桨轴频f_shaft周期调制。在工程上这个调制过程通常写成x(t) A0 * [1 m * cos(2*pi*f_shaft*t)] * n_band(t)其中n_band(t)是空化噪声的窄带分量m是调制深度。对目标识别来说f_shaft就是我们要找到的“指纹”。直接对x(t)做FFT调制边带被淹没在宽带噪声中但把宽带信号包络取出来再做FFTf_shaft就出现在低频段且背景平坦、峰位清晰。直接取绝对值也能得到包络但那不是解析意义上的包络。|sin|会产生强谐波对带通滤波后的信号做平方又会引入直流偏置并放大倍频。希尔伯特变换的优势在于构造解析信号z(t) x(t) j * H{x(t)}其中虚部是实部的正交版本。包络A(t) |z(t)|在数学上是唯一的最小相位包络不会因为载波频率选取不当而额外产生调制分量。在MATLAB里hilbert(x)返回的是解析信号不是纯希尔伯特变换本身。它先对x做FFT把负频率置零再逆FFT所以abs(hilbert(x))就是包络。需要注意hilbert默认按整段数据计算边界会有Gibbs振荡。我一般会先把数据分段再做包络或者用envelope函数替代避免首尾异常尖峰污染后续谱平均。2.2 调制谱参数频段、窗长、帧数与频率分辨率DEMON分析的参数不像普通频谱分析那样只选FFT点数。第一个要选的是带通频段[fL, fH]。辐射噪声的宽带连续谱在低频段能量高但螺旋桨调制通常要在几百Hz到几千Hz的频段内观察因为那里空化噪声的调制深度更明显。带通选得过窄包络里会保留载波泄漏选得过宽其他机械噪声源会混进来。第二个关键参数是包络后做FFT的帧长T_frame。调制频率一般在几Hz到几十Hz频率分辨率由df 1 / T_frame决定而不是由原始采样率决定。T_frame取1秒以上才能分辨1Hz以下的差别。常见参数如下参数典型值作用与影响分析频段 [fL, fH]500–3000 Hz决定进入解调的频带需避开强线谱帧长 T_frame0.5–2 s调制谱分辨率 df1/T_frame帧重叠率50%增加平均次数降低包络谱方差包络降采样倍数8–20降低FFT点数缩短计算时间谱平均次数20–50稳定轴频峰抑制随机毛刺包络信号的有效采样率等于降采样后的fs_env。理论上可观测的最大调制频率是fs_env/2降采样倍数不能过大。比如原始采样率50kHz降16倍后fs_env 3125 Hz足够容纳轴频几十倍频的观测范围。如果轴频只有几Hz把帧长放到2sdf 0.5 Hz能清晰分离轴频与电源50Hz干扰。包络降采样用resample会默认做抗混叠滤波这一点比直接抽取安全。实测中直接写env(1:16:end)会出现伪峰因为包络中仍可能有高频残差被折叠到低频段。所以代码里应优先使用resample(env, 1, 16)而不是手动抽取。2.3 能量解调、绝对值解调与希尔伯特解调的统一框架三种解调方法可以用同一个流程描述先对x做带通滤波得到x_bp再做非线性变换最后去直流。工程上常用的三种方式是能量解调y x_bp.^2绝对值解调y abs(x_bp)以及希尔伯特解调y abs(hilbert(x_bp))。从频域看平方等效于信号与自身卷积会把原带通信号的频谱搬移到两倍频和零频绝对值操作会产生基频、倍频和直流叠加希尔伯特解调则先移频再低通理论上不产生额外倍频。这解释了为什么希尔伯特解调得到的调制谱更干净它保留了窄带信号的包络信息只输出真正的调制频率分量。但希尔伯特解调对带通滤波的质量更敏感。如果频带内混入两个不同调制状态的强线谱解析包络会把它们“拍频”成一条虚假谱线。所以进入hilbert之前带通滤波器的阻带衰减建议做到40dB以上。用designfilt设计巴特沃斯时阶数至少6阶必要时8阶而不是随手用MATLAB默认的低阶滤波器。3. P340217.m单条辐射噪声信号的DEMON处理流程3.1 数据加载与预处理P340217.m.txt是文本格式的MATLAB脚本数据一般放在同目录的txt或dat文件里。从命名规律看P340217应该是某个航次记录编号。读取这类数据时常见做法是先看文件前几行确定列布局再决定用load还是importdata。我一般会先把数据转成列向量并去除直流偏置因为后续的希尔伯特包络需要一个零均值输入。raw load(P340217_data.txt); % 每行一个采样点也可能是多列 if size(raw, 2) 1 x raw(:, 1); % 取水听器通道 else x raw; end x x - mean(x); % 去直流否则包络谱0Hz处会出现大峰 fs 50000; % 实际脚本中需要按记录文件的采样率修改这里按文本读取时如果数据量超过500MB建议改用memmapfile避免一次性占用过多内存。去直流不能省因为希尔伯特包络本身会产生直流分量原始信号若有直流偏置这个分量会被额外放大最终在调制谱低频端形成很强的斜坡。带通滤波器用designfilt设计通常选择巴特沃斯或椭圆滤波器。阶数过高会带来时域振铃过低则带外衰减不足。下面的代码使用6阶巴特沃斯带通阻带衰减约30dB足够做初步解调若后续发现带外干扰明显再把阶数提高到8。fL 500; fH 3000; bpf designfilt(bandpassiir, FilterOrder, 6, ... HalfPowerFrequency1, fL, ... HalfPowerFrequency2, fH, ... SampleRate, fs); x_bp filter(bpf, x);filter沿时间方向滤波边界会有瞬态响应。实际处理时前1000点通常需要丢弃或者在滤波前对数据前后各补一段镜像滤波后截掉避免包络首尾失真。3.2 希尔伯特包络解调的实现核心代码只有三行细节都在包络的后续处理上。先取解析信号幅度再去直流再降采样。降采样放在去直流之后可以让resample内部的低通滤波器不被直流偏置影响。env_raw abs(hilbert(x_bp)); % 解析信号幅度 env env_raw - mean(env_raw); % 去直流保留调制分量 env resample(env, 1, 16); % 降采样到 fs/16 fs_env fs / 16;hilbert内部对整段数据做一次FFT数据很长时计算开销很大。如果记录超过30秒建议把x_bp切成互不重叠的块逐块计算包络拼回后再降采样。否则一次hilbert调用可能会卡住MATLAB主线程很多秒。降采样倍数的选择要看原始采样率。50kHz采样率下降16倍得到fs_env 3125Hz包络谱最高可以观察到1562Hz的调制频率对螺旋桨轴频几十次谐波都足够。如果轴频本身只有几Hz可以把降采样倍数提高到32进一步降低后续帧处理的维度。3.3 谱平均与峰值提取包络谱的标准做法是分帧后加窗再平均。直接对整段包络做FFT数据非平稳成分会把谱峰拉宽用pwelch平均则损失每一帧内调制相位信息。所以这里手动切帧帧长2s重叠1s逐帧做幅度谱后再平均保留线谱幅度用于峰值比较。frame_len round(fs_env * 2); % 2秒一帧 hop round(frame_len / 2); win hann(frame_len, periodic); n_frames floor((length(env) - frame_len) / hop) 1; spec_sum zeros(frame_len, 1); for k 1:n_frames idx (k-1)*hop (1:frame_len).; seg env(idx) .* win; seg seg - mean(seg); S abs(fft(seg)); % 幅度谱 spec_sum spec_sum S; end spec_avg spec_sum / n_frames; f_mod (0:frame_len-1). * fs_env / frame_len; % 只观察0-100Hz调制频段 idx_mod find(f_mod 0 f_mod 100); [pks, locs] findpeaks(spec_avg(idx_mod), f_mod(idx_mod), ... MinPeakHeight, max(spec_avg(idx_mod))*0.3, ... MinPeakDistance, 0.5);hann(frame_len, periodic)是DEMON分析里常用的窗周期汉宁窗旁瓣更低能减少谱泄漏对相邻频点的影响。逐帧减mean(seg)是必要的因为每帧包络均值不同整体减去全局均值会留下帧间直流起伏。MinPeakDistance设为0.5Hz防止同一个谱峰因窗函数旁瓣被拆成两个峰。MinPeakHeight设为最大峰高的30%这个阈值比较激进只适合初步筛选。后续应当结合谐波关系确认叶频等于轴频乘以叶片数一般在包络谱中能看到轴频、二倍轴频和叶片数倍频才算真正锁定了目标。3.4 常见错误与工程建议现象可能原因解决方式包络谱0Hz处有巨峰未去直流或帧内未减均值滤波、解调、加窗前后各减一次均值轴频峰宽远超1Hz帧长过短或数据有明显频漂加长帧长检查整段转速是否稳定高频处出现等间距伪峰降采样前混叠使用resample而非抽取相同参数下不同数据峰位跳变带通频段选错拍频产生虚假峰扫描带通范围交叉验证我处理P340217这类数据时会把3.1到3.3写成一个函数输入[fL, fH, T_frame, decim]输出包络谱和峰值表。这样在后续对比不同频段时不需要反复复制脚本代码也能避免改了参数忘记恢复的尴尬。4. sanzhongduibiP34.m三种解调方法的对比实验4.1 对比方案设计sanzhongduibiP34.m这个脚本名已经说明了任务对P34号记录做三种解调方法对比。对比的前提是必须保证三个分支使用完全相同的带通滤波器、帧长、重叠率、窗函数和降采样倍数唯一不同的只有包络提取的公式。否则任何谱图差异都可能是参数不一致导致的不能归因于方法本身。方法包络公式优点风险能量解调y x_bp.^2实现简单噪声功率压缩轴频倍频明显直流偏置强绝对值解调y abs(x_bp)不涉及复计算速度快频谱有奇次谐波基频幅度不稳希尔伯特解调y abs(hilbert(x_bp))包络干净谱线锐利对带外泄漏敏感边界有振荡在MATLAB早期版本中hilbert对长序列的计算开销很大很多老代码选择平方解调。但现在的机器上希尔伯特解调不再有性能瓶颈除非数据超过几百MB。推荐把三种方法都实现一遍用同一组帧参数输出三张包络谱再叠加对比。4.2 三种方法的MATLAB实现与结果判读在sanzhongduibiP34.m里三类解调核心代码通常只有一行差异。为了可读性我把它们拆成三个变量后续完全复用同一段分帧平均代码。env_energy x_bp .^ 2; % 能量解调 env_abs abs(x_bp); % 绝对值解调 env_hilbert abs(hilbert(x_bp)); % 希尔伯特解调后续对每个env_*执行相同的去直流、降采样、分帧平均流程得到三条包络谱。判读时先找四个频点轴频f_s、叶频z*f_sz为叶片数、轴频倍频以及50Hz电源干扰。在P340217这类近场记录中轴频谱峰通常出现在10–25Hz之间倍数关系清晰。三种方法的差异可以用一句话概括能量解调在2*f_s处会出现接近甚至高于基频的谱峰绝对值解调在0.5*f_s附近有杂散分量希尔伯特解调只在f_s及其整数倍处有峰且基频峰高度显著高于背景。有一点很容易忽略三种方法的包络谱纵轴单位不一致直接画在同一张图上时能量解调的大直流分量会把其他谱压成一条水平线。正确的做法是三条谱各自除以自己的最大值再做叠加。4.3 性能差异与适用场景就P340217这段信号来说若带通选在500–3000Hz三种方法都能看到轴频区别在谱峰锐度。希尔伯特解调的半峰宽最窄对轴频估计的稳定性最好。特别是当轴频只有6Hz左右时能量解调的倍频峰会落在12Hz与某些水下机械噪声分量混叠容易误判希尔伯特解调则能把基频孤立出来。如果信号带内混有齿轮箱或泵的强线谱绝对值解调反而有优势。它对相位不敏感不会像希尔伯特解调那样因为两线谱的非线性相加产生拍频。所以我的习惯是先用希尔伯特解调跑一遍如果发现轴频峰旁边出现等间距杂峰再切换能量解调做交叉验证。脚本最后如果保存频谱图注意把三个子图的纵轴范围统一为0–1的归一化幅度。不统一的图会误导后续人工判读。实际项目中我还会把三种方法识别出的轴频写成一个CSV再和AIS或目标运动参数比对验证哪一个峰对应的真实目标。5. 选频段、定帧长让希尔伯特解调在实海数据上稳定复现5.1 带通范围与螺旋桨噪声特征实海数据的带通范围不是固定的。低速漂航目标和高速航行目标的空化噪声频谱重心不同一般中低速目标的信号能量集中在几百Hz到2kHz。拿到新数据后我不会直接套用500–3000Hz而是先画原始信号频谱找能量集中的连续谱区域并避开已知单频干扰。选定后做一次窄带扫描从[fL, fH] [400, 1000]开始每次增加500Hz带宽比较相同参数下希尔伯特包络谱的轴频峰高。峰高最高且背景平坦的频段就是最优带通。5.2 帧长与重叠率的影响帧长直接决定轴频分辨率和谱平均次数。在P340217数据上帧长1s时轴频峰宽约1Hz已经够用若轴频低于5Hz需把帧长加到2s甚至4s。帧长增加会减少可平均帧数10s数据、2s帧长、50%重叠可以平均9帧1s帧长可以平均19帧。我的优先策略是先用最短帧长跑通确认轴频大致位置再逐步加长帧长。如果加长后峰高反而下降说明信号存在慢漂移此时应减少重叠率而不是继续加长。重叠率超过75%时相邻帧高度相关平均带来的增益趋近于零。5.3 用仿真信号校准后处理流程在跑实海数据前最好先用参数完全可控的仿真信号验证整条链路。常见做法是用窄带噪声乘以正弦调制模拟已知轴频的辐射噪声。下面的代码生成轴频8Hz、叶片数4的仿真信号输入到第3章的流程后应该在8Hz和32Hz处看到明显谱峰。fs 50000; t (0:fs*20-1). / fs; noise randn(size(t)); x_sim filter(bpf, noise); % 使用与真实数据相同的带通滤波器 mod_depth 0.6; env_sim 1 mod_depth * sin(2*pi*8*t) .* (0.6 0.4*sin(2*pi*2*t)); x_sim x_sim .* env_sim;把x_sim代入3.2节的希尔伯特解调流程如果主峰不是8Hz优先检查降采样后的频率轴是否算错。特别是fs_env fs / decim这里如果漏写括号轴频位置会偏移很大。校准通过后再切换P340217实际数据能节省大量排错时间。校准图里还应记录轴频峰相对于背景的比值。若比值小于3说明信号太弱或带通选偏需要回到5.1节重新扫描频段。这种量化检查比肉眼看图可靠得多也能在换一条新海试数据时快速判断处理参数是否需要重调。本文还有配套的精品资源点击获取