简介Demon谱分析是水声信号处理中的一种重要方法常用来揭示幅度调制与包络结构适用于舰船辐射噪声、海洋环境噪声等非平稳信号的分析与特征提取。这份资料以Matlab仿真为主线面向水声工程科研人员、信号处理学习者以及水下目标探测相关工程师提供了从数据读取、预处理、包络提取到谱计算与可视化的完整算法实现。压缩包共包含4个文件其中有可运行的Matlab脚本、参数与中间结果数据文件以及一段真实实验音频整体约23.67MB结构紧凑便于逐项学习。目前已有1824人浏览学习作为专题资源具备一定参考价值。内容涉及类似“带大船”实测场景的音频和不同频段的处理数据配合源码可复现完整流程帮助理解调制结构分解、频带选择及特征参数获取的思路同时代码留有改造空间可直接用于科研分析或课程设计。1. 先弄清楚DEMON谱分析在解决什么问题从水声噪声里挤出螺旋桨转速做水声信号处理的人几乎都碰到过同一个需求拿到一段含舰船或水下航行器噪声的录音想知道它的螺旋桨转速、叶片数甚至据此做目标识别。直接对原始信号做FFT往往什么也看不出来因为螺旋桨噪声是宽带信号能量铺在几百赫兹到几千赫兹。demon.zip这类资料包之所以在从业者里流传本质是它把“demon谱分析”的完整流程——仿真信号生成、包络检波、谱线提取——打包成了能直接复现的东西。但你不必等拿到它才动手这里就用一段可运行的仿真代码把整个流程走通。这个标题里的“水声信号处理”和“水声信号”两个关键词指的正是这类受海洋信道强烈干扰的低频宽带信号而DEMON谱分析是目前从水声信号里提取螺旋桨周期特征最成熟的手段之一。2. 为什么螺旋桨噪声里藏着转速DEMON谱的物理模型与仿真依据2.1 空化噪声的调制机理轴频和叶频凭什么混在宽带信号里螺旋桨高速旋转时叶片背面压力降低水发生空化产生大量气泡。这些气泡随压力波动不断生成和溃灭辐射出连续的宽带噪声频率范围可以从几十赫兹一直延伸到几十千赫兹没有明显线谱。直接做FFT只能看到一条随频率缓慢变化的功率谱包络轴频的信息完全被淹没这就是为什么普通频谱分析在这个场景里失效。关键在“周期性门控”每个叶片扫过同一角度时空化强度都经历一次相似的变化于是整段宽带噪声的幅度被周期性地增强再减弱。对一台叶片数为Z的螺旋桨幅度调制里同时包含轴频f_s每转一周所有叶片共同形成一次强度起伏和叶频Z·f_s每个叶片各自通过某一固定方位造成的周期。商船螺旋桨多见4到7叶轴频往往只有2到15赫兹叶频则在几十赫兹以内。例如轴频2.5Hz、六叶桨的叶频就是15Hz二次谐波落在30Hz。这个频段远低于载波所在的宽带频段所以整个信号在频域上呈现为“高频宽带载波低频幅度调制”的复合结构。DEMON谱分析Detection of Envelope Modulation On Noise噪声包络调制检测正是为这种结构设计的先把某个窄带载波滤出来再做包络检波最后对包络做FFT。包络谱上会出现轴频和叶频的谱线转速和叶片数从而可以直接读出。它的核心假设是“噪声的幅度调制信息集中在低频”这个假设在大多数螺旋桨噪声场景成立但在多途干涉严重或目标远离时可能被破坏这一点留到避坑章节专门讲。2.2 仿真信号模型不是伪正弦而是“窄带噪声乘周期包络”仿真的第一课是建模方式。很多初学者生成仿真信号时直接写一堆正弦函数x(t)A1·sin(2πf_s·t)A2·sin(2πZf_s·t)noise。这种做法做频谱分析没毛病但做DEMON谱分析会给你错误的信心因为真实螺旋桨信号是“被周期调制的宽带噪声”不是“被噪声污染的周期信号”。正确模型如下x(t) [1 m_s·cos(2πf_s·t) m_b·cos(2πZf_s·t)] · n_bp(t) v(t)其中n_bp(t)是经过带通滤波的高斯白噪声代表空化噪声的窄带分量v(t)是加性环境噪声。m_s和m_b是调制深度取值通常在0.3到0.8之间。若调制深度接近1包络会出现过调实际检波后谐波分量变得非常丰富反而干扰基频识别若小于0.3谱峰就不够突出容易被噪声底淹没。真实数据里调制深度还会随航速、深度变化仿真阶段把m_s取0.5、m_b取0.4是比较稳的起点。为什么要坚持“窄带噪声乘周期包络”这个结构原因有两点其一真实信号的载波部分具有随机相位和随机起伏用正弦载波会低估滤波器群延迟对包络的影响其二Hilbert变换对窄带噪声做包络检波时瞬时幅度本身就带随机波动这正好用来检验后级低通平滑参数是否合理。仿真不是越漂亮越好而是要保留真实场景里那些“让算法不舒服”的随机性。下面这张表给出了后续代码统一使用的仿真参数方便你在跑通后逐项改动看效果。参数仿真设定值物理含义轴频 f_s2.5 Hz螺旋桨每秒转2.5圈对应150 rpm叶片数 Z6常见商船螺旋桨叶片数叶频 Z·f_s15 Hz每片桨叶周期性扫过形成的调制频率载波频段600–1000 Hz空化噪声中能量较集中的一段窄带调制深度 m_s / m_b0.5 / 0.4幅度调制的强弱实测一般在0.3–0.8信噪比约20 dB加性噪声相对调制信号的强度2.3 仿真在整个分析流程中的位置先有标准答案才能谈误差DEMON谱分析流程可以拆成四步带通滤波、包络检波、低通平滑、FFT谱分析。每一步都有参数而真实录音里没有标准答案你只能凭经验判断谱峰是否可信。仿真的价值在于它是带标准答案的轴频、叶频、调制深度都由你设定跑完流程后把谱峰位置与设定值对比误差直接量化。这也是demon.zip这类资料的价值所在——它通常把可运行的仿真脚本、样本信号和说明文档聚在一起让你先在干净环境里调通流程再换真实数据。但注意仿真只能标定算法流程和参数范围不能替代实测验证。它不含多途干涉、海洋环境噪声的非平稳变化、传感器与船体振动耦合这些因素会让实测包络谱比仿真脏得多。所以正确的使用方式是先在仿真信号上把流程跑稳记录每种参数下的谱图形态再拿这些经验去解读实测数据。下一章的代码就是在仿真信号上把完整流程走一遍并给出可重复的输出和判断标准。3. 用Python从零跑通DEMON谱分析仿真信号、包络检波与谱线输出3.1 生成仿真信号轴频2.5Hz、叶频15Hz、载波800Hzimport numpy as np from scipy.signal import butter, sosfilt fs 48000 # 采样率覆盖600-1000Hz载波留足滤波余量 T 30.0 # 信号时长保证能切出多帧做验证 t np.arange(0, T, 1/fs) shaft_freq 2.5 # 轴频对应螺旋桨转速150 rpm blade_num 6 # 叶片数 blade_freq shaft_freq * blade_num # 叶频 15 Hz mod_depth_s 0.5 # 轴频调制深度 mod_depth_b 0.4 # 叶频调制深度 # 生成宽带噪声再带通成窄带载波 rng np.random.default_rng(42) noise rng.standard_normal(len(t)) sos_bp butter(4, [600, 1000], btypebandpass, fsfs, outputsos) carrier sosfilt(sos_bp, noise) # 窄带噪声载波 env_mod 1 mod_depth_s * np.cos(2*np.pi*shaft_freq*t) \ mod_depth_b * np.cos(2*np.pi*blade_freq*t) signal carrier * env_mod # 载波被周期包络调制 signal signal 0.1 * rng.standard_normal(len(t)) # 加性底噪第16到17行把白噪声带通到600–1000Hz得到窄带载波第19到20行构造周期包络再用它去调制载波。注意顺序和物理过程的对应先有载波再有包络相乘最后才加环境噪声。调制深度都小于1包络不会过零检波后不会出现严重谐波。随机种子固定为42方便你调参数时对比前后谱图变化。采样率为什么取48000Hz载波中心频率800Hz理论上16kHz采样率就够但带通滤波器的过渡带需要空间采样率太低会让滤波器阶数变得很高相位失真加重。取48kHz还有个好处即便以后想换用更高频段的实测数据例如几千赫兹的流噪声采样率仍然够用不用重写代码。信号时长30秒是为了后续能切出5段6秒数据做多帧验证这也是实测中比较常见的录音长度。3.2 带通滤波与包络检波Hilbert变换的正确姿势from scipy.signal import hilbert, butter, sosfilt # 第一步接收端再滤一次模拟带通处理 sos_rx butter(4, [500, 1200], btypebandpass, fsfs, outputsos) x_bp sosfilt(sos_rx, signal) # 第二步Hilbert变换取包络 analytic hilbert(x_bp) envelope np.abs(analytic) # 第三步去直流并低通平滑 envelope envelope - np.mean(envelope) sos_lp butter(6, 200, btypelow, fsfs, outputsos) env_smooth sosfilt(sos_lp, envelope)接收端的带通带宽故意比发射端宽设成500–1200Hz而不是完全匹配的600–1000Hz。这是模拟滤波器失配实际接收系统的通带不可能和发射信号完全对齐留出余量才能看出算法对失配的容忍度。如果接收端带通和发射端完全一致谱图会过于理想实测时反而容易翻车。Hilbert变换作用于x_bp再用np.abs取模得到瞬时幅度就是包络。这里有个容易犯的错直接对原始宽频信号做Hilbert。那样做包络里会混入载波频段的残余波动后级低通如果截止频率设得不够低包络谱里会多出一堆高频假峰。必须先带通、再Hilbert、再低通顺序不能乱。最后低通截止200Hz是为了保留0–200Hz的调制谱内容同时把载波频段的波动压下去。3.3 对包络做FFT窗函数、谱线与结果解读N 16384 win np.hanning(N) start int(fs * 20) # 跳过滤波器和Hilbert的起始瞬态 segment env_smooth[start:startN] * win spec np.fft.rfft(segment) freqs np.fft.rfftfreq(N, 1/fs) power np.abs(spec) ** 2 # 只看0-30Hz的调制谱区域 lo, hi 0, 30 mask (freqs lo) (freqs hi) for f, p in zip(freqs[mask], power[mask]): if p 1e-3: # 只打印明显谱峰 print(f{f:.3f} Hz power{p:.2e})FFT点数用16384对应的频率分辨率是fs/N 48000/16384 ≈ 2.93Hz。轴频2.5Hz在2.93Hz分辨率下会偏到2.93Hz附近叶频15Hz则相对准一些。所以只看原始FFT峰值位置轴频肯定是有偏的后面需要用抛物线插值修正。窗函数选Hanning主瓣宽约四个频率分辨率旁瓣衰减较快对找峰更稳。start从20秒处开始切是因为滤波器和Hilbert变换在起始段有瞬态响应直接取前几秒数据会把边界失真带进谱里。怎么判断跑出来的谱对不对打印结果里应该能看到15Hz附近的强峰以及2.5Hz附近某个偏了几分之一赫兹的峰。如果只有15Hz的峰而没有轴频峰说明轴频调制深度太低或搜索区间太窄如果0–3Hz内出现一片连续高能量说明低通截止设得太高把载波残余放进了包络谱。这些都是后文避坑章节的常见问题先有个印象。3.4 关键参数速查表及其中的理解要点参数取值选择理由采样率48000 Hz覆盖载波频段留足滤波器过渡带发射端带通600–1000 Hz模拟空化噪声中的能量集中窄带接收端带通500–1200 Hz故意放宽测试滤波器失配容忍度信号时长30 s可切5段6秒做多帧验证FFT点数16384分辨率约2.93Hz接近但不足需插值包络低通截止200 Hz调制谱能量集中在几十赫兹内留足余量窗函数Hanning主瓣宽度可控旁瓣衰减适中提示FFT分辨率是本流程最容易被忽视的坑。轴频只有2.5Hz分辨率却只有2.93Hz谱峰必然偏移。下一章专门讲怎么把这个误差降下来以及滤波器上下限、低通截止、谱峰判定各自该怎么定。4. 参数决定成败带通上下限、包络平滑与谱分辨率怎么调4.1 带通滤波上下限按载波频段选不按轴频选带通滤波的目标是“隔离出一段信噪比较高的载波窄带”这段窄带的中心频率和带宽由实际接收频段决定和轴频没有直接关系。但初学者最常见的误用是看到轴频2.5Hz就把带通设成0到10Hz结果滤完只剩低频环境噪声包络谱里全是慢起伏的假峰轴频特征反而被淹没。选带通范围时先看接收信号的总功率谱找到能量比较集中的凸起频段。螺旋桨空化噪声的能量通常比较平坦但不同航速和深度下会有些偏移常见做法是在预估频段内用扫频方式找到包络谱峰最清晰的一段。带宽方面一般取中心频率的20%到40%作为起步值中心频率800Hz时带宽取160到320Hz对应上例中的600–1000Hz带宽400Hz算是偏宽的设置已经可以工作如果实测谱峰糊优先把带宽收窄到700–900Hz试试。带宽太窄也有问题窄带内噪声功率低Hilbert检波后包络的随机起伏变大相当于信噪比下降。带宽太宽则可能把两个调制机制不同的频段混在一起包络谱互相叠加峰就糊了。所以带通参数要看包络谱结果迭代不要一次定死。我一般会写一个简单的参数扫描脚本把中心频率从500Hz到2000Hz、带宽从10%到50%各跑一遍对比轴频峰高和噪声底再选最优组合。4.2 包络检波后的低通截止别把几十赫兹的调制谱滤没了包络信号仍然带着载波频段的波动必须低通。但低通截止频率选多少直接决定你能不能看到叶频。假设轴频2.5Hz、叶片数6叶频15Hz加上二次谐波30Hz低通截止至少放到40到60Hz。上面代码里用200Hz是为了稳妥但如果信号里存在强低频晃动比如海面波浪起伏200Hz会把不少干扰放进包络谱。经验做法是先看包络信号的功率谱密度确认调制能量集中在哪个频段再定截止频率不要一上来就定死。观察点的原则是低通截止应至少是叶频的三倍给谐波留空间但不要超过载波频段的十分之一否则载波残余压不干净。以15Hz叶频为例截止取45到60Hz比较合适如果后续发现轴频谐波丰富再放宽到80Hz。低通阶数同样影响结果。阶数越高阻带衰减越快但群延迟越大谱峰位置会偏移。6阶Butterworth是折中群延迟在200Hz截止时大约几十毫秒对2.5Hz轴频来说相位影响很小。如果使用FIR滤波器要注意群延迟是常数可以在时间轴上统一补偿比IIR更好处理但计算量更大。4.3 FFT点数、窗函数与频率分辨率轴频2.5Hz需要多长的数据频率分辨率公式是Δf fs / N。要用FFT分辨2.5Hz和5Hz两个相邻谱峰轴频和它的二次谐波Δf至少要小于两个频率间距的一半也就是1.25Hz对应N ≥ 38400在48kHz采样率下约0.8秒数据。但实际谱峰旁边有旁瓣还有噪声底把Δf压到目标频率的十分之一以下才安心也就是N ≥ 192000约4秒数据。如果信号足够长优先增加参与FFT的数据长度而不是补零。补零只能让频谱曲线看起来更平滑并不能把两个本就叠在一起的谱峰分开。很多人在短数据上补零到几百万点以为分辨率提高了其实真实分辨能力一点没变只是“插值显示”的假象。正确做法是把N取到信号实际长度的80%以上剩下的部分用窗函数平滑截断。窗函数的影响也要算进去。Hanning窗的主瓣宽度约4Δf这意味着即便Δf已经压到1Hz两个相距1.5Hz的谱峰也会因为主瓣重叠而难以分离。这时候要用更窄主瓣的窗比如Blackman-Harris但它旁瓣衰减更慢需要更长数据来补偿。实际工程里我优先保证数据量再选窗最后才谈分辨率。4.4 谱峰判定与信噪比阈值为什么固定阈值是玄学很多同学直接在包络谱里找全局最大值然后认为它是轴频。这在干净仿真里没问题在实测里经常翻车因为叶频调制深度常常比轴频大全局最大峰往往是叶频而不是轴频。正确做法是先用先验限制搜索频段船舶轴频通常2到20Hz先在这个区间里找候选峰再把峰高和该区间局部噪声底作比较。局部信噪比用这个公式算SNR_peak 10·log10(P_peak / median(P_local))。P_local取峰值周围±1Hz邻域内的功率值用中位数估计噪声底比全局均值稳健得多。峰值超过噪声底6dB是最低门槛10dB以上才比较可信。为什么说固定阈值是玄学因为不同航速、不同距离下包络谱的绝对功率变化极大一个在A场景合理的阈值在B场景可能把全部峰都滤掉或者全部放进来。只有相对局部噪声底的比值才有跨场景意义。航速变化还会让轴频缓慢漂移。所以连续几帧的估计值应该平滑变化单帧突跳的峰多半是干扰。自动处理系统里我一般把上一帧的轴频估计作为当前帧搜索范围的中心只在±0.3Hz范围内搜索这样既排除远处假峰又能跟踪转速缓慢变化。这个思路比任何阈值都有效。5. 避坑DEMON谱分析常见的翻车现场与排查思路5.1 轴频谱线翻倍半波检波惹的祸现象仿真里设定轴频2.5Hz输出谱峰却在5Hz叶频也变成原来的两倍15Hz变成了30Hz。原因早期的模拟DEMON谱分析常用二极管检波本质是半波整流负半周被切掉包络变成“整流后的正弦”基频因此翻倍。数字实现里如果对包络再做了一次绝对值后取低通或者用了非相干的平方律检波但没有开方也会出现同样的翻倍效果。解决检查包络提取链路。用Hilbert变换后取模是标准做法确认没有额外再取一次绝对值。观察原始包络波形如果负半周完全消失波形有明显的“半正弦”形状就可以判定翻倍来源。这个现象也可以反向利用如果轴频谱线消失只有叶频用叶片数先验去除后可以反推轴频。5.2 谱线糊成一团带通窗口太宽或低通截止太高现象0到30Hz范围内没有明显的独立谱峰而是一连串密集小峰噪声底整体抬高看起来像一片丘陵。原因接收端带通范围太宽包络里混入多段不同调制特征的载波各段包络互相叠加或者包络低通截止设到500Hz以上把载波频段的残余波动一起放进了包络谱。解决把带通收窄到载波能量最集中的区间例如从500–1200Hz改成700–900Hz同时把包络低通截止降到60Hz左右再试。判断标准是包络信号的频谱里载波频段成分应该比调制频段成分低20dB以上。如果改完仍然糊先怀疑噪声底太高而不是算法问题把信号时长拉长到60秒以上再做多帧平均。5.3 谱峰算出来总漂移FFT分辨率不够时先插值再怀疑算法现象多次试验中轴频估计值在2.5Hz附近来回跳2.4、2.7、2.9都有但设定值一直没变。原因FFT频率分辨率只有2.93Hz真实频率落在两条谱线之间时峰值位置取决于窗函数的形状和噪声的随机起伏估计结果直接在量化步长上跳变。这不是算法错了是分辨率不够。解决先用抛物线插值把峰值频率修正到亚分辨率级别具体代码在下一章。如果插值后仍然漂就增大N或用Zoom FFT把目标频段单独细化。不要用补零冒充分辨率补零只改变显示密度不改变谱线之间的真实分辨能力。5.4 仿真能出谱实测却找不到峰多途效应与海洋环境干扰现象同样的参数仿真信号包络谱里轴频清清楚楚换成实测录音后峰消失或位置乱跳。原因海洋信道多途传播导致同一信号的多个时延副本叠加包络的周期性被破坏海洋环境噪声风浪、生物、航船是非平稳的可能在某几秒内把信号完全盖住。解决把单次FFT改成多帧平均例如取10段4秒数据分别做包络谱后中值平均同时记录每帧的局部信噪比低于6dB的帧直接丢弃不要硬猜。还可以用窄带跟踪在上一次估计轴频附近的±0.3Hz区间内滑动搜索峰值而不是每次全频段扫。多途造成的包络相位畸变很难完全消除但通过多帧中值平均可以把随机部分压下去。5.5 Hilbert变换边界跳变包络检波的边缘失真现象包络信号开头几十毫秒有一个很大的尖峰导致第一帧FFT里出现异常高的低频分量轴频峰反而被淹没。原因Hilbert变换是全局运算信号截断边缘会产生边界效应带通滤波器在起始段也有暂态响应两者叠加让包络前端失真。这个问题在做逐帧处理时特别明显因为每一帧的起点都是截断点。解决生成仿真信号时多留一段长度分析时把前2到3秒丢弃对实测数据每帧开头丢弃滤波器群延迟加暂态时间的总和通常取0.5秒以上。另一个办法是先把整段长信号做Hilbert变换再分段做后续处理这样每帧只有一次边界效应而不是每帧都有。5.6 叶频被当成轴频先验信息缺失的工程后果现象输出谱峰集中在15Hz看起来干净漂亮但实测的轴频其实只有2.5Hz转速估算直接差了6倍后续目标识别全部错位。原因包络谱里叶频的调制深度往往比轴频大全局找最大峰时容易锁定叶频。如果没有轴频范围的先验约束算法会把叶频当作主轴频而且谱峰越干净越难察觉。解决用轴频先验范围限定搜索区间例如2到20Hz在识别到候选峰后检查它是否存在约整数倍关系、周围是否还有更低的谱峰。自动系统里可以融合叶片数知识同一目标连续帧的轴频应平滑变化而叶频可能随负载变化跳变。另外如果转速估出来是整数赫兹或特别规整的数值也要多留个心眼真实轴频极少恰好落在整数上。6. 把DEMON谱做成自动工具谱峰插值、转速估计与结果验证6.1 用抛物线插值把轴频估计精度提高一个量级idx np.argmax(power[mask]) idx_abs np.where(mask)[0][idx] f_bin freqs[idx_abs] p_left, p_center, p_right power[idx_abs-1], power[idx_abs], power[idx_abs1] delta 0.5 * (p_left - p_right) / (p_left - 2*p_center p_right) f_est f_bin delta * (fs / N) rpm f_est * 60 print(f轴频估计: {f_est:.3f} Hz, 转速约 {rpm:.1f} rpm)抛物线插值的前提是峰值两侧的功率值满足p_center大于左右两边如果峰值正好落在第一条谱线上需要给左侧补一个很小的量再算。Hanning窗下插值精度可以做到优于十分之一的频率分辨率也就是0.3Hz以内。但要注意插值只修正周期图峰值的位置不会把两个本就叠在一起的谱峰分开所以它不能替代增加数据长度。6.2 多帧中值平均与置信度检查怎么判断结果可信把30秒信号切成5段6秒每段独立做一次完整的带通、包络、FFT流程得到5个轴频估计。取中值作为最终结果同时看5个值之间的极差极差小于0.3Hz说明整个过程稳定可信极差大于0.5Hz就说明数据里存在非平稳干扰需要重查滤波参数或丢帧。另一条检查是转速合理性水面船常见轴速60到300rpm1到5Hz水下航行器更低算出一个60Hz的“轴频”就要报警那多半是把叶频或电路工频干扰当成了目标。我自己的习惯是每次跑完DEMON谱分析都把原始信号段、带通后波形、包络波形和最终包络谱存成PNG和npy方便事后回查。早年吃过没有留中间结果的亏调完参数以为效果变好了换了数据却发现之前能出峰的谱图现在出不来又舍不得删新参数最后只能靠git提交记录找回旧版本。所以如果你打算照着上面的流程长期用建议从一开始就把中间结果存下来否则后面排查“哪个环节变了”会花掉大量时间。希望帮到你。本文还有配套的精品资源点击获取