简介双谱分析是信号处理中的高级技术主要用于检测非线性、非高斯噪声环境下的信号特征在机械设备故障诊断中应用广泛。这份资源以“对信号作双谱分析”为核心主题面向信号处理初学者、算法研究者以及设备维护工程师帮助用户从频域层面揭示信号中隐藏的非线性依赖关系适用于振动、声音等故障信号的离线特征分析。资源以rar压缩包形式提供共包含两个文件均为.m格式的Matlab脚本包体积仅3KB。文件数虽少但浓缩了双谱估计的核心算法例如基于傅里叶变换的互功率谱计算以及针对特定信号的分析程序。使用者既可直接调用函数处理自有数据也能结合注释学习算法细节便于在Matlab环境下快速运行验证。目前已有2652人学习过这份资源属于轻量实用的参考代码。对希望掌握双谱分析计算流程、绘制二维/三维双谱图并尝试将双谱方法应用到旋转机械故障特征提取中的工程师来说可作为入门和移植用例便于缩短代码开发时间提升信号特征识别效率。1. 双谱分析为什么功率谱看不见的东西它能看到做信号处理的人迟早会撞上这类场景。某段实测振动信号功率谱上只有一片鼓包怎么看都分不出成分换成双谱分析图上却出现一个干净利落的峰位置恰好落在某两个频率之和上。这不是玄学——双谱分析拿信号的三阶累积量做二维傅里叶变换保留普通功率谱丢掉的相位信息专门用来发现非线性系统的二次相位耦合而且在理论上对高斯噪声免疫。适合谁轴承齿轮故障诊断、心电脑电分析、语音声学、通信信号识别凡是手里信号可能经过非线性系统、又泡在噪声里的都值得试。下面的内容按“数学原理 → Python 实现 → 参数调节 → 避坑 → 实战技巧”展开每一步都给你能直接改跑的代码和判断依据。2. 双谱的数学原理与三条估计路径从三阶累积量到二维傅里叶变换2.1 三阶累积量为什么能“看见”相位耦合普通功率谱是信号自相关的傅里叶变换本质上是二阶统计量它只回答“哪些频率有能量”不回答“这些频率之间有没有关系”。双谱换了个思路从三阶累积量出发C3(τ1, τ2) E{ x(t) · x(tτ1) · x(tτ2) }对 C3 做二维傅里叶变换就得到双谱 B(f1, f2)。这个定义里藏着两个关键性质。第一高斯过程的三阶矩恒为零所以理论上高斯噪声在双谱里完全不贡献——这是双谱抗噪的立身之本但也有限定条件有限样本下仍有残差且非高斯噪声会破坏这个性质。第二三阶累积量包含三个时刻的乘积频域对应的是三个频率分量的乘积信号相位信息在这里没有被抵消。常见例子是二次相位耦合QPC。设信号由三个正弦组成频率分别是 f1、f2、f3f1f2且相位满足 φ3 φ1 φ2。功率谱只在 f1、f2、f3 处各有一条谱线完全无法告诉你 f3 是独立源还是非线性耦合的结果。而双谱在 (f1, f2) 处的值正比于 E{ X(f1)·X(f2)·X*(f1f2) }相位关系固定时这个期望非零当 φ3 与 φ1φ2 无关时期望趋于零。这就是“看得见”的本质双谱是复值谱模值反映耦合强度相位反映耦合相位关系。很多人把 FFT 当黑匣子用到了双谱这里就必须关心相位了。普通功率谱里 X(f)·X*(f) 把相位直接消掉双谱则是三个谱的乘积相位关系成了核心信息。理解到这一层后面调参数、判读峰值才有依据。2.2 三种估计方法对比直接法、间接法与参数化法双谱实现路径有三条工程上最常见的做法是直接法。直接法先把信号分段对每段做 FFT然后对 X(f1)·X(f2)·X*(f1f2) 做平均间接法则先估计三阶累积量再对累积量做二维 FFT参数化法用 AR 模型系数构造双谱短数据下表现更好但模型阶数敏感。方法核心思路优点缺点直接法分段 FFT 后频域三元乘积平均实现简单、与现有 FFT 流程衔接好、易向量化分段不足时方差大间接法先估三阶累积量再做二维 FFT便于加窗平滑、理论上更直观计算量大延迟范围选择影响结果AR 参数化法用 AR 模型系数构造双谱短数据、低信噪比下表现较好阶数选择敏感偏差不易控制我一般优先用直接法。原因是它直接复用你手里已经调好的 FFT 和分段流程计算可以完全向量化几秒跑完一段 16 秒的实测信号。间接法适合做学术验证比如想仔细控制累积量域的窗函数形状参数化法在数据特别短、比如只有几百个点时值得试但阶数需要多轮实验不适合快速批量分析。2.3 工程上为什么普遍用分段平均直接法直接法如果不分段单段数据的双谱估计方差很大而且这个方差不会随数据点数增加而收敛——这是三阶统计量的固有特性。分段平均就是 Welch 法在高阶谱里的延伸每段独立估计然后对所有段取平均。段数越多估计越一致代价是每段变短、频率分辨率下降。分段平均还有一个实际好处可以直接套用已有的分帧参数。比如你处理音频或振动信号时已经在用 50% 重叠的汉宁窗双谱分析可以直接沿用不需要额外摸索一套流程。我做过的一次诊断项目里原始数据经过带通滤波后先跑一版功率谱定位主要频带再在感兴趣的频段内做双谱整个流程只在原来的分帧参数上改了段长十分钟就出图。需要提醒的是直接法的频域乘积 X(f1)·X(f2)·X*(f1f2) 包含了和频关系 f3 f1 f2。实际信号里如果存在差频非线性比如 f3 f2 - f1对应的峰出现在双谱的另一个象限。所以判读时不能只盯着正频率区域看否则容易漏掉真正的耦合关系。3. 用 Python 从零实现双谱估计可直接改跑的代码3.1 构造一个带二次相位耦合的测试信号验证算法之前先造一个已知答案的信号。采样率 200 Hz时长 16 秒包含 20 Hz 与 30 Hz 两个基频一个 50 Hz 的二次相位耦合项相位取 0等于 20 Hz 与 30 Hz 相位之和再加入一个相位随机的 76 Hz 独立分量和标准差 0.6 的高斯噪声。import numpy as np rng np.random.default_rng(42) fs 200.0 N_total int(fs * 16) # 16 秒数据 t np.arange(N_total) / fs # 二次相位耦合50 Hz 的相位 20 Hz 相位 30 Hz 相位 0 sig ( np.cos(2 * np.pi * 20 * t) np.cos(2 * np.pi * 30 * t) np.cos(2 * np.pi * 50 * t) # 耦合项相位固定为 0 np.cos(2 * np.pi * 76 * t rng.uniform(0, 2 * np.pi)) # 独立项相位随机 ) # 叠加高斯噪声检验双谱对高斯噪声的抑制能力 sig sig 0.6 * rng.standard_normal(N_total)逻辑说明20 Hz 和 30 Hz 是驱动项50 Hz 是它们的和频且相位锁定。76 Hz 分量虽然频率独立但相位是随机的理论上不应在双谱上形成稳定峰。高斯噪声理论上对三阶累积量零贡献。这样一个信号能同时验证“真实峰是否出现”和“伪峰是否被压住”两个问题。参数说明噪声标准差 0.6 相对信号幅值 1.0 来说不算小功率谱上 50 Hz 峰仍然可见但双谱的峰应当比功率谱干净得多。如果你想测试算法极限可以把噪声标准差加大到 1.5 以上观察双谱峰如何被有限样本残差淹没。3.2 核心估计函数分段、加窗、频域三元乘积平均这里给出向量化的直接法实现。分段长度 256 点50% 重叠汉宁窗每段 FFT 之后做幅值归一化再计算所有合法频率组合的三元乘积并平均。def estimate_bispectrum(x, seg_len256, overlap0.5): 直接法双谱估计。 返回B 为 (seg_len//21, seg_len//21) 的复数组 freqs 为对应的频率轴。 if x.ndim ! 1: raise ValueError(输入必须是一维信号) step int(seg_len * (1 - overlap)) n_seg (len(x) - seg_len) // step 1 if n_seg 4: raise ValueError(分段太少建议增加数据长度或减小 seg_len) Nf seg_len // 2 1 win np.hanning(seg_len) # 构造频率索引矩阵k3 k1 k2只保留不超出 Nyquist 的部分 k1 np.arange(Nf)[:, None] k2 np.arange(Nf)[None, :] k3 k1 k2 mask k3 Nf k3_valid np.where(mask, k3, 0) B np.zeros((Nf, Nf), dtypecomplex) for i in range(n_seg): seg x[i * step : i * step seg_len] seg seg - seg.mean() # 去均值消除直流贡献 X np.fft.rfft(seg * win) / win.sum() # 加窗并做幅值归一化 prod X[k1] * X[k2] * np.conj(X[k3_valid]) B prod * mask B / n_seg freqs np.arange(Nf) * (fs / seg_len) # 频率轴 return B, freqs B, freqs estimate_bispectrum(sig, seg_len256, overlap0.5)逻辑说明每段信号先去均值避免直流分量在低频谱区制造伪脊加窗后除以窗能量和保证幅值尺度与矩形窗一致这一步在定量比较时很重要。三元乘积 X(f1)·X(f2)·X*(f1f2) 是双谱的核心只有当 f3 与 f1f2 的相位关系稳定时跨段平均才不会互相抵消。mask 把超出 Nyquist 的频率组合置零本质上双谱有效区域是 f1 ≥ 0、f2 ≥ 0、f1f2 ≤ fs/2 的下三角。参数说明seg_len256 时频率分辨率是 200/256 ≈ 0.78 Hz对 20/30/50 Hz 这种间隔足够overlap0.5 让 16 秒数据产生约 24 段方差控制有余量。如果你的信号只有 2 秒建议把 seg_len 降到 128否则分段数不足伪峰问题会很严重。3.3 结果可视化与峰值判读双谱是二维复值函数实际使用中看幅值图。画图时用 pcolormesh 直接映射幅值并做归一化避免动态范围过大把弱峰压没。import matplotlib.pyplot as plt Bm np.abs(B) Bm / Bm.max() fig, ax plt.subplots(figsize(6, 5)) mesh ax.pcolormesh(freqs, freqs, Bm, shadingauto, cmapviridis) ax.set_xlabel(f1 (Hz)) ax.set_ylabel(f2 (Hz)) ax.set_title(Bispectrum magnitude (normalized)) fig.colorbar(mesh) plt.tight_layout() plt.show()正常情况下图上会在 (20, 30) 和 (30, 20) 两个位置出现明显峰因为双谱关于对角线对称。50 Hz 本身不会独立出现成峰它只以坐标和的形式参与。76 Hz 处即使功率谱有谱线双谱这里也不该有稳定峰。判读时先找“两个频率之和等于第三个频率”的坐标点再看峰是否落在物理上合理的频率组合上。很多初学者看到图上有亮斑就报耦合忽略了亮斑可能是直流泄漏或分段不足造成的伪峰下面的避坑章节会专门展开。4. 双谱分析四个必调参数分段长度、重叠率、窗函数与归一化4.1 分段长度频率分辨率与统计方差的跷跷板分段长度 seg_len 是双谱分析里最关键的参数。它同时决定频率分辨率和独立分段数seg_len 越大频率分辨率越高Δf fs / seg_len但分段数越少估计方差越大seg_len 越小段数越多方差降下来了但频率分辨率变差两个靠近的耦合峰可能糊成一团。我的经验区间是 128 到 1024 点。采样率 200 Hz 时128 点对应 0.78 Hz 分辨率、约 1.56 Hz 的峰展宽适合看窄带耦合1024 点分辨率 0.2 Hz但 16 秒数据只有约 7 段伪峰风险明显上升。实际调参时先跑 seg_len256看峰是否清晰如果峰宽得厉害再加大如果方差大、背景不平再减小。频率轴间隔与 FFT 点数绑定。直接法里 rfft 的输出点数就是 seg_len//2 1不需要单独设 FFT 长度。有些间接法实现允许对二维累积量补零再做 FFT那会让频率轴变密、峰更平滑但不会增加真实分辨率。我一般把段长直接当 FFT 长度用省事且不容易出错。4.2 重叠率段数越多越好吗重叠率直接控制分段数。0% 重叠时段数最少方差最大50% 到 75% 重叠是工程上最常用的区间。一个常见的直觉误区是“重叠越多段数越多估计越稳”实际并非线性关系——相邻段共享大量数据它们之间高度相关有效独立段数远小于表面段数。重叠率超过 75% 后增加段数带来的方差收益几乎可以忽略计算量却成倍上涨。我在对双谱做阈值判断时会按独立段数来估算置信水平。独立段数近似为 N_total / seg_len / (1 - overlap)段数少于 10 时双谱背景噪声会出现明显的随机起伏这时候任何峰都要多问一句换个随机种子还在不在。重叠率的另一个作用是在分段长度较短时补救段数不足比如只有 4 秒数据、seg_len256 时50% 重叠能凑出约 7 段勉强够看趋势但不适合下结论。4.3 窗函数别用矩形窗也别迷信高阶窗窗函数对双谱的影响比对功率谱更复杂。加窗作用于每个分段的时域数据再经 FFT 进入三元乘积所以旁瓣泄漏会被三次方放大。矩形窗的旁瓣在某些频率组合下会产生虚假的三频乘积这是我实际踩过的坑。汉宁窗是默认选择旁瓣衰减和主瓣宽度平衡较好Hamming 窗类似但旁瓣稍高一点Blackman 窗旁瓣更低但主瓣更宽会抹掉相邻的耦合峰。窗函数还有一个容易被忽略的细节加窗后每段 FFT 的幅值会损失。定量比较不同信号的双谱幅值时必须做能量归一化。第 3 章代码里已经写了/ win.sum()如果你手边的实现没做这一步同一段信号加不同窗会得到量级不一样的幅值这不是算法问题是归一化缺失。4.4 归一化双谱幅值不能直接跨数据比较直接法得到的双谱幅值正比于信号幅度的三次方。同一个耦合机制信号幅度翻倍双谱幅值涨八倍。这带来两个实际问题一是同一组数据里强分量旁边的弱耦合峰会被动态范围压掉二是不同采集条件的数据之间幅值完全不可比。我一般先对原始信号做标准化也就是除以标准差然后再进双谱估计。这样至少能把同一个量级的数据放在一起看。如果要做严格一点的耦合强度量化需要进一步计算双相干系数把双谱幅值除以三个对应频率功率谱乘积的平方根结果落在 0 到 1 之间这个我在最后一章展开。记住一点未经归一化的双谱幅值只能用来找峰的位置不适合作为耦合强度的度量。参数建议范围影响备注分段长度128~1024 点频率分辨率与方差先试 256按峰形调整重叠率50%~75%有效独立段数超 75% 收益很小窗函数汉宁 / Hamming旁瓣泄漏被三次方放大矩形窗慎用幅值归一化信号标准化 / 窗能量补偿幅值可比性定量分析必须做5. 双谱分析避坑指南四个踩坑记录与排查思路5.1 伪峰把独立频率当成相位耦合现象双谱图上出现一个明显的峰坐标恰好满足 f3 f1 f2于是下结论说存在二次相位耦合。但把信号拆开看f3 是独立源产生的相位与 f1、f2 没有任何关系。原因有限样本下X(f1)·X(f2)·X*(f1f2) 的平均值不会严格为零。分段数少于 10 段时随机相位凑出来的伪峰幅值可以达到真实耦合峰的 30% 到 50%。我在验证算法时用 76 Hz 独立分量就复现过这个现象。解决先用第 3 章的测试信号确认算法正常再对真实数据增加分段数到 20 段以上。如果数据不够长把重叠率提到 75%。仍然存疑时做随机化检验把频谱相位打乱后重算双谱重复 50 次取 95% 分位数作为阈值峰值低于阈值的都当伪峰处理。5.2 直流与趋势项对角线方向的亮脊现象双谱图沿主对角线方向出现一条亮带几乎覆盖所有频率组合把真正想看的小峰全压住了。信号时域上看着有缓慢漂移均值不为零。原因三阶累积量对均值非常敏感。信号里的直流分量和线性趋势在频域表现为零频附近的强分量它们与其他频率的三元乘积在低频谱区形成一条脊。很多人只对原始数据做了零均值化忘记去趋势漂移项依然在贡献能量。解决分段之前额外做一次 detrend。NumPy 里 scipy.signal.detrend 就能处理线性趋势做完之后再减均值。如果数据里有明显的工频干扰或基频泄露还可以先做个高通滤波把 5 Hz 以下的内容去掉再进双谱估计。5.3 高重叠率没有换来方差改善现象把重叠率从 50% 调到 90%分段数翻了几倍但双谱背景噪声几乎没有下降计算时间倒是涨了不少。原因相邻段高度相关90% 重叠时每段只新引入 10% 的数据新增段携带的独立信息很少。直接法双谱估计是分段平均方差改善由独立段数决定而不是表面段数决定。有效独立段数大致是 N_total / seg_len / (1 - overlap)overlap 接近 1 时这个数增长极慢。解决把重叠率限制在 75% 以内。想真正降方差优先加长原始数据或者缩短分段长度这两条路都比无限提高重叠率有效。检查方差是否改善时看背景噪声的绝对水平别看峰高——峰高对平均本来就不敏感。5.4 双谱幅值跨数据不可比现象两组数据同样存在 QPC第一组双谱峰高 0.8第二组峰高只有 0.1误判为耦合强度差八倍。实际上只是第二组信号整体幅值小了一半。原因双谱幅值正比于信号幅度的三次方未做任何归一化时幅值差异完全被采集增益、传感器灵敏度等非物理因素主导。直接拿双谱幅值做数据间比较翻车是必然的。解决先用信号标准差做标准化再计算双相干系数。双相干系数的分母包含对应三个频率的功率谱等于把幅度依赖同时消除结果落在 0 到 1 之间才能跨数据集比较。具体公式和代码见下一章。6. 从双谱到诊断结论切片、幅值比与双相干系数三个实战技巧6.1 频率切片把二维平面降成一维曲线双谱是二维的批量分析时逐图看太慢。固定一个频率沿着另一个方向取切片就能把二维图变成一维曲线。比如固定 f1 20 Hz沿 f2 方向看 30 Hz 处有没有峰这就是一次针对性的耦合检查。对大量样本做同一频率对的切片还能画出耦合强度随样本的变化曲线。代码上就是Bm[:, 20]这样一行索引的事但这一步能把双谱从“画图看”变成“算指标看”。对角切片 B(f, f) 对应 f 与自身耦合形成 2f 和频适用于检测二次谐波类非线性。实际诊断里固定频率切片比对角切片更常用因为故障特征频率往往是已知的直接查对应坐标点最直接。6.2 双相干系数带阈值的耦合强度双相干系数把双谱幅值除以对应频率功率谱乘积的平方根是所有实测数据诊断时最值得用的归一化指标b(f1, f2) |B(f1, f2)| / sqrt(P(f1) · P(f2) · P(f1f2))它的取值在 0 到 1 之间越接近 1 说明相位耦合越稳定。工程上我习惯以 0.3 作为“存在明显耦合”的底线低于这个值的峰基本不可信。这个指标的另一个价值是抗幅值干扰不同采集条件下得到的 b 值可以直接放在同一张图里比较不需要再做任何额外归一化。6.3 短时双谱看耦合随时间的变化双谱平均把整段时间的耦合信息揉在一起耦合只出现在某几秒时会被稀释。短时双谱的思路是用滑动窗口沿时间轴一帧一帧算双谱只提取目标频率对的幅值得到一条“耦合强度随时间变化”的曲线。窗口长度可以取前面定的 seg_len步长取窗口的一半。这本质上是在双谱基础上做了个时频分析对非平稳信号特别有用。我早期在一批轴承数据上直接跑全时双谱结果峰值一直若隐若现后来改成短时双谱才发现耦合只在特定转速区间出现。那次之后我的习惯是任何数据都先跑一版短时双谱确认耦合在时间轴上确实稳定再做全时段平均。这个习惯帮我省掉了不少误判。希望帮到你。本文还有配套的精品资源点击获取