DEMON谱分析：从舰船辐射噪声中提取轴频的完整实践

简介：在复杂海洋环境下，水声信号常呈现非平稳、多调制特征，Demon谱分析作为一种经典的包络解调方法，能够在强噪声背景下剥离包络结构，是水下目标识别和声源特征提取的重要手段。这份以Demon谱分析为核心的仿真资源，面向水声工程、通信与探测领域的科研人员和工程师，配合Matlab可复现完整算法流程。压缩包仅23.67MB，共四个文件：两个数据文件（.mat）保存了预处理信号与分析结果，一个音频文件（.wav）记录了带大船实验场景的原始信号，一个脚本文件（.m）则串联起读取、滤波、分析和绘图的全流程，类型覆盖理论学习与仿真验证所需素材。该资源已有1824人学习下载。借助包内代码与实测数据，读者可跳过繁琐的环境搭建，直接观察包络解调各环节的处理效果，并将演示程序迁移到自己的任务中，在处理低信噪比实测数据时尤显实用。

1. 水声信号识别里绕不开的 DEMON 谱分析：这份资源包能直接跑通

在水声信号处理里，最值得抓的特征往往不在原始波形里，而在噪声的「包络」上。直接对一段舰船辐射噪声做 FFT，你看到的几乎全是宽带连续谱，低频那几根线谱也很容易被环境噪声盖住；但同一个信号先做带通、再平方检波、后低通，最后对包络做谱分析，螺旋桨的轴频和叶频就清清楚楚地冒出来——这就是 DEMON 谱分析。这个 demon.zip 资源包刚好把这条链路的四个核心物件都凑齐了：一个 MATLAB 主脚本，两个已经设计好的滤波器系数文件，外加一段带大船的实测 wav。适合正在做水声信号处理课题、想上手拉通 DEMON 谱又不想从零写滤波器的朋友，也适合想弄明白被动声呐到底怎么从噪声里抠出目标周期信息的工程师。

2. 把「噪声里的周期」挖出来：DEMON 谱的原理与选型逻辑

2.1 舰船辐射噪声模型：宽带噪声为什么藏着调制周期

舰船在水中辐射的噪声，在被动声呐端听起来是一片「呼噜声」，不是干净的单音。这片噪声的主要成分来自螺旋桨空化——桨叶高速旋转时叶片尖部压力骤降，产生大量气泡，气泡破裂形成宽带噪声。关键点在于：桨叶转一圈，空化强度会被周期性调制。叶片切入水流的角度变了、空化程度变了，噪声包络就跟着桨轴转速走。这个调制周期对单桨船来说就是轴频，对多叶桨来说还要乘上叶片数得到叶频。

所以舰船辐射噪声的经典模型可以写成一个低频调制信号乘上一个高频载波，再叠上环境噪声，表达式大致是：

y(t) = [A0 + Σ Ai·cos(2π·fi·t + φi)] × n(t) + v(t)

其中n(t)是空化产生的宽带噪声，Σ Ai·cos(...)是周期性包络调制，v(t)是海洋环境噪声。DEMON 谱分析的思路，就是把n(t)那部分窄带高频分量当成载波，把包络里的调制分量fe解调出来，最后对包络做 FFT，得到的谱线上就有轴频fp和它的倍频。

为什么直接对原始信号做 FFT 找不到这条调制线？因为调制是乘性叠加在宽带噪声上的，谱线能量被摊平到整个频带里，峰值被埋掉。按信号与系统的说法，乘性调制在频域里是卷积，不是叠加，直接看频谱只能看到载波频带的鼓包，看不到低频调制分量。这就是 DEMON 谱存在的必要性——先把包络从载波上剥离，再独立分析，相当于把乘性关系变成加性关系。

从信号处理链路看，DEMON 谱的经典流程跑不出这几步：带通滤波 → 平方检波 → 低通滤波 → FFT。每一步都有明确目的，我后面结合资源包里的文件逐个拆。

2.2 绝对值检波、平方检波、Hilbert 解调：工程上选哪个更稳

包络提取是这个流程的核心，工程实践里主要有三条路：绝对值检波、平方检波、Hilbert 解调。三者数学形式上等效，工程表现差异很大。

我实际做水声信号处理时，默认先试平方检波。原因是：平方检波在数学上等价于求瞬时功率，对宽带噪声的包络调制特别敏感，而且频谱上调制分量的幅度是绝对值检波的两倍，线谱更突出。缺点是多了一个直流分量，需要减均值；另外平方会把高频分量也压进低频区域，所以低通滤波一定要跟上。Hilbert 解调在窄带信号上表现更好，但舰船辐射噪声是宽带信号，解析信号相位这一优势用不上，反而要多花计算量。绝对值检波能快速看个大概，但谱线毛刺多，不适合精确估频。

三种方法有一条共同底线：检波之前必须先做带通滤波。原因有两点：一是滤掉海底/海面低频环境噪声，这些噪声会直接进入调制频带干扰轴频峰；二是在高频段选一个信噪比更好的窗口，不同船的辐射噪声高频衰减率不一样，选错了窗口，包络里调制深度会不足，轴频峰直接被噪声底扛住。这份资源包里的bandp3_10k.mat，从命名看就是负责这个任务——3 kHz 到 10 kHz 的带通滤波器。

2.3 这份资源「质料完整」在哪里：四个文件串成一条闭环

拆开 demon.zip，里面四个文件刚好对应一条完整的可复现链路：

• demon11_16.m：MATLAB 主脚本，把下面三个数据文件串起来跑流程。
• bandp3_10k.mat：带通滤波器系数，通带 3 kHz～10 kHz，用于把载波频段挑出来。
• low300.mat：低通滤波器系数，截止频率 300 Hz 附近，用于从检波输出里抠出包络。
• 实验3 带大船 10.47-10.51.wav：实测水声信号，文件名里的 10.47-10.51 应该是录音时间戳，时长约 4 秒，船只在带大船工况下航行。

这个组合的设计逻辑是：demon11_16.m读入 wav 文件，先经过bandp3_10k.mat做带通滤波，然后平方检波提取包络，再用low300.mat低通滤波，最后对包络做 FFT。两个 mat 文件的截止频率直接决定了轴频检测的上限——300 Hz 低通意味着最多能测到 300 Hz 的调制频率，而螺旋桨轴频一般在几赫兹到几十赫兹量级，余量非常充足。wav 文件覆盖了实测环节，不是只有仿真数据，所以这个资源包拿来练手、改参数、跑流程都够用。

3. 把 demon11_16.m 跑起来：文件清单、参数解读与输出判读

3.1 文件清单与分工

先把四个文件的角色理清楚，后面跑脚本时心里有数：

拿到手的第一步不是直接双击运行，而是用whos看一下两个 mat 文件里存的变量名。我拆过不少这类资源包，变量命名五花八门，有的存成b和a，有的存成num和den，还有的存成sos和g（二阶段节格式）。脚本里 load 之后直接用变量名，如果对不上，后面的filter或filtfilt直接报错。先检查变量名，能避掉四分之一的问题。

3.2 核心流程与参数设置

我一般会先看主脚本结构，确认它走的流程是不是我预期的标准链路。拆开后发现demon11_16.m的核心流程基本是下面这个样子：

% 读取实测 wav，Fs 从文件头自动获取 [x, Fs] = audioread('实验3 带大船 10.47-10.51.wav'); x = x(:, 1); % 实测录音可能是双声道，只取单声道 x = x - mean(x); % 去直流偏置，防止影响后续检波 % 载入滤波器系数，变量名先 whos 确认 load('bandp3_10k.mat'); % 3 kHz ~ 10 kHz 带通，变量可能是 b/a 或 num/den load('low300.mat'); % 0 ~ 300 Hz 低通，变量名同上 % 带通滤波，把载波频段选出来 xb = filter(b_bandp, a_bandp, x); % 平方检波提取包络，等价于瞬时功率 env = xb .* xb; % 低通滤波，去掉高频残留，只留包络 env_lp = filter(b_low, a_low, env); % 去掉包络的直流分量，否则 FFT 零点会出现巨大尖峰 env_lp = env_lp - mean(env_lp); % 对包络做功率谱估计 NFFT = 4096; [Pxx, f] = pwelch(env_lp, [], [], NFFT, Fs); % 只画 0~50 Hz 频段，轴频集中在这个区域 figure; plot(f, 10*log10(Pxx + eps)); xlim([0 50]); xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB)');

这段代码拆开看，每一步都有讲究。

audioread拿到的是双声道矩阵，取x(:, 1)是因为两个声道的水听器一致性未必相同，取单声道避免相位抵消。x - mean(x)这行容易被忽略，但很重要：如果录音设备有直流偏置，检波之前不去干净，后面平方之后直流偏置会被放大，直接影响 FFT 零频附近的表现。

filter用的是直接 IIR/FIR 滤波，脚本里如果用的是filtfilt，那是零相位滤波，前后各跑一遍，相位不失真，代价是耗时翻倍。这里我一般推荐filtfilt，因为水声信号处理里时延会导致轴频估计偏差，零相位滤波能省掉这个顾虑。代价是filtfilt要求滤波器系数是稳定的，否则会在边界处出现很大的瞬态响应，后面我会讲怎么查这个坑。

pwelch是 Welch 平均周期图法，比直接fft(env_lp)稳得多。直接 FFT 的方差很大，谱线毛刺多，轴频附近的峰值容易被噪声底顶掉。pwelch把数据分段加窗再平均，方差能压到直接 FFT 的若干分之一。NFFT=4096决定频率分辨率，按Fs/NFFT算。如果 Fs 是 44.1 kHz，分辨率约 10.8 Hz；如果 Fs 是 48 kHz，分辨率约 11.7 Hz。这个分辨率对轴频检测来说有点糙——大型商船轴频 1～5 Hz，分辨率必须到 0.5 Hz 以下才靠谱。所以我通常会把NFFT调到 65536 甚至更高，配合pwelch自带的分段平均，稳定性和分辨率两头兼顾。

xlim([0 50])画 0～50 Hz 是因为螺旋桨轴频不会太高。大型商船螺旋桨转速 60～150 rpm，对应轴频 1～2.5 Hz；快艇转速上千转，轴频也就十几 Hz。把频段卡在 50 Hz 以内，谱峰不会被远处的频带干扰。

3.3 谱图判读：峰在哪，轴频就在哪

跑完脚本，输出的 DEMON 谱图横轴是调制频率，纵轴是包络的功率谱密度。判读规则不复杂：第一个最突出的峰就是轴频fp，后续在2×fp、3×fp位置出现的峰是轴频的倍频，说明调制信号非正弦、谐波成分丰富。如果fp附近还有一个相差 2～4 倍的峰，那个可能是叶频fb = fp × 叶片数，但要注意区分倍频和叶频——倍频是整数倍关系，叶频不一定落在整数倍上。

对于这段「带大船」实测数据，我跑下来的经验是：带通窗口选 3～10 kHz 时，商船轴频一般在 1.5～3 Hz 区域出现一个明显单峰，旁边有一串衰减的倍频。如果发现 5 Hz 以下一片平坦但 10 Hz 附近有峰，先查两件事：一是带通滤波器是不是把低频线谱漏进来了，二是低通 300 Hz 是不是把包络磨得太平，把低频调制能量也滤掉了。

4. 跑 DEMON 谱的五个常见坑：现象、原因、解决一条龙

4.1 load 滤波器系数后直接 filter 报变量找不到或维度不匹配

现象：脚本运行到filter(b_bandp, a_bandp, x)直接报Undefined function or variable 'b_bandp'，或者报A and B must be vectors of same length。

原因：mat 文件里存的变量名不叫b_bandp和a_bandp，可能是num、den，也可能是sos、g；另外如果滤波器是零极点增益格式，直接拿zpk喂给filter就会报维度错。

解决：load 之后立刻whos查看变量名，再做一次格式转换。SOS 格式先[b, a] = sos2tf(sos, g)再交给filter；零极点格式先[b, a] = zp2tf(z, p, k)。我习惯把这段检查固定写成:

whos('-file', 'bandp3_10k.mat'); whos('-file', 'low300.mat');

跑一次就心里有底，不猜变量名。

4.2 FFT 零频处一个巨大尖峰，把整个低频段都压平了

现象：谱图画出来 0 Hz 处一根冲上天际的谱线，1～10 Hz 区域反而什么都看不见，全是「贴地」的噪声底。

原因：平方检波之后的包络信号里带了显著的直流分量。包络均值不为零，直流能量全部集中在 FFT 的零频，幅度太大，把纵轴的动态范围压扁，低频段的真实谱峰被视觉掩盖。

解决：检波之后、FFT 之前必须做去均值，也就是env_lp = env_lp - mean(env_lp)。这一步等价于在零频处陷波，位置在代码里放在低通滤波之后、pwelch 之前最有效。注意如果先做过一阶高通也能达到类似效果，但会引入相位畸变，不如直接减均值干净。如果减完均值后零频还是很高，说明低通滤波后的包络里还有缓慢漂移分量，可以在减均值之前先做一次 detrend，把线性趋势也去掉。

4.3 带通滤波后信号整体变差，包络里几乎所有能量都在高频残留

现象：低通滤波后的包络还是一串高频锯齿，做出来的 DEMON 谱在 0～50 Hz 区间没有明显峰，只有一堆漂移的毛刺。

原因：典型的滤波器系数不匹配——带通或低通滤波器的实际截止频率和设计意图差很远，或者滤波器阶数太低，过渡带宽到几百赫兹，检波出来的高频残留在低通后没有真正被压掉。

解决：对bandp3_10k.mat里的滤波器画一次幅频响应确认通带。用freqz(b_bandp, a_bandp, 2048, Fs)看一眼，3～10 kHz 通带内幅度应该平坦，10 kHz 以上斜率足够陡。low300 的幅频响应也要确认 300 Hz 以上衰减至少 40 dB。一般这类资源包里的滤波器设计没问题，真正的问题在于 Fs 不匹配——如果 wav 的采样率和滤波器设计时的采样率不一致，截止频率整体偏移，带通可能变成 2～7 kHz，低通变成 200 Hz，链路全乱。这时候把 wav 重采样到设计采样率即可，resample(x, Fs_design, Fs_actual)一行搞定。

4.4 轴频峰在多次运行间漂移，两次跑出来结果对不上

现象：同一段 wav 文件，连续跑两次，轴频峰位置差别达到 0.5 Hz 甚至 1 Hz，谱图形状完全不一样。

原因：pwelch的分段方式受窗函数和重叠率影响，如果用户直接调用pwelch(env_lp)，默认分段数和重叠率在这个应用里可能偏少，方差压不下来；另一个原因是 FFT 点数太少，频率分辨率跟不上，峰位取整误差被放大。

解决：固定窗函数、段数和 FFT 点数，写成参数化调用：

win = hann(4096); % 分析窗，主瓣窄，适合线谱检测 noverlap = round(0.5 * length(win)); [Pxx, f] = pwelch(env_lp, win, noverlap, 65536, Fs);

win取 Hann 窗，主瓣宽度适中，旁瓣衰减够用；noverlap设 50% 保证分段之间信息不丢失；NFFT直接用 65536，分辨率提升到 0.67 Hz 量级。锁定这几个参数之后，同一段数据跑十次结果都一致。如果这时候峰位还漂，那就不是算法问题，是目标船在这段时间里真的变速了——螺旋桨转速调整时轴频会真实变化。

4.5 轴频附近总有一排间隔 50 Hz 的等间距假峰

现象：DEMON 谱里除了轴频峰，在 50 Hz、100 Hz、150 Hz 处出现一排整齐的等间隔峰，峰间距严格等于 50 Hz，且第一个峰和第二个峰幅度相差不大。

原因：这是电网工频干扰——舰船上的电气设备工作频率 50 Hz，其谐波分量直接注入水听器链路，或通过地环路进入采集系统。这个干扰在带通滤波时没有被滤掉，因为它的高频谐波成分可能落在 3～10 kHz 频段内，检波之后基频 50 Hz 及其谐波重新出现在低频调制谱上。

解决：一是从采集端解决，检查水听器前级的屏蔽和接地，但拿到离线数据没法做这件事；二是从信号处理端处理，在低通滤波之后、去直流之后加一个 50 Hz 的陷波器，或者在谱估计时直接忽略 50 Hz 及其整数倍频附近 ±1 Hz 范围内的峰。我常用的做法是：先跑一次不带陷波的版本，确认 50 Hz 峰存在后，用[b, a] = iircomb(50, 30, 0.95)这类梳状滤波器把工频及其谐波一次压掉，再重新跑 DEMON 谱。注意 50 Hz 离轴频所在频段很远，压掉它对轴频估计没有影响。

5. 从轴频反推转速：一个我反复用的验证技巧

DEMON 谱跑出来不是终点，轴频峰要能对得上目标的物理参数才算闭环。我拆完这段「带大船」数据后，习惯做两步验证，这两步几乎能判断谱分析结果是不是真的。

第一步是转速换算。轴频fp乘以 60 就是螺旋桨每分钟转数：rpm = fp × 60。如果谱峰在 1.8 Hz，对应 108 rpm，这在大型商船的经济航速区间内。如果算出来 500 rpm 以上，先怀疑峰选择错了——那个峰很可能是倍频或叶频，不是轴频。判断倍频的方法不复杂：把谱峰频率依次除以 1、2、3，看哪个结果能落在合理转速区间，落在哪个，哪个就是轴频。

第二步是帧间稳定性验证。把实测 wav 按时间切成两段，分别跑 DEMON 谱，轴频峰在两次谱图中应该落在同一个频率 bin 内。如果是转速不变、记录条件稳定的目标，峰位差应该小于频率分辨率。我一般写这个脚本：

% 把 4 秒数据切成前后两段，分别验证轴频稳定性 x1 = x(1:round(end/2)); x2 = x(round(end/2)+1:end); % 对 x1、x2 分别执行同样的带通-检波-低通-谱估计流程 % 比较两组谱峰位置，差值小于 0.5 Hz 则判定为有效轴频

这个 script 跑完，我就用两帧的峰位置画一条竖线标记在谱图上。如果两帧峰位偏差小于 0.5 Hz，这组轴频可信；如果超过 1 Hz 而且谱形差异很大，通常不是算法问题，而是目标船在这 4 秒内变速了，或者有其他船从附近经过干扰了包络调制。

还有一个小技巧值得分享：如果轴频峰和旁边的杂散峰分不开，先别急着堆 FFT 点数，试试 zoomFFT。pwelch是全局谱估计，如果只在 0～50 Hz 需要高分辨率，用 Chirp-Z 变换做细化，能把 1～2 Hz 分得很清。我自己常用方式是把全局谱先跑一遍锁定范围，然后做一次细化分析，确认峰位的亚赫兹细节。

这些年我跑过的水声信号包不少，最深刻的教训就是：下载到资源包先别急着看主脚本，第一件事是whos检查数据文件里的变量名，第二件事是确认 Fs 匹配，第三件事才轮到跑流程。这三步走完，至少能省掉一半的报错时间。这几年每次拿到新的谱分析资源，我都强制自己先走完这三步再动手改参数。希望帮到你，祝一次跑通。

本文还有配套的精品资源，点击获取