简介：一份高阶谱分析HOSA完整工具箱资源（兼容MATLAB R12），面向信号处理、机械故障诊断、生物医学与通信工程领域的研究者。该工具箱通过计算高阶累积量、矩函数及相关函数，支持一至四阶谱（含双谱、三谱）计算、多变量联合分析、数据预处理与图形界面，可揭示非线性非高斯信号中的相位耦合、二次/三次非线性特征。包内共147个文件，以.m源码为主，另有.mat数据、.dat测试数据、.c底层程序及说明文档，总大小仅2.74MB，携带方便。目前已有112人学习下载。借助太阳黑子、猞猁种群等内置实测数据，读者可快速复现双谱估计、谐波恢复、非最小相位系统辨识等典型流程，并灵活迁移至脑电信号非线性分析、机械振动故障检测、非高斯噪声下的通信信号处理等场景。
1. hosa工具箱与3/2谱：一段噪声信号里藏着的耦合关系，为什么双谱能看见
功率谱看上去干干净净的信号，不一定真没问题。之前我接过一段设备振动数据，频谱上只有几个清晰尖峰，怎么看都正常，可设备就是间歇性出毛病；后来有人提醒我试试双谱，结果在频率分量互相“纠缠”的位置冒出一个明显峰。这就是标题里 3/2谱（双谱）的意义：它不是玄学，而是用三阶统计量在二维频域里找频率分量之间的耦合关系。hosa工具箱（hosa toolbox）是干这件事最顺手的老牌工具，网上常见的 hosa.rar 解压出来就能在 Matlab/Octave 里用。这篇笔记讲清楚它是什么、最小命令怎么跑、参数怎么调、坑在哪儿，适合做故障诊断、水声、生物电信号和雷达信号分析的工程师。
2. 看懂3/2谱：高阶谱的数学逻辑与hosa toolbox能干的三件事
2.1 从三阶累积量到双谱：3/2谱这个名字怎么来的
功率谱是二阶统计量的傅里叶变换，它只回答“信号里有哪些频率分量、各自多强”，完全不关心这些分量之间有没有相位关系。可工程信号里大量有用的特征正好藏在相位关系里：两个频率分量的相位被“锁”在第三个分量上，这在功率谱里一点痕迹都没有。要看到这层关系，就得往三阶统计量上走。
定义上，零均值信号的三阶累积量写成：
C3(m, n) = E{ x(t) · x(t+m) · x(t+n) }
对它做二维傅里叶变换，得到的就是双谱：
B(f1, f2) = Σm Σn C3(m, n) · exp(-j2π(f1·m + f2·n))
双谱是复数，自变量从一维频率变成了二维频率平面（f1, f2）。功率谱是一维的，来源于二阶量；双谱是二维的，来源于三阶量，于是有些中文笔记直接叫它“3/2谱”。这个叫法不算严谨，但很好记：三阶统计量对应二维谱，三对二，就成了 3/2。
高斯噪声的三阶累积量理论上是零，所以只要信号里没有真正的三阶耦合成分，双谱上就不该出峰。这一点比功率谱“抗噪”得多——功率谱里高斯噪声还有一条底噪基线，双谱里高斯噪声直接被平均没了。这也是很多故障特征在功率谱上看不到、在双谱上一看一个准的根本原因。
2.2 幅值、相位、双相干：双谱输出里每一样能拿来干什么
hosa工具箱算出来的双谱是个二维复数矩阵。工程上不要直接盯着复数看，先拆成幅值和相位两部分。幅值反映耦合强度：某个频率组合（f1, f2）上幅值越大，说明这三个频率分量的相位耦合越稳定、越强。相位部分则记录耦合的具体形态，比如二次相位耦合的方向是非线性系统产生的还是纯随机碰出来的。
还有一个更稳的指标叫双相干（bicoherence），它把双谱幅值做了归一化，数值落在 0 到 1 之间。接近 1 说明这个频率组合的耦合关系稳定可信，接近 0 说明那个峰可能只是随机相位凑出来的假象。做诊断时我一般先看双相干确认“是不是真的”，再用双谱幅值看“有多强”。
双谱既然是二维的，那它就能沿某个方向“切”出一维线来用。最常见的是固定 f2 看 f1 的响应，或者沿 f1=f2 的对角线切一个对角切片。这样做的好处是把二维图像压成一维曲线，量化阈值、自动寻峰都变得简单。标题里的 3/2谱如果出现在某篇笔记的代码注释里，多半指的就是这条切片。
2.3 什么时候值得掏HOSA：高斯噪声、二次相位耦合与非线性检测
我一般只在三种情况下掏 hosa工具箱。第一，高斯噪声背景下的弱信号检测，双谱能把噪声压掉，信噪比远好于功率谱；第二，怀疑信号里有二次相位耦合，比如机械故障时出现的调制边带、载波与故障特征频率之间的耦合，双谱直接在对应频率坐标上出峰；第三，鉴别系统非线性，线性系统不会产生新的频率耦合，双谱上出现意料之外的峰，往往说明通道里存在非线性环节。
反过来，有些场景不该用。信号长度太短，分段后凑不出足够多的观测样本，双谱的方差会压不下来，峰全是毛刺；信号本身是平稳线性谱，功率谱已经能说明问题，再上HOSA属于白费力气。最典型的误用是把双谱当“更好的功率谱”，一上来就全频段扫图，结果既看不出峰又不知道怎么解释。正确做法是先有具体怀疑的耦合关系，再让双谱去验证它。
3. 装上hosa工具箱并跑出第一张双谱：解压、路径与最小命令
3.1 解压、目录结构与路径配置
你手里的 hosa.rar 解压后，核心是一批 .m 源文件，部分版本还带编译好的 mex 二进制文件，文件名结尾常见 .mexw64、.mexa64 之类。这些 .m 文件就是工具箱的全部逻辑，mex 文件只是一些核心函数的加速版本。理解这一点很重要，后面避坑章节会讲：mex 文件在旧版编译的可不一定能在新环境里跑。
拿到解压目录后，第一步是把目录加进 Matlab 搜索路径。我一般用绝对路径，避免相对路径在切换工作目录后失效：
addpath('D:/tools/hosa'); % 换成你自己的解压绝对路径 savepath; % 保存路径，下次启动仍有效 which bispecd % 检查是否找到主函数addpath只对当前会话生效，savepath把路径写入 Matlab 的路径配置文件，这样重启以后还能用。如果savepath提示没有写入权限，那就不要挣扎，直接在用户启动脚本里加一行addpath('D:/tools/hosa')，每次启动自动加载。which bispecd是验证手段，它能确认 Matlab 实际找到的是不是 hosa 目录里的那个函数——这一步能提前暴露很多路径冲突问题。
3.2 最小可用流程：造一段二次相位耦合信号，画出双谱
为了不拿真实数据试错，我习惯先造一段已知结构的仿真信号，跑通流程再换真数据。下面这段代码造了一个最经典的二次相位耦合信号，第三个频率分量的相位是前两个的和，双谱一定会在对应坐标上出峰，等于给自己留了一个可校验的答案。
fs = 1000; % 采样率 1000 Hz t = (0:4095)' / fs; % 4 秒信号 fa = 50; fb = 90; % 两个基频分量 phi1 = 0.3; phi2 = 0.7; % 二次相位耦合：第三个分量的频率 = fa+fb，相位 = phi1+phi2 x = cos(2*pi*fa*t + phi1) + cos(2*pi*fb*t + phi2) ... + cos(2*pi*(fa+fb)*t + phi1 + phi2) ... + 0.2 * randn(size(t)); % 加一点高斯噪声 % 直接法双谱估计 nfft = 256; % FFT 点数 wind = 'hamming'; % 窗类型 nsamp = 256; % 每段长度 overlap = 0.5; % 50% 重叠 B = bispecd(x, nfft, wind, nsamp, overlap); % 画出双谱幅值 figure; mesh(abs(B)); xlabel('f1 下标'); ylabel('f2 下标'); zlabel('|B(f1,f2)|');这里bispecd走的是直接法路径：分段、加窗、做 FFT、对每段谱做三重乘积后平均。信号特意加了高斯噪声，功率谱上会多一条噪声底，但双谱上噪声会被平均掉，最后的峰只来自那三个相位锁定的余弦分量。nfft=256决定频率轴的分辨率，nsamp=256表示每一段取 256 个点，overlap=0.5让相邻段有一半重叠，相当于在不加长数据的前提下把段数翻倍，压低估计方差。
3.3 函数家族怎么分工：从累积量到双谱再到三谱
hosa工具箱不是只有一个函数，而是一族。我整理了一份常用分工表，照着选函数基本不会走错：
| 函数 | 全称含义 | 典型输出 | 适用场景 |
|---|---|---|---|
| cumest / cum3est | 累积量估计 | 三阶累积量矩阵 | 先算累积量，再做后续处理 |
| bispecd | 直接法双谱估计 | nfft×nfft 复数矩阵 | 数据较长，追求速度，最常用 |
| bispeci | 间接法双谱估计 | nfft×nfft 复数矩阵 | 噪声较重，先平滑累积量再变换 |
| trispecd | 三谱估计 | 四维复数数组 | 需要看四阶统计量时 |
| bicoher | 双相干估计 | 归一化 0~1 矩阵 | 判断双谱峰是否真实可信 |
直接法和间接法的选择，是刚上手最容易纠结的点。bispecd直接对分段数据做 FFT 再三重相乘平均，速度快，适合数据量大、先把结果跑出来看看的场景；bispeci先估计三阶累积量，再做二维傅里叶变换，因为累积量估计过程本身有平滑效果，所以噪声重的情况下峰更干净，代价是慢一点、参数多一个最大延迟长度。
实际项目里我通常先用bispecd出图，确认有峰之后再用bispeci复核一遍。两个结果互相对得上，这个峰才敢写进报告。trispecd用到得少，但碰上信号四阶耦合、双谱上没反应的情况，它就是你最后一张牌。
4. 把双谱调出可用的峰：bispecd参数、频率轴与降维切片
4.1 直接法参数搭配逻辑：分辨率与方差怎么权衡
bispecd的四个参数看起来简单，搭配起来全是经验。nfft决定频率轴的采样点数，也决定画出来的图有多少个像素点；nsamp是每段长度，直接决定频率分辨率，段越长分辨率越高；overlap增加段与段之间的重叠，等效增加统计平均次数，压低方差；wind控制加窗形态，抑制频谱泄漏。写在这里的参数表是我常用的起点：
| 参数 | 作用 | 起点建议 | 调参方向 |
|---|---|---|---|
| nfft | 频率轴长度 | 256 或 512 | 峰变宽时增大 |
| nsamp | 每段长度 | 与 nfft 相等 | 频率分辨率不够时增大 |
| overlap | 段间重叠比例 | 0.5 | 图太毛糙时提高到 0.75 |
| wind | 窗类型 | 'hamming' | 泄漏严重时换窗 |
最核心的权衡发生在nsamp和方差之间。固定总长度时，段切得越短，可平均的段数越多，方差越小，但每段频率分辨率越差，双谱峰会被抹宽；段切得越长，峰越尖锐，但段数少，方差压不住，图上一片刺。我的经验是先用nsamp=nfft=256、overlap=0.5跑一版，如果峰糊成一团就增大nfft，如果图变成刺猬就增大overlap。
一个很容易忽略的细节是：信号要先做零均值化再做双谱估计。hosa工具箱内部一般会处理，但如果你先对数据做了带通滤波，滤完一定要重新去均值，否则直流泄漏会在双谱低频区域制造一片假峰。去均值就一行x = x - mean(x)，不亏。
4.2 频率轴约定与峰位置换算：别把坐标看反
bispecd返回的 B 是一个 nfft×nfft 的复数矩阵，行和列分别对应两个频率轴，但不同版本的 hosa工具箱对频率归一化的约定不完全一样。常见的有两种：一种把采样频率当成 1，矩阵下标直接乘采样率就是物理频率；另一种把奈奎斯特频率当成 1，下标要乘采样率的一半。这个差异足以让刚上手的人对着一个“位置不对”的峰怀疑人生。
我处理这个问题的方法很简单——用 3.2 节里那段已知频率的仿真信号去校准。fa=50Hz、fb=90Hz、采样率=1000Hz，跑完双谱后用find定位最大峰的行列下标，把下标与理论频率做对比，一次就能确定当前版本的坐标约定：
B_abs = abs(B); [kr, kc] = find(B_abs == max(B_abs(:)), 1); % 第一种约定：下标/nfft*fs f1_est = (kr - 1) / nfft * fs; f2_est = (kc - 1) / nfft * fs; % 如果算出来是 25/45 这种半数值，说明约定是把奈奎斯特当 1，各乘 2 即可这里find定位的是双谱幅值最大的位置，理论上应该落在 (50, 90) 或 (90, 50) 附近，因为双谱对 f1 和 f2 是对称的，两处都会有峰。校准完之后，把换算公式写死在处理脚本开头，后面换真实数据时就不会再迷茫。坐标如果差了正好一半，多半就是归一化约定问题，不是信号问题。
4.3 从二维谱里提一维特征：固定频率切片与对角切片
二维双谱图画出来好看，但做成自动诊断特征就麻烦了——阈值怎么定、峰值怎么全自动提取，都不好办。工程落地上我几乎总是把双谱降维成一维曲线再用。两种切片最常用：固定 f2 在某一个特征频率，看 f1 方向的响应；或者沿 f1=f2 对角线切一刀，看倍频耦合的整体分布。
下面这段代码展示了固定频率切片怎么用，它跟在 3.2 节代码后面直接运行：
% 固定 f2 在 90 Hz 对应的下标：由校准阶段确定 kc_ref kc_ref = round(90 / fs * nfft) + 1; % 按采样率为 1 的约定估算 slice = abs(B(:, kc_ref)); % 取整列，得到一维谱线 % 换算 f1 轴 f1_axis = (0:nfft-1) / nfft * fs; plot(f1_axis, slice); xlabel('f1 (Hz)'); ylabel('切片幅值'); grid on;固定频率切片拿到的这组数，天然就是“某个特征频率和其他所有频率的耦合强度曲线”。做齿轮箱或轴承诊断时，我常把故障特征频率定在 f2 上，然后看 f1 轴出现在哪些频率上有峰，峰出现的位置就能反推是哪两个分量在耦合。这样做的好处是：二维图里的主观判断变成了一维曲线上的固定阈值比较，自动化程度高很多。
对角切片则适合做“整体耦合程度”的粗筛查。把对角线上的幅值累加成一个标量，或者画成一条随频率变化的曲线，哪段频率区间耦合多，一眼就能看出来。它牺牲了频率组合的细节，换来的是稳定性和可解释性。这两个切片方向配合使用，差不多能覆盖 80% 的故障筛查需求。
5. hosa工具箱避坑指南：5个让新手翻车的细节
5.1 坑：addpath 之后重启 Matlab 又找不到函数
现象：第一次addpath能用，重启后which bispecd报“未找到”，脚本直接中断。
原因：savepath没有真正生效。常见情况是 Matlab 安装在系统受保护目录里，普通用户没有写路径配置文件的权限；另一种情况是你在脚本里加了addpath，但脚本退出后路径自然消失。
解决：先执行savepath看有没有报错。报权限错就用用户启动脚本方案，在 startup.m 里写addpath('D:/tools/hosa')，每次启动自动加。加完用which bispecd -all确认找到的是 hosa 目录下那个函数，而不是别的目录里的同名文件。
5.2 坑：新版 Matlab/Octave 加载不了旧 mex 文件
现象：跑某个函数时报“无法加载 mex 文件”，或者提示“不是有效的 Win64 应用程序”。
原因：hosa.rar 里自带的 mex 二进制是在旧版环境或 32 位系统下编译的，新版 Matlab 需要的 .mexw64 格式对不上；在 Octave 里还需要对应平台的 .mex 文件，混用基本必炸。
解决：先找同名 .m 文件替代。hosa工具箱的大部分核心函数都有纯 .m 版本，去掉对 mex 的依赖也能跑，只是慢一点。如果某个函数只有 mex 没有源码，就得找原始源码重新编译；在 Octave 下可以用mkoctfile重新编译。实在不行就换一个以 .m 为主的版本来用，别和二进制文件较劲。
5.3 坑：双谱图密密麻麻全是刺，找不到峰
现象：mesh(abs(B))出来一片毛刺，看起来像随机噪声，和功率谱比没有任何优势。
原因：双谱是三阶统计量的平均结果，它比功率谱需要多得多的样本去压低方差。总数据长度固定时，段切得太长导致段数太少，或者overlap太低，统计平均次数不够，方差压不下来。
解决：先把nsamp降到 128 或 256，增加分段数；再把overlap提高到 0.75，让同一段数据被多次利用；如果整段数据本身太短，前面所有参数都没用，必须先采集更长数据或者降采样把点数压缩到目标频带里。我调参数有个习惯动作：先用仿真信号把参数调到能干净出峰，再保持同样的参数组合上真实信号，这样至少排除“参数不合适”这个变量。
5.4 坑：峰的坐标和理论频率对不上，还差一半
现象：双谱峰明显存在，但坐标换算出来是 25Hz 和 45Hz，而理论值是 50Hz 和 90Hz，刚好差一半。
原因：频率归一化约定不同。有的 hosa版本把采样频率当 1，有的把奈奎斯特频率当 1，后者算出来的物理频率就少一半。坐标方向也可能反，两个频率轴谁对应行、谁对应列，不同版本里不一定一致。
解决：不要对着文档猜，直接用已知频率的二次相位耦合仿真信号校准。跑完用find定位最大峰下标，分别按两种约定换算，看哪个和理论值吻合。校准过后把这个换算关系固定到一个函数里，以后所有数据都走同一个转换。
5.5 坑：和别的工具箱或同名函数冲突
现象：调用cumest或bispecd时报错，错误信息指向的路径不是 hosa 目录，而是另一个工具箱或你自己脚本里定义的同名 m 文件。
原因：Matlab 按搜索路径顺序找函数，同名函数先被找到的先执行。老工具箱函数名很短很短，像cumest这种名字很容易和别人脚本里的函数撞车。
解决：用which bispecd -all列出所有同名版本，确认优先级。把 hosa 目录用addpath(..., '-begin')加到搜索路径最前面，或者干脆用完整路径调用run('D:/tools/hosa/bispecd.m')这种方式绕开路径搜索。项目目录里如果有自己写的同名脚本，尽早改名，别和高阶谱函数较劲。
6. 把hosa双谱用成诊断工具：先校准、再切片、后下结论
6.1 校准先行：每个新环境都要跑一遍自检信号
把 hosa 当工具箱用和当诊断工具用，是两回事。前者只需要函数能跑，后者需要一整套验证闭环。我的固定动作是：每换一台机器、每换一个版本、每新装一次 Matlab，都先跑一遍 3.2 节那段二次相位耦合仿真信号，确认三个东西——函数能正常返回、频率轴换算正确、给定参数组合能干净出峰。这三个确认完了，真实数据的结论才站得住。
6.2 工程验证套路：功率谱打底、双谱定位、切片量化
现在我处理一段新信号的标准套路分四步。第一步跑功率谱，把信号里占主导的频率分量搞清楚，确定可疑的频率组合；第二步用bispecd出双谱图，看可疑组合处有没有峰；第三步用bicoher确认峰是否可信，避免随机相位碰出来的假峰混进来；第四步把固定频率切片做成曲线，设置阈值、提取峰位置，输出可量化的特征值。这套流程下来，双谱不再是“看起来很厉害”的图，而是能写进报告、能对比不同工况的定量指标。
6.3 一个诚实的前提与一个偷懒外的坚持
说句实在话，双谱不是万能的。我见过不少信号，功率谱上什么问题都能看到，双谱反而平平无奇——这种情况就别硬上高阶谱了。真正值得用双谱的，是那些“频率成分都清楚、但就是解释不了故障为什么发生”的信号。耦合关系是 HOSA 独有的信息维度，其他工具给不了。我现在拿到一段陌生信号，第一步永远先跑功率谱，怀疑有耦合再掏 hosa，这个顺序省了我大量无效计算。也希望这套做法帮你在新项目里少走几步弯路，直接落在能用的双谱结果上。
本文还有配套的精品资源，点击获取
