★ FEATURED ARTICLE 简介在复杂海洋环境下水声信号常呈现非平稳、多调制特征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 解调。三者数学形式上等效工程表现差异很大。方法实现复杂度输出特性适用场景绝对值检波最低一行代码输出含直流偏置需去均值快速预览信号较强时可用平方检波低乘法即可输出含直流偏置但调制分量倍频清晰工程上最常用SNR 较低时优于绝对值Hilbert 解调较高需要hilbert()输出解析信号幅度相位信息完整窄带信号、需要瞬时相位时我实际做水声信号处理时默认先试平方检波。原因是平方检波在数学上等价于求瞬时功率对宽带噪声的包络调制特别敏感而且频谱上调制分量的幅度是绝对值检波的两倍线谱更突出。缺点是多了一个直流分量需要减均值另外平方会把高频分量也压进低频区域所以低通滤波一定要跟上。Hilbert 解调在窄带信号上表现更好但舰船辐射噪声是宽带信号解析信号相位这一优势用不上反而要多花计算量。绝对值检波能快速看个大概但谱线毛刺多不适合精确估频。三种方法有一条共同底线检波之前必须先做带通滤波。原因有两点一是滤掉海底/海面低频环境噪声这些噪声会直接进入调制频带干扰轴频峰二是在高频段选一个信噪比更好的窗口不同船的辐射噪声高频衰减率不一样选错了窗口包络里调制深度会不足轴频峰直接被噪声底扛住。这份资源包里的bandp3_10k.mat从命名看就是负责这个任务——3 kHz 到 10 kHz 的带通滤波器。2.3 这份资源「质料完整」在哪里四个文件串成一条闭环拆开 demon.zip里面四个文件刚好对应一条完整的可复现链路demon11_16.mMATLAB 主脚本把下面三个数据文件串起来跑流程。bandp3_10k.mat带通滤波器系数通带 3 kHz10 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 文件清单与分工先把四个文件的角色理清楚后面跑脚本时心里有数文件类型在链路里的角色demon11_16.mMATLAB 脚本主控流程读 wav、调滤波器、检波、FFT、画图bandp3_10k.mat数据文件存放带通滤波器系数通带 3 kHz10 kHzlow300.mat数据文件存放低通滤波器系数截止约 300 Hz实验3 带大船 10.47-10.51.wav数据文件实测舰船辐射噪声用于谱分析拿到手的第一步不是直接双击运行而是用whos看一下两个 mat 文件里存的变量名。我拆过不少这类资源包变量命名五花八门有的存成b和a有的存成num和den还有的存成sos和g二阶段节格式。脚本里 load 之后直接用变量名如果对不上后面的filter或filtfilt直接报错。先检查变量名能避掉四分之一的问题。3.2 核心流程与参数设置我一般会先看主脚本结构确认它走的流程是不是我预期的标准链路。拆开后发现demon11_16.m的核心流程基本是下面这个样子% 读取实测 wavFs 从文件头自动获取 [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 的若干分之一。NFFT4096决定频率分辨率按Fs/NFFT算。如果 Fs 是 44.1 kHz分辨率约 10.8 Hz如果 Fs 是 48 kHz分辨率约 11.7 Hz。这个分辨率对轴频检测来说有点糙——大型商船轴频 15 Hz分辨率必须到 0.5 Hz 以下才靠谱。所以我通常会把NFFT调到 65536 甚至更高配合pwelch自带的分段平均稳定性和分辨率两头兼顾。xlim([0 50])画 050 Hz 是因为螺旋桨轴频不会太高。大型商船螺旋桨转速 60150 rpm对应轴频 12.5 Hz快艇转速上千转轴频也就十几 Hz。把频段卡在 50 Hz 以内谱峰不会被远处的频带干扰。3.3 谱图判读峰在哪轴频就在哪跑完脚本输出的 DEMON 谱图横轴是调制频率纵轴是包络的功率谱密度。判读规则不复杂第一个最突出的峰就是轴频fp后续在2×fp、3×fp位置出现的峰是轴频的倍频说明调制信号非正弦、谐波成分丰富。如果fp附近还有一个相差 24 倍的峰那个可能是叶频fb fp × 叶片数但要注意区分倍频和叶频——倍频是整数倍关系叶频不一定落在整数倍上。对于这段「带大船」实测数据我跑下来的经验是带通窗口选 310 kHz 时商船轴频一般在 1.53 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 处一根冲上天际的谱线110 Hz 区域反而什么都看不见全是「贴地」的噪声底。原因平方检波之后的包络信号里带了显著的直流分量。包络均值不为零直流能量全部集中在 FFT 的零频幅度太大把纵轴的动态范围压扁低频段的真实谱峰被视觉掩盖。解决检波之后、FFT 之前必须做去均值也就是env_lp env_lp - mean(env_lp)。这一步等价于在零频处陷波位置在代码里放在低通滤波之后、pwelch 之前最有效。注意如果先做过一阶高通也能达到类似效果但会引入相位畸变不如直接减均值干净。如果减完均值后零频还是很高说明低通滤波后的包络里还有缓慢漂移分量可以在减均值之前先做一次 detrend把线性趋势也去掉。4.3 带通滤波后信号整体变差包络里几乎所有能量都在高频残留现象低通滤波后的包络还是一串高频锯齿做出来的 DEMON 谱在 050 Hz 区间没有明显峰只有一堆漂移的毛刺。原因典型的滤波器系数不匹配——带通或低通滤波器的实际截止频率和设计意图差很远或者滤波器阶数太低过渡带宽到几百赫兹检波出来的高频残留在低通后没有真正被压掉。解决对bandp3_10k.mat里的滤波器画一次幅频响应确认通带。用freqz(b_bandp, a_bandp, 2048, Fs)看一眼310 kHz 通带内幅度应该平坦10 kHz 以上斜率足够陡。low300 的幅频响应也要确认 300 Hz 以上衰减至少 40 dB。一般这类资源包里的滤波器设计没问题真正的问题在于 Fs 不匹配——如果 wav 的采样率和滤波器设计时的采样率不一致截止频率整体偏移带通可能变成 27 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其谐波分量直接注入水听器链路或通过地环路进入采集系统。这个干扰在带通滤波时没有被滤掉因为它的高频谐波成分可能落在 310 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是全局谱估计如果只在 050 Hz 需要高分辨率用 Chirp-Z 变换做细化能把 12 Hz 分得很清。我自己常用方式是把全局谱先跑一遍锁定范围然后做一次细化分析确认峰位的亚赫兹细节。这些年我跑过的水声信号包不少最深刻的教训就是下载到资源包先别急着看主脚本第一件事是whos检查数据文件里的变量名第二件事是确认 Fs 匹配第三件事才轮到跑流程。这三步走完至少能省掉一半的报错时间。这几年每次拿到新的谱分析资源我都强制自己先走完这三步再动手改参数。希望帮到你祝一次跑通。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
