简介双谱分析是一种高级信号处理技术擅长从非线性、非高斯信号中提取隐藏特征广泛应用于机械故障诊断与设备状态监测。压缩包内提供了一套基于 MATLAB 的双谱分析实现包含 2 个 .m 文件包体仅 3KB代码精简、依赖少便于快速阅读和二次开发。两个脚本分别对应双谱计算与西储大学故障数据分析前者完成从傅里叶变换到互功率谱、双谱结果输出的核心链路后者针对具体故障数据提供应用示例并可通过替换输入数据适应更多故障类型帮助理解双谱图所呈现的频率间非线性依赖及其与设备异常的关系。对于正在学习高阶谱分析、或希望在滚动轴承故障诊断中引入双谱方法的工程师这套资源提供了可直接运行的参考实现。目前已有 2645 人学习代码结构紧凑适合教学自学也可为工程应用提供基础是从理论到实战的快捷起点。1. 先把“双谱”这件事说人话做信号处理这些年我越来越觉得光盯着功率谱看远远不够。你手上有一段信号如果用功率谱分析本质上是把信号拆成不同频率成分看每个频率上“有多少能量”。这个思路在很多场景下够用比如看个频谱峰值、做个滤波设计、估计个信噪比。但一旦遇到非平稳、非线性或者带有相位耦合的信号功率谱就非常容易翻车——因为它把相位信息全部丢掉了。双谱分析恰恰补上了这个短板。它是高阶谱分析里最常用的一种本质上是把信号的三阶累积量做傅里叶变换得到的结果既能反映频率成分又能反映频率之间的相位关系。说直白点功率谱告诉你“哪些频率有能量”双谱告诉你“这些频率之间是不是存在某种调制或耦合关系”。这对于机械故障诊断、水声信号检测、生物电信号研究这类需要捕捉非线性特征的场景价值非常大。这篇文章我会围绕“对信号作双谱分析”这条主线从数学原理讲到可落地的 Python 实现再搭配实际案例和踩坑经验。无论你是刚接触高阶谱分析的入门者还是已经在工程里被非线性问题折腾得头疼的从业者都可以照着这篇文章把双谱分析的流程完整搭起来。2. 双谱的数学定义与核心物理含义2.1 从三阶累积量到双谱定义并没有那么吓人如果你去看早期的文献双谱的定义通常写成这样B(f1, f2) ∑∑ C3(τ1, τ2) · e^{-j2π(f1τ1 f2τ2)}其中 C3(τ1, τ2) 是信号的三阶累积量。看着有点抽象但拆开理解就简单了。一阶统计量是均值二阶统计量是方差和自相关刻画的是“信号的波动幅度”三阶统计量就是三阶累积量刻画的是“信号偏离正态分布的程度”。为什么三阶累积量能刻画出非线性我习惯用一个类比你把一个线性系统想象成一张平整的桌子输入正弦波输出还是同频率的正弦波只是幅度和相位变了。但如果系统里面存在非线性比如一个轻微打滑的轴承输入一个频率成分输出里就会出现谐波、和频、差频。这些新产生的频率分量之间不是独立存在的它们有固定的相位关系。功率谱把这种关系当作噪声处理掉了而三阶累积量能保留这种关系双谱则把这种关系呈现在频率平面上。2.2 双谱的对称性不用算全平面双谱一个非常实用的性质是它的对称性。因为它是由三阶累积量做二维傅里叶变换得到的所以它在 (f1, f2) 平面上具有多种对称关系。实际计算时只需要计算一个三角形区域就能覆盖全部信息这个区域通常称为“主区域”或“非冗余区域”。以采样率 fs 归一化后的频率为例双谱的计算范围通常限定在f2 0f1 f2f1 f2 fs / 2这个三角形的边界条件本质上是跟采样定理和双谱自身的周期性质有关的。你可能觉得这是数学细节但真正写代码的时候正确地处理这个区域能省下接近 83% 的计算量。很多人第一次算双谱吭哧吭哧把整个二维平面都算一遍结果画出来的图一大半都是没用的重复信息后面我会详细说怎么规避。2.3 为什么功率谱做不到相位信息到底有多重要我再用一个更直白的例子说明。假设你有一个信号由三个频率分量组成f1、f2以及 f1f2。如果这三个分量的初始相位满足某种约束关系比如 φ3 φ1 φ2这在信号处理里叫二次相位耦合。这时候双谱在 (f1, f2) 处会出现一个明显的峰值说明系统内部存在非线性相互作用。但如果你只看功率谱你只会看到 f1、f2、f1f2 三个频率处有峰完全不知道它们之间还有这层关系。换句话说二次相位耦合这个现象在功率谱里跟三个独立的频率成分看起来一模一样只有双谱能区分。这就是双谱分析的核心价值所在它在频域里保留了相位信息让“频率之间存在关联”这件事变得可见、可测量。3. 双谱估计的完整实操流程3.1 经典方法分段平均与周期图法实际工程里我们不可能拿到理论上的三阶累积量只能用有限长的观测数据去估计。最常用的方法是“分段平均周期图法”思路跟功率谱估计里的 Welch 法很像只是从二阶推广到了三阶。具体步骤分四步把长度为 N 的信号分成 K 段每段 M 个点段与段之间可以有一定重叠。对每段数据去均值必要时加窗函数然后做傅里叶变换得到每段的频谱 X_k(f)。对每段频谱计算三重乘积 X_k(f1) · X_k(f2) · X_k*(f1f2)。对 K 段的三重乘积取平均得到双谱估计。写成代码就是import numpy as np def bispectrum_estimation(x, nfft256, windowhann, overlap0.5): 分段平均法估计双谱 x: 输入信号 nfft: FFT点数 overlap: 段间重叠率 # 分段参数 seg_len nfft step int(seg_len * (1 - overlap)) n_seg (len(x) - seg_len) // step 1 # 加窗 if window hann: win np.hanning(seg_len) else: win np.ones(seg_len) # 用于累加三重乘积 bispec np.zeros((nfft, nfft), dtypecomplex) for i in range(n_seg): seg x[i*step : i*step seg_len] seg seg - np.mean(seg) # 去均值 seg seg * win # 做FFT X np.fft.fft(seg, nfft) # 三重乘积并累加 for f1 in range(nfft//4, nfft//2): for f2 in range(nfft//4, f11): f3 f1 f2 if f3 nfft: bispec[f1, f2] X[f1] * X[f2] * np.conj(X[f3]) # 取平均 bispec / n_seg return bispec这段代码虽然能用但双层循环在 Python 里跑起来非常慢。实际项目里我一般会先用 numpy 的广播特性或者干脆用矩阵运算来加速更高效的做法后面“计算量优化”小节会讲。3.2 参数怎么选四个关键决定成败第一个是分段长度。分段长度决定了频率分辨率。假设采样率是 1000 HzFFT 点数是 256那么频率分辨率大约是 3.9 Hz。如果你的信号里有间隔很窄的谱峰分辨率不够就分不开。但分段越长分段数越少双谱估计的方差会变大这就是一个典型的偏差-方差权衡。我的经验是先根据目标频率分辨率倒推 FFT 点数再根据总数据长度确认分段数是否能达到足够平均次数。第二个是重叠率。重叠率太低分段数少估计不稳定重叠率太高段与段之间相关性增大虽然分段数多了但有效自由度并不会等比例提升。一般取 50% 到 75% 之间比较合理。第三个是窗函数。加窗是为了抑制频谱泄漏。汉宁窗是默认选择但如果你的信号里存在强线谱成分可以考虑用凯塞窗旁瓣衰减更好控制。需要特别注意窗函数会引入幅度修正双谱计算关注的是相对关系而非绝对幅度所以一般不用做幅度修正但如果你后续要比较不同信号的绝对双谱值就得把窗函数的能量归一化考虑进去。第四个是去均值。很多人忽略这一步其实非常关键。直流分量会对低频区域的双谱估计产生严重影响。如果信号里有缓慢变化的趋势项最好先做一次高通滤波或者去趋势处理再进入双谱分析流程。3.3 基于 HOSA 工具箱验证结果自己写代码容易出错尤其是索引映射这种细节。我建议在正式分析之前先用现成的工具箱验证一下你的实现是否正确。MATLAB 里有高阶谱分析工具箱HOSAPython 里虽然没有官方标准库但有一些第三方实现可以参考。验证方法很简单构造一个已知具有二次相位耦合的信号比如fs 1000 t np.arange(0, 1, 1/fs) # f150Hz, f2120Hz, 且存在相位耦合 x np.sin(2*np.pi*50*t) np.sin(2*np.pi*120*t) np.sin(2*np.pi*170*t 0.5)然后跑你的双谱估计代码看 (50, 120) 这个坐标处是否出现明显的峰。如果峰的位置对不上优先检查频率轴映射关系再看三重乘积里第三个频率是不是 f1f2。用这个办法我每次在新的环境里重新实现双谱算法都能在十分钟内确认正确性强烈建议你也这么做。4. 双谱结果怎么读峰值、切片与相位判据4.1 双谱幅值图看图说话的核心方法双谱分析做完最直接的可视化方式是画双谱幅值等高线图或伪彩色图。坐标轴是 f1 和 f2颜色代表双谱幅值。由于前面提到的对称性你只需要关注主三角区域内的内容。在幅值图里峰值出现在 (f1, f2) 处意味着该频率对存在明显的非线性耦合。这里有个容易犯的错误峰值出现的位置要结合 f1f2 是否也在频谱上有能量来判断。如果 f1f2 处根本没有信号成分那这个峰值很可能是估计误差造成的不能当作有效耦合证据。我一般会同时绘制三张图双谱幅值等高线图、沿固定 f2 切片的幅值曲线、以及对数尺度下的幅值图。等高线图适合看整体分布切片适合精确定位峰值频率对数图适合观察动态范围大的弱耦合成分。4.2 双谱相位另一个信息维度很多人算完双谱只看幅值忽略了相位信息。双谱的相位本身也有物理意义在某些场景下甚至比幅值更有用。例如在分析脑电信号时不同脑区之间的相位耦合模式可能与认知状态密切相关在机械故障诊断中双谱相位的稳定性也可以作为故障特征。具体判读时如果双谱相位在某个频率对附近保持近似恒定说明该频率对之间存在稳定的相位关系如果相位杂乱无章说明即使幅值上有峰耦合也可能是随机的。一个实用的量化指标是“双谱相干系数”bicoherence它把双谱幅值做了归一化取值在 0 到 1 之间越接近 1 说明耦合越强。这个指标比直接用双谱幅值更稳健推荐在工程分析中优先使用。4.3 二次相位耦合的定量判据判断是否存在二次相位耦合最直接的是看双谱相干系数。它的定义是bic(f1, f2) |B(f1, f2)| / sqrt(P(f1) · P(f2) · P(f1f2))其中 P(f) 是功率谱。这个公式很好理解分母相当于把三个频率成分各自的能量做了归一化分子是三个频率成分协同作用的强度。如果三个频率成分只是恰好同时存在但彼此独立分子会比较小bic 趋近于 0如果存在真正的相位耦合bic 会明显偏大。实际使用中我会设置一个阈值比如 bic 0.3同时对应频率对在幅值图上也有峰值才判定为存在二次相位耦合。这个阈值不是绝对的要根据信噪比和数据量调整但至少给了我们一个客观的判据而不是光靠肉眼看图。5. 双谱分析在哪些场景下真正有用5.1 旋转机械故障诊断从振动信号里找“非线性指纹”在滚动轴承和齿轮箱的振动信号里早期故障往往表现为非线性特征。比如轴承出现局部剥落时振动信号中会产生冲击成分这些冲击会在频域激发一系列谐波和边频带而且各成分之间往往存在相位耦合关系。这种耦合在正常状态下几乎不存在但在故障状态下会变得非常明显。我做一个实际项目时需要区分正常轴承和轻微磨损轴承。单纯看功率谱两者的差异很小只是某些频段幅值略有抬升很难给出稳定判据。后来改用双谱分析在特定频率对 (fc, fr) 处正常轴承的双谱幅值几乎为零而磨损轴承出现了一个非常明显的峰。这个特征稳定复现最终成了诊断模型的一个关键输入特征。5.2 水声信号检测从强噪声背景里捞弱信号水声环境里背景噪声极其复杂不仅有环境噪声还有各种生物、船只辐射的干扰。如果目标信号和干扰都属于高斯分布那么理论上它们的三阶累积量为零双谱中不会留下痕迹。但实际的水声信号往往是非高斯的双谱就能提供额外的辨识维度。在被动声呐场景中螺旋桨空化噪声会表现出明显的非线性调制特征。双谱分析可以提取这些特征用来区分不同类型的辐射噪声源。这个方向上双谱虽然不能完全替代传统的谱分析但作为辅助特征已经表现出了很强的实用性。5.3 生物电信号研究脑电与心电中的耦合分析脑电信号的非高斯性和非线性早就被广泛证实。双谱分析在研究脑电的节律耦合、不同脑区之间的同步性、以及某些神经系统疾病时提供了很独特的信息。比如在癫痫发作期的脑电信号中特定频段之间的相位耦合模式会发生明显改变双谱可以捕捉到这种变化。心电信号的研究中双谱分析被用来分析心率变异性中的非线性成分。传统的 HRV 频谱分析主要看低频和高频功率的比值但双谱能进一步揭示自主神经调节中的非线性交互作用这对某些疾病的早期诊断有潜在价值。需要注意的是生物电信号的信噪比通常较低而且容易受到肌电伪迹、工频干扰等影响。因此在实际分析前必须做好预处理带通滤波、剔除伪迹段、必要时使用独立成分分析分离干扰。6. 计算效率与代码优化技巧6.1 为什么直接双层循环会慢到怀疑人生如果你用前面那段朴素的 Python 代码处理一段 10 秒、采样率 10 kHz 的信号FFT 点数选 256那么 f1 和 f2 的循环大概有一万多次迭代每次迭代还要做复数乘法和数组索引。实测下来跑一次要好几秒。如果参数搜索、多通道分析或者实时性要求高的场景这个速度完全不能用。优化思路主要是两条一是利用双谱的对称性只计算主三角区域二是用矩阵运算替代 Python 循环。6.2 用矩阵运算改写核心计算以 FFT 点数 N256 为例我们构造一个二维频率索引网格利用 numpy 的广播机制一次性完成三重乘积计算然后只保留主三角区域。这样可以轻松获得一两个数量级的加速。def bispectrum_fast(x, nfft256): # 分段去均值加窗逻辑略 X np.fft.fft(x, nfft) # 构造频率索引网格 f np.arange(nfft) f1, f2 np.meshgrid(f, f, indexingij) f3 f1 f2 # 只计算有效区域 valid (f3 nfft) (f2 0) (f1 f2) # 三重乘积矩阵 bispec np.zeros((nfft, nfft), dtypecomplex) triple X[f1] * X[f2] * np.conj(X[f3]) bispec[valid] triple[valid] return bispec注意这里只演示了单段的计算。实际分段平均时可以提前把每段的 FFT 结果存成一个大矩阵然后在矩阵层面做平均进一步减少开销。6.3 降采样、频段裁剪与矢量化如果信号带宽本身不大可以先做带通滤波再降采样这样同样的数据长度下FFT 点数可以选小一点计算量能大幅度降低。但降采样前一定要做好抗混叠滤波不然高频成分折叠到低频区域会污染双谱结果。另一个技巧是频段裁剪。如果你只关心某个频段内的耦合关系比如机械故障诊断中的 0 到 500 Hz可以把 FFT 结果中不需要的频点直接置零或者只保留感兴趣的频率索引进行计算。虽然理论上不能完全消除泄漏影响但实践中配合窗函数能获得满意的结果。还有一点值得提如果数据真的很长可以考虑把分段平均的循环放到 numba 或者 Cython 里加速。我试过用 numba 把朴素的循环版本加速基本能达到接近 C 语言的性能代码改动也少适合不想折腾矩阵化的人。7. 实操中的常见问题与避坑指南7.1 数据长度不够双谱变成一张“麻子脸”新手最容易遇到的问题是数据太短分段平均次数不足双谱估计结果方差极大画出来全是随机噪声般的尖峰。解决这个问题没有捷径只有两个方向一是尽量获取更长的数据二是降低频率分辨率需求缩短 FFT 点数以增加分段数。一个实用经验是对于双谱分析至少保证 20 段以上的平均次数。如果数据实在有限可以适当增加重叠率到 75%这样在总数据长度不变的情况下能多出一些段数。但要注意重叠率太高段之间的相关性增加实际改善有限要清楚这个边际效应。7.2 噪声和干扰太大把假峰当成真峰双谱的一个特性是对高斯噪声不敏感因为高斯信号的三阶累积量理论为零。但现实中的噪声往往不是理想高斯的比如脉冲噪声、工频干扰以及信号本身泄露出来的强线谱都可能在双谱里制造假峰。我踩过的一个坑是信号的基频分量很强即使加窗后仍有泄漏导致在双谱上出现一条“脊线”差点被误判为故障特征。后来我养成了一个习惯在做双谱分析之前先对信号做一次窄带滤波把已知的强线谱成分抑制掉再进行分析。另外使用双谱相干系数替代双谱幅值作为判据也能有效抑制这类由能量差异带来的假象。7.3 频率轴标错了结果全对坐标全错这是我在实现过程中犯过最隐蔽的错误。FFT 的输出频率轴从 0 到 fs但双谱的定义通常使用归一化频率0 到 0.5 或 0 到 1 之间如果映射关系搞错会导致双谱峰出现在错误的位置上。强烈建议在验证阶段就使用已知二次相位耦合的合成信号把频率轴校准好再处理真实数据。索引映射看起来是小问题实际排错时非常折磨人。7.4 不知道用什么工具先验证再投入如果你只是想快速验证双谱分析对某个问题是否有用不建议一开始就自己造轮子。先在 MATLAB 的 HOSA 工具箱或者 Python 的第三方实现上跑通流程确认特征提取思路可行再考虑针对场景做代码优化和定制。这样可以避免把大量时间花在调试底层实现上。提示我在实际项目中采用的流程是“合成信号验证算法 → 真实数据功能验证 → 特征有效性评估 → 工程化优化”这个顺序基本可以保证每一步都建立在可靠的基础上。8. 个人实操体会与下一步想做的事做了不少双谱分析的工程应用后我最大的体会是双谱分析不是万能的但它确实提供了一个功率谱完全给不到的视角。关键是要想清楚你要解决什么问题——如果只是看频率成分和能量分布功率谱足够了没必要上双谱但你一旦怀疑系统存在非线性耦合或者需要在低信噪比下提取非高斯信号特征双谱就值得纳入工具箱。从实操角度我给后来者三个建议。第一永远先做合成信号验证这是检验实现正确性最有效的途径。第二不要只看双谱幅值图配合双谱相干系数才能得到更可靠的结构性判断。第三双谱分析的计算开销确实不小但通过分段参数调整和矩阵化优化绝大多数场景都可以在可接受的时间内完成。我目前正在尝试把双谱特征跟深度学习模型结合起来用双谱图作为神经网络输入用于更复杂的模式识别任务。这个方向还处在探索阶段但初步结果已经显示出一些有意思的可能性。后续有进展了再单独写一篇分享。本文还有配套的精品资源点击获取