DEMON谱分析原理与Matlab实现：从水声噪声中提取螺旋桨轴频叶频

做水声信号处理的同行应该都遇到过这种事：手里拿到一段实测噪声，直接做FFT看频谱，结果是一片宽宽带带的山包，线谱被噪声盖得严严实实，有用的信息好像全丢了。但这条船到底什么状态、螺旋桨转多快，其实都藏在这段噪声的“包络”里。所谓DEMON（Detection 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和LOFAR（Low Frequency Analysis Recording）搞混。两者的关系我用最粗的方式概括：

• LOFAR谱是对原始信号直接做高分辨率功率谱，看的是频域上的窄带线谱，比如机械振动传递到水中的轴频线谱、齿轮啮合线谱。
• DEMON谱是“先解调、再做谱分析”，看的是宽带噪声包络里蕴含的调制频率。

打个不一定严谨但容易理解的比方：LOFAR看的是“声音本身的音高”，DEMON看的是“声音响度的起伏节奏”。在船舶辐射噪声分析里，这两者是互补的。LOFAR可能在低频段看到轴频的线谱，但很容易被环境噪声淹没；DEMON则利用了空化宽带噪声能量大的特点，抗干扰能力更强。实际工程中常把两张谱结合起来判断目标状态。

1.3 为什么“解调”这一步是核心

直接对x(t)做FFT，那些调制信息是看不见的，因为n(t)是宽带随机噪声，频谱上本身没有突出的离散谱线。只有把“载波”n(t)去掉，才能留下m(t)。怎么去载波？最朴素、也是工程上最常用的手段就是平方检波：信号平方后，载波能量被搬移，包络信息里包含的调制频率就会以“差拍”形式出现在低频段。这是整个DEMON分析最基本也最容易被忽略的原理。理解了这一步，后面的滤波器设计、降采样、谱估计就都有依据了。

2. 方案设计与参数选取，为什么这套流程能work

2.1 标准DEMON处理流程拆解

一套经典的DEMON谱分析流程大概是这样的：

1. 对原始信号做带通滤波，把无用频段滤掉，只保留空化噪声能量集中的频段；
2. 对滤波后的信号做平方检波（或者用希尔伯特变换求包络）；
3. 对检波后的信号再做低通滤波，只留下调制包络的低频分量；
4. 低通后因为频率范围大幅缩小，可以先降采样，减少后续FFT计算量；
5. 对包络信号去掉直流分量，然后分段加窗、做功率谱估计（比如Welch平均法）；
6. 在得到的DEMON谱中搜索峰值，结合谐波关系识别轴频和叶频。

这套流程看起来简单，实际调试时每一步都有讲究。接下来我逐个讲参数怎么定、为什么这么定。

2.2 滤波器频段和采样率怎么选

先说带通滤波。空化噪声的能量分布范围很宽，不同航速、不同水深、不同螺旋桨类型，能量集中的频段都不一样。仿真实测中选频段要看两个条件：一是这个频段内噪声能量要够强，这样调制信号的信噪比才高；二是要避开明显的强线谱干扰，比如机械噪声的窄带线谱。我在仿真里用的3k-10kHz是常见水声频段，但拿到实际信号时不要照抄，最好先画一版宽带功率谱看一眼能量大致分布在哪个范围再定。

然后是低通滤波器的截止频率。调制信号里我们要关注的是轴频和叶频。船舶螺旋桨轴频一般很低，从1Hz到几十Hz，叶频也就是几十到一两百赫兹的量级。低通截止频率取到2000Hz已经非常宽松了。取宽一点的好处是滤波器阶数不用很高，相位畸变小，坏处是降采样倍数会受限。如果低通截止设为2kHz，按奈奎斯特定律，降采样到4kHz就能保证不混叠。这一步定下来，后面所有参数就有锚点了。

2.3 FFT点数与频率分辨率的关系

包络谱要能分辨出轴频和叶频，频率分辨率就得够细。Welch法里频率分辨率公式是：

Δf = fs_dem / Nfft

比如降采样率fs_dem=4000Hz，Nfft取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;

这里的调制深度m=0.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谱上看到轴频、叶频的同时，还会看到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 = ii+1: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.5Hz，10秒数据只有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谱，看到那个小小的轴频峰安静地立在那里，那一刻你会觉得前面踩过的坑都值了。