做水声信号处理的同行应该都遇到过这种事手里拿到一段实测噪声直接做FFT看频谱结果是一片宽宽带带的山包线谱被噪声盖得严严实实有用的信息好像全丢了。但这条船到底什么状态、螺旋桨转多快其实都藏在这段噪声的“包络”里。所谓DEMONDetection of Envelope Modulation On Noise噪声包络调制检测谱分析干的事就是把宽带噪声的幅度包络提取出来再做一次谱分析从中读出螺旋桨的轴频、叶频和叶片数。这是一套非常经典、而且在水声工程里至今仍然高频使用的解调分析流程。这篇文章我打算从原理讲到Matlab实现把DEMON谱分析从输入一段水声信号到输出一张可读的调制谱全流程拆开附带可以直接跑的完整源码和参数计算过程。适合刚入门水声信号处理、或者正在做船舶辐射噪声分析但没理清包络检波细节的同学参考。写完你可以直接把手里的实测.wav数据丢进代码里看看能不能读出螺旋桨的调制特征。1. 从“听声辨船”说起DEMON谱到底在分析什么1.1 螺旋桨空化噪声是怎么来的先捋一捋物理背景不然后面看谱图容易看懵。船舶螺旋桨在高速旋转时桨叶背面的压力会降到水蒸气饱和蒸气压以下水中生成大量空泡这些空泡不断生成、破裂辐射出非常宽的连续噪声。这个噪声覆盖几百赫兹到几十千赫兹听起来就是哗哗的“流噪声”。最关键的一点是空化噪声的强度不是恒定的螺旋桨每转一圈桨叶依次扫过某些位置空化强度就会周期性地起伏。于是宽带噪声的幅度上被叠加了一层周期性“调制”这个调制频率就是螺旋桨的轴频每秒转数乘以叶片数也就是叶频。所以水听器接收到的信号可以简化成这样的模型x(t) [A0 m(t)] * n(t)其中n(t)是宽带空化噪声m(t)是周期性调制函数它包含轴频的基频和一系列谐波尤其是叶频成分。DEMON谱的目标就是把这个藏在宽带噪声里的周期调制m(t)提取出来并做频率分析。这句话很重要——DEMON分析的不是“窄带线谱本身”而是宽带噪声幅度包络的“周期性”。1.2 DEMON谱和LOFAR谱的分工很多初学者会把DEMON和LOFARLow Frequency Analysis Recording搞混。两者的关系我用最粗的方式概括LOFAR谱是对原始信号直接做高分辨率功率谱看的是频域上的窄带线谱比如机械振动传递到水中的轴频线谱、齿轮啮合线谱。DEMON谱是“先解调、再做谱分析”看的是宽带噪声包络里蕴含的调制频率。打个不一定严谨但容易理解的比方LOFAR看的是“声音本身的音高”DEMON看的是“声音响度的起伏节奏”。在船舶辐射噪声分析里这两者是互补的。LOFAR可能在低频段看到轴频的线谱但很容易被环境噪声淹没DEMON则利用了空化宽带噪声能量大的特点抗干扰能力更强。实际工程中常把两张谱结合起来判断目标状态。1.3 为什么“解调”这一步是核心直接对x(t)做FFT那些调制信息是看不见的因为n(t)是宽带随机噪声频谱上本身没有突出的离散谱线。只有把“载波”n(t)去掉才能留下m(t)。怎么去载波最朴素、也是工程上最常用的手段就是平方检波信号平方后载波能量被搬移包络信息里包含的调制频率就会以“差拍”形式出现在低频段。这是整个DEMON分析最基本也最容易被忽略的原理。理解了这一步后面的滤波器设计、降采样、谱估计就都有依据了。2. 方案设计与参数选取为什么这套流程能work2.1 标准DEMON处理流程拆解一套经典的DEMON谱分析流程大概是这样的对原始信号做带通滤波把无用频段滤掉只保留空化噪声能量集中的频段对滤波后的信号做平方检波或者用希尔伯特变换求包络对检波后的信号再做低通滤波只留下调制包络的低频分量低通后因为频率范围大幅缩小可以先降采样减少后续FFT计算量对包络信号去掉直流分量然后分段加窗、做功率谱估计比如Welch平均法在得到的DEMON谱中搜索峰值结合谐波关系识别轴频和叶频。这套流程看起来简单实际调试时每一步都有讲究。接下来我逐个讲参数怎么定、为什么这么定。2.2 滤波器频段和采样率怎么选先说带通滤波。空化噪声的能量分布范围很宽不同航速、不同水深、不同螺旋桨类型能量集中的频段都不一样。仿真实测中选频段要看两个条件一是这个频段内噪声能量要够强这样调制信号的信噪比才高二是要避开明显的强线谱干扰比如机械噪声的窄带线谱。我在仿真里用的3k-10kHz是常见水声频段但拿到实际信号时不要照抄最好先画一版宽带功率谱看一眼能量大致分布在哪个范围再定。然后是低通滤波器的截止频率。调制信号里我们要关注的是轴频和叶频。船舶螺旋桨轴频一般很低从1Hz到几十Hz叶频也就是几十到一两百赫兹的量级。低通截止频率取到2000Hz已经非常宽松了。取宽一点的好处是滤波器阶数不用很高相位畸变小坏处是降采样倍数会受限。如果低通截止设为2kHz按奈奎斯特定律降采样到4kHz就能保证不混叠。这一步定下来后面所有参数就有锚点了。2.3 FFT点数与频率分辨率的关系包络谱要能分辨出轴频和叶频频率分辨率就得够细。Welch法里频率分辨率公式是Δf fs_dem / Nfft比如降采样率fs_dem4000HzNfft取4096那么分辨率差不多是0.98Hz。这个分辨率要分辨8.2Hz的轴频、41Hz的叶频是绰绰有余。假如你的目标螺旋桨转速极低轴频只有2Hz那Nfft要超过8000点相当于窗长2秒以上。有时候你会发现轴频峰值在谱图上看不见先别怀疑算法去算一下分辨率够不够。有个细节是加窗。Welch法分段时必须加窗一般用汉明窗或者汉宁窗。窗函数选择会影响谱泄漏和主瓣宽度对单频调制信号来说差别不大。我这里用hamming比较通用。FFT点数一般取2的幂便于FFT加速但不是必须的非2幂也能算。3. Matlab源码逐步拆解从仿真信号到DEMON谱完整实现3.1 先生成一段带调制特性的水声信号在做真实信号之前我强烈建议先用仿真信号把流程跑通。这样每个参数坏了你都知道是怎么回事因为真实信号你是不知道标准答案的。我这里的仿真思路是生成白噪声用带通滤波器把它做成带宽受限的“空化噪声”再用一个含轴频和叶频的周期性包络去乘它。clear; clc; close all; %% 参数设置 fs 50000; % 原始采样率 50kHz T 10; % 信号时长 10秒 N fs * T; t (0:N-1) / fs; %% 螺旋桨参数标准答案 shaft_freq 8.2; % 轴频 8.2 Hz blade_num 5; % 叶片数 5 blade_freq shaft_freq * blade_num; % 叶频 41 Hz %% 调制包络1 0.6*[cos(轴频) cos(叶频)] m 0.6; % 调制深度 mod_env 1 m * (cos(2*pi*shaft_freq*t) cos(2*pi*blade_freq*t)); %% 宽带空化噪声对白噪声做带通滤波频段3k-10kHz [b_carrier, a_carrier] butter(6, [3000 10000]/(fs/2), bandpass); carrier filter(b_carrier, a_carrier, randn(1, N)); %% 合成水声信号 signal mod_env .* carrier;这里的调制深度m0.6是什么意思就是包络起伏幅度是载波平均幅度的60%在实测中这个值通常不会太高低的时候可能只有0.2甚至更小。调制深度越低包络谱的峰值就越不明显这是DEMON分析的一个天然痛点。仿真里用0.6是为了先让你看清谱线长什么样。3.2 带通滤波与包络检波的核心代码接下里进入DEMON主流程。注意我这里带通滤波和仿真里的载波生成用了相同频段这是故意为之——仿真信号本来就是在这个频段内生成的。真正分析实测信号时带通频段要靠前面的功率谱预估来定。滤波后用平方检波再做一次低通提取出低频包络。%% DEMON分析主流程 % 1. 带通滤波锁定空化噪声主要频段 [b_bp, a_bp] butter(4, [3000 10000]/(fs/2), bandpass); x_bp filter(b_bp, a_bp, signal); % 2. 平方检波 x_sq x_bp .^ 2; % 3. 低通滤波保留调制包络截止2000Hz [b_lp, a_lp] butter(4, 2000/(fs/2), low); x_env filter(b_lp, a_lp, x_sq);平方检波之后信号里混着直流分量、低频调制分量、还有高频残余。低通滤波一方面把二次检波产生的高频项滤掉另一方面把信号带宽压到2kHz以下方便后面降采样。如果你不低通就直接降采样高频分量折叠回低频段会在DEMON谱里造成一堆假的谱峰查起来非常头疼。3.3 降采样与Welch谱估计这时信号仍然在50kHz采样率上但有效带宽只有2kHz完全没有必要保留50k的采样率。降采样到4kHz数据量缩小12.5倍FFT算起来快得多而且不影响结果。我用resample完成降采样它内部带了抗混叠滤波比直接抽值靠谱。% 4. 降采样到4000Hz fs_dem 4000; x_dem resample(x_env, fs_dem, fs); % 5. 去掉直流分量 x_dem x_dem - mean(x_dem); % 6. Welch功率谱估计 nfft 4096; % 频率分辨率 4000/4096 ≈ 0.98Hz win hamming(nfft); noverlap 2048; % 50%重叠 [Pxx, f_dem] pwelch(x_dem, win, noverlap, nfft, fs_dem); %% 绘制DEMON谱 figure(Color, w); plot(f_dem, 10*log10(Pxx), LineWidth, 1.2); xlim([0 150]); xlabel(调制频率 (Hz)); ylabel(功率谱密度 (dB)); title(DEMON谱仿真信号轴频8.2Hz叶频41Hz); grid on;运行之后你应该能在8.2Hz处看到一根非常明显的谱峰41Hz处看到第二根谱峰16.4Hz左右还有一根轴频的二次谐波。多数情况下轴频基频是包络谱里最高的峰叶频次之。但注意平方检波会产生谐波和交叉项比如5次谐波、轴频与叶频的和差项这些在谱图上都会出现并不代表真实的螺旋桨特征。你看图的时候要会判断真正的轴频/叶频峰值会有谐波关系而且频率都比较规整。3.4 封装成可直接复用的函数工程上不建议把一长串脚本到处复制我习惯把DEMON分析封装成一个函数。参数用结构体传入方便批量处理和调参。function [freq, demon_spectrum] demon_spectrum_analysis(x, fs, params) % DEMON谱分析主函数 % 输入 % x : 输入水声信号行向量 % fs : 原始采样率 % params : 结构体包含以下字段 % freq_range - 带通滤波频段 [fL fH] % dem_fs - 降采样率 % nfft - FFT点数 % noverlap - 重叠点数 % 输出 % freq : 频率轴Hz % demon_spectrum : DEMON功率谱线性刻度 if nargin 3 || isempty(params) params.freq_range [3000 10000]; params.dem_fs 4000; params.nfft 4096; params.noverlap 2048; end % 1. 带通滤波 [b_bp, a_bp] butter(4, params.freq_range/(fs/2), bandpass); x_bp filter(b_bp, a_bp, x); % 2. 平方检波 x_sq x_bp .^ 2; % 3. 低通滤波截止频率设为降采样率的一半 f_lp params.dem_fs / 2; [b_lp, a_lp] butter(4, f_lp/(fs/2), low); x_env filter(b_lp, a_lp, x_sq); % 4. 降采样 x_dem resample(x_env, params.dem_fs, fs); x_dem x_dem - mean(x_dem); % 5. Welch功率谱 win hamming(params.nfft); [demon_spectrum, freq] pwelch(x_dem, win, params.noverlap, params.nfft, params.dem_fs); end这个函数用起来很简洁我后面所有实验都是基于这个函数改的。你拿到实测wav文件后用audioread读进来然后把采样率和信号丢进去基本就能出一版结果。当然真实数据的带通频段、降采样率肯定要手动调。4. 常见问题与调试实录峰值没了调制线找不着4.1 问题速查表我在给项目调试DEMON谱时踩过不少坑整理成一张速查表按“现象-原因-解法”的顺序来方便你对照排查现象可能原因处理方法DEMON谱上全是低频缓慢起伏没有尖峰带通滤波频段选错调制信号被滤掉先画原始信号宽带谱看能量集中在哪再设freq_range轴频处看不到峰但叶频能看到FFT分辨率不够轴频低于Δf增大nfft或降低降采样率轴频低于2Hz时考虑nfft8000峰值一大堆不知道哪个是轴频平方检波产生了大量谐波和交叉项寻找候选峰之间的最小公约数频率结合叶片数先验判断谱图右侧有一片高能量平台低通滤波不够狠没有滤干净高频残余降低低通截止频率或提高滤波器阶数0Hz附近能量极大把低频峰盖住去直流没做好x_dem x_dem - mean(x_dem)一定要在谱估计之前做最容易迷惑人的就是“峰值一大堆”这种情况。平方检波本质是非线性运算它会将调制包络里的各频率分量两两组合成和差项所以你在DEMON谱上看到轴频、叶频的同时还会看到2倍轴频、轴频叶频、叶频-轴频之类的峰。不要一看到峰就都当成螺旋桨特征要学会用谐波关系的约束去筛选。4.2 如何稳定锁定轴频和叶频工程上我一般不会只靠眼睛看图会在程序里写一个自动峰值检出的步骤。思路不复杂先用findpeaks找出功率谱里幅度最高的几个候选峰然后对所有候选峰做两两比值判断看它们是否满足整数倍关系。满足整数倍关系的最小频率基本就是轴频。% 峰值检测在dB域找最高的10个峰 Pxx_db 10*log10(Pxx); [~, locs] findpeaks(Pxx_db, SortStr, descend, NPeaks, 10); cand_freqs sort(f_dem(locs)); % 谐波关系校验检查候选峰之间是否存在整数倍关系 for ii 1:length(cand_freqs) for jj ii1:length(cand_freqs) if cand_freqs(ii) 0.5 continue; end ratio cand_freqs(jj) / cand_freqs(ii); % 比值接近整数认为满足谐波关系 if abs(ratio - round(ratio)) 0.03 fprintf(候选轴频 %.2f Hz对应谐波 %.2f Hz\n, ... cand_freqs(ii), cand_freqs(jj)); end end end这个比例阈值0.03是经验值。对于低速转动的螺旋桨轴频可能只有几赫兹频率分辨率有限比值偏差会大一点。阈值给得太小会漏检给太大会误报需要根据你的Nfft和信号时长动态调整。另外findpeaks是Signal Processing Toolbox里的函数如果没有工具箱可以用简单的滑动窗极大值检测代替——遍历频谱每个点和左右相邻的k个点比较比所有邻居都大就记为候选峰。4.3 从仿真到实测拿到真实数据会遇到的额外问题仿真信号跑通之后第一次处理实测数据时大概率会被现实教育一遍。常见的情况有这么几类第一实测水声信号里不止一条船可能有多个辐射源调制包络会互相叠加。这种情况下DEMON谱会出现好几簇峰值一簇对应一条船的螺旋桨特征。想区分它们不容易通常需要结合LOFAR谱和波束形成做空间滤波先让目标信号在时域上干净一些再进DEMON。第二环境噪声里的瞬态脉冲比如生物click声、雨噪声、冰裂声这些脉冲在平方检波后会产生宽频段的冲击拉高整个包络谱的背景把周期调制峰盖住。处理办法是在进入DEMON之前先做一个“去瞬态”预处理把幅度远超背景的样本点用中值替代或整段剔除。第三数据长度不够。包络谱要看到轴频至少需要包含几个轴频周期。轴频8.2Hz时周期约0.12秒10秒数据有80多个周期Welch平均后谱线很稳。但如果轴频只有1.5Hz10秒数据只有15个周期谱估计方差就会偏大。采样时间能长尽量长或者减少Welch平均段数来换取频率分辨率。5. 还能怎么继续拓展这个DEMON工程5.1 自动识别与谐波校验的工程化思路上面给的谐波校验其实还比较粗糙工程上有两个更稳的做法。第一个是在频域做“梳状滤波器扫描”假设一个候选轴频f0计算f0、2f0、3f0……每个位置的功率和作为“该轴频可信度”。扫描f0从0.5Hz到20Hz可信度最高的f0就被判为轴频。这种方法利用了多次谐波的相关性比只看单峰稳健得多。第二个是对数域谱背景归一化。用滑动中值滤波估计整个频谱的背景噪声再用每个频点的功率除以背景值得到一个“谱显著度”曲线。在显著度曲线上做峰值检测能有效压制环境噪声引起的伪峰。特别是海况不好、背景噪声很重的时候这一招能救回很多原本要被淹没的弱调制峰。5.2 与机器学习结合的方向最近几年有很多人把DEMON谱当作特征输入分类器用于水下目标识别。但我想提醒一句DEMON谱本身对噪声非常敏感同一个目标在不同航速下轴频叶频都会变直接拿全频段向量去训练很容易过拟合。更合理的做法是先自动提取出轴频、叶频、峰值显著度、调制深度这几个标量特征再做分类。标量特征的物理含义明确泛化能力比原始谱向量好得多。如果你感兴趣可以在这个demo基础上加一个特征提取环节。5.3 实时处理要点如果想把DEMON分析做成流式实时处理遇到的主要矛盾是“更新速率”和“频率分辨率”不可兼得。要看到轴频频谱窗至少要有两三个轴频周期的长度这就决定了最短延迟。解决办法是分双链路一条链路用短窗做快速更新监控调制谱的大致变化另一条链路用长窗做精细DEMON谱得到稳定的轴频叶频估计。这种架构在工程里比较常见。另外实时场景下滤波器的计算量也要考虑。butter滤波器是IIR型阶数不高时运算量很小。如果前端信号流是16kHz采样率但带通在3k-10kHz那数据里其实含有超过奈奎斯特频率的成分必须先用抗混叠低通把采样率降下来再进流程。我见过有人直接把16k采样率的信号拿来做DEMON然后发现带通上限设为8k以下参数怎么调都不对就是这个原因。写在最后一点实际体会整套DEMON流程写下来代码量其实不超过80行但每行背后都有物理含义和工程取舍。我在调试时最深的体会是不要一上来就追求“花哨”算法先把带通、检波、低通、谱估计这条基础链路每一级的输出波形都画出来看一遍。哪个环节信号爆了、哪个环节谱线被滤掉了一眼就能看出来比反复调一个神秘参数高效得多。另外仿真信号只能帮你验证流程真正的问题永远在实测数据里。把这条链路跑熟之后建议尽快拿一段真实水声录音来试试。哪怕结果不完美、峰值不明显也比仿真跑一百遍学到的东西多。如果你调通了对着弹出来的DEMON谱看到那个小小的轴频峰安静地立在那里那一刻你会觉得前面踩过的坑都值了。