简介本资源为基于Matlab实现的小波相干分析代码包面向本科、硕士及科研人员进行信号处理与教研学习可解决多变量非平稳信号时频相关性分析、相位关系提取等典型问题适用于智能优化、神经网络预测、图像处理等仿真场景。压缩包共109个文件约3.08MB以35个m脚本文件为核心配合15个md说明、12个txt文档、11个png结果图及11个sample样例数据另含html、pdf等辅助资料覆盖代码、数据与运行结果便于对照理解算法流程。资源已在Matlab2014与2019a环境下验证内含运行结果可帮助读者快速复现小波相干谱、交叉小波变换及显著性检验等关键步骤并借助样例数据与说明文档排查常见报错。目前已有142人学习下载适合需要将小波相干方法落地到科研仿真与论文实验中的读者参考使用。1. 从一段非平稳信号说起wavelet-coherence 到底解决什么问题手里有一段两个通道同时采的振动信号或者两路脑电、两组气象序列你想知道它们是不是一起动的。最直接的做法是算相关系数但相关系数只给一个数把时间信息全抹平了。更要命的是很多真实信号是非平稳的——两台设备在某个时间段突然同步共振过了这段又各走各的这种什么时候同步、在哪个频段同步的问题普通相关和 FFT 相干都答不上来。wavelet-coherence小波相干就是冲这个来的它把信号先做连续小波变换再在时频平面上算两路信号的相干谱最后给你一张横轴时间、纵轴频率、颜色表示相干强度的图。matlab 里这套东西有现成函数但真正落地时参数怎么设、边界怎么处理、显著性怎么判才是决定结果能不能信的关键。这篇就按能跑通、能复现、能避坑的路子把 wavelet-coherence 的 matlab 实现讲透适合做信号分析、故障诊断、神经科学、地球物理这类需要看两路信号时频耦合关系的人。2. wavelet-coherence 的数学骨架与 matlab 落地选型2.1 从连续小波变换到交叉小波谱小波相干不是凭空冒出来的它建立在连续小波变换CWT之上。对信号 x(t)其 CWT 定义为W_x(a,b) ∫ x(t) · (1/√a) · ψ*((t-b)/a) dt其中 a 是尺度对应频率b 是平移对应时间ψ 是小波母函数* 表示复共轭。matlab 里cwt函数直接给你 W_x但要注意它默认返回的是尺度域转频率得自己算。有了两路信号的 CWT交叉小波谱就是 W_xy W_x · W_y*。小波相干定义为R²(a,b) |S(a⁻¹ · W_xy)|² / ( S(a⁻¹ · |W_x|²) · S(a⁻¹ · |W_y|²) )这里的 S 是平滑算子作用在时间-尺度两个方向上。这个平滑是整套方法里最容易被忽略、也最容易翻车的地方——不平滑的话相干值恒等于 1因为单点复数的模平方除以自身模平方就是 1。所以平滑窗的选择直接决定相干图有没有意义。我一般用 Torrence Compo 那套经典实现平滑窗在时间方向用 Gaussian尺度方向用滑动平均两个方向的窗长要匹配小波的时频分辨率。matlab 官方没有内置 wavelet coherence 函数网上流传的wtc.m是 Grinsted 工具箱里的这是最常见的做法。2.2 小波母函数怎么选Morlet 不是唯一答案选母函数是第一个决策点。常见选项母函数适用场景时频分辨率特点Morlet振荡信号、耦合分析频率分辨率好时间分辨率中等Paul瞬态、突变检测时间分辨率好DOG边缘、奇异点时间分辨率极好频率差Morse可调参数灵活通过 β、γ 调时频权衡做相干分析绝大多数情况选 Morlet因为你要看的是两路信号在某个频段的相位关系频率分辨率比时间分辨率更重要。Morlet 的关键参数是中心频率 ω0一般取 6这是 Torrence Compo 推荐的平衡点。ω0 越大频率分辨率越高但时间定位越糊ω0 越小反过来。我试过 ω06 和 ω012 对比后者在低频段相干带更清晰但高频瞬态耦合会被抹掉。matlab 里如果用cwt可以指定amorMorlet 的近似或者自己构造。但要做相干建议直接用 Grinsted 工具箱的wt.m因为它和后面的wtc.m平滑逻辑是配套的混用不同来源的 CWT 和相干函数平滑窗对不上结果会偏。2.3 尺度序列与频率映射别让频率轴骗了你CWT 的尺度 a 和实际频率 f 的关系是 f fc / (a · Δt)其中 fc 是母函数的中心频率Δt 是采样间隔。Morlet 的 fc ≈ ω0 / (2π)当 ω0 较大时近似成立。很多人直接拿尺度当频率画图结果频率轴是错的低频段被压缩、高频段被拉伸看着相干带位置全偏。正确做法是先定频率向量再反算尺度。比如采样率 100 Hz想看 0.5–20 Hz就先生成对数分布的频率向量再用 a fc / (f · Δt) 得到尺度。matlab 代码fs 100; % 采样率 Hz dt 1/fs; freq logspace(log10(0.5), log10(20), 60); % 关心的频率 omega0 6; % Morlet 中心频率 fc omega0 / (2*pi); % 母函数中心频率近似 scale fc ./ (freq * dt); % 尺度序列这段的关键是freq用对数分布因为小波在频率域本身就是对数分辨的线性分布会导致高频段尺度点过密、低频段过疏。scale算出来后传给wt.m返回的相干矩阵行对应尺度、列对应时间画图时纵轴用freq而不是scale。提示如果 ω0 取得小比如小于 5fc ≈ ω0/(2π) 这个近似误差会变大最好用数值方法求母函数的实际峰值频率。3. 用 matlab 跑通 wavelet-coherence 的最小可复现流程3.1 造一段已知耦合的合成信号做验证上手第一件事不是拿真实数据跑而是造一段你已知答案的信号验证代码对不对。经典做法两个信号共享一个 10 Hz 分量但只在中间一段时间耦合其余时间各自独立。fs 100; T 20; t 0:1/fs:T-1/fs; N length(t); % 公共分量只在 8-14 秒出现 common sin(2*pi*10*t) .* (t 8 t 14); % 信号1公共分量 5Hz 独立分量 噪声 x common 0.8*sin(2*pi*5*t) 0.3*randn(1,N); % 信号2公共分量 3Hz 独立分量 噪声 y common 0.8*sin(2*pi*3*t) 0.3*randn(1,N);这段代码的意图很明确common是两路共享的 10 Hz 振荡只在 8–14 秒存在其余时间 x 和 y 各有自己的频率成分。跑完相干分析后你应该在 10 Hz、8–14 秒区域看到高相干其他区域相干低。如果结果不是这样说明平滑、尺度或频率映射有问题。randn加的噪声幅度 0.3 是经验值太小看不出抗噪性太大把耦合淹没0.3 左右能兼顾。3.2 调用 wtc 计算相干并画图假设你已经把 Grinsted 工具箱的wt.m、wtc.m、smooth.m等放在路径里主流程如下% 参数设置 omega0 6; freq logspace(log10(0.5), log10(20), 60); scale (omega0/(2*pi)) ./ (freq * (1/fs)); % 计算小波变换 [Wx, period] wt(x, scale, morlet, omega0); [Wy, ~] wt(y, scale, morlet, omega0); % 计算小波相干 [WCT, coi, sig95] wtc(Wx, Wy, mcc, 0); % 画图 figure; imagesc(t, freq, abs(WCT)); set(gca, YDir, normal, YScale, log); hold on; plot(t, 1./coi, w--, LineWidth, 1.5); % 影响锥 contour(t, freq, sig95, [1 1], k, LineWidth, 1); % 95% 显著线 colormap(jet); colorbar; xlabel(时间 (s)); ylabel(频率 (Hz)); title(小波相干谱);wtc的第三个参数mcc表示用 Monte Carlo 方法估计显著性最后一个参数 0 是不做相位显示要相位就设 1。coi是影响锥cone of influence在边界处小波变换受边缘效应影响coi 以外的区域不可信必须用白虚线标出来。sig95是 95% 显著性等值线只有落在里面的高相干区才值得讨论。注意wtc返回的WCT是复数取abs得到相干幅值 R²取angle得到相位差。画相干图用abs分析领先滞后关系用angle。3.3 影响锥与显著性检验两个不能省的步骤影响锥COI是很多人图省事直接忽略的东西。小波在时间边界处会看到零填充导致边界附近的能量虚高相干值也跟着虚高。COI 的边界大致是 e-folding 时间对 Morlet 约等于 √2 · scale。matlab 里wtc直接返回 coi你只要画出来并声明coi 以外不解释就行。显著性检验更关键。小波相干的背景谱不是均匀的红色噪声AR1下相干值的分布和白色噪声完全不同。wtc用 Monte Carlo 模拟生成显著性阈值你需要传入两路信号的 AR1 系数。如果懒得估可以用mcc让它自己估但样本短的时候估不准。我的习惯是先用arfit或简单的 Yule-Walker 估出 AR1 系数再传给wtc% 估计 AR1 系数 ar1_x aryule(x - mean(x), 1); ar1_y aryule(y - mean(y), 1); alpha_x -ar1_x(2); % aryule 返回 [1, -a1] alpha_y -ar1_y(2); [WCT, coi, sig95] wtc(Wx, Wy, mcc, 0, alpha_x, alpha_y);aryule返回的第二个系数是 -a1取负号才是 AR1 系数。这一步不做显著性线会偏松或偏紧导致你误判耦合区域。合成信号验证时如果 10 Hz 耦合区没被 sig95 圈住八成是 AR1 系数没传对。4. 真实数据落地时的参数调优与结果解读4.1 平滑窗长度相干图的分辨率旋钮平滑窗是 wavelet-coherence 里最玄学的参数。窗太短相干值噪声大满图都是假的高相干斑点窗太长真实的短时耦合被抹平相干带糊成一片。Grinsted 工具箱默认在时间方向用 Gaussian 窗尺度方向用滑动平均窗长和小波尺度挂钩。我的经验法则时间方向平滑窗长取 1–2 倍小波尺度对应的时间跨度尺度方向取 0.5–1 倍尺度。具体调的时候先跑默认值看合成信号验证图如果 10 Hz 耦合区的相干带比预期宽就缩短时间窗如果噪声区出现大量孤立高相干点就加长尺度窗。matlab 里wtc的平滑是内置的改窗长要动smooth.m里的参数或者自己重写平滑。更灵活的做法是自己实现平滑% 自定义平滑时间方向 Gaussian尺度方向滑动平均 function S my_smooth(W, scale, dt, t_win, s_win) % t_win: 时间窗长秒s_win: 尺度窗长尺度单位 [ns, nt] size(W); S zeros(ns, nt); for i 1:ns % 尺度方向滑动平均 lo max(1, i - round(s_win/2)); hi min(ns, i round(s_win/2)); Ws mean(W(lo:hi, :), 1); % 时间方向 Gaussian tw round(t_win / dt); g exp(-(-tw:tw).^2 / (2*(tw/2)^2)); g g / sum(g); S(i, :) conv(Ws, g, same); end end这段代码把两个方向的平滑拆开t_win和s_win可以独立调。conv用same保持长度不变边界处会有衰减配合 COI 一起看就行。调参时建议固定一个、动另一个别同时改否则你不知道是哪个起了作用。4.2 相位分析谁领先谁滞后相干幅值只告诉你耦合强不强相位告诉你怎么耦合。angle(WCT)给出两路信号在时频点的相位差正值表示 x 领先 y负值反之。但相位在相干值低的地方没有意义所以解读相位前先用相干幅值做掩膜只保留 R² 0.5 的区域。画相位箭头图是常见做法在相干高的网格点上画箭头箭头方向表示相位差。matlab 里用quiver% 相位箭头只在相干 0.5 处画 [R2, ~] deal(abs(WCT), angle(WCT)); mask R2 0.5; [ti, fi] meshgrid(1:length(t), 1:length(freq)); ti ti(mask); fi fi(mask); u cos(phase(mask)); v sin(phase(mask)); quiver(t(ti), freq(fi), u, v, 0.5, k);箭头指向右表示同相相位差 0指向左表示反相相位差 π指向上或下表示 90° 相位差。这套解读在分析两个振荡器耦合、脑区同步、气候指数遥相关时都用得上。注意箭头密度别太高否则糊成一团一般每隔 3–5 个网格点画一个。4.3 多尺度耦合的频段划分真实信号的耦合往往不是单一频段而是多个频段同时存在。比如机械振动里低频的转频和高频的啮合频率可能同时耦合。这时候别指望一张相干图看全要分频段统计。做法是把频率轴分成几个 band对每个 band 内的相干值做时间平均得到频段-时间的耦合强度曲线bands [0.5 2; 2 8; 8 15; 15 20]; % 频段边界 for k 1:size(bands,1) idx freq bands(k,1) freq bands(k,2); band_coh(k, :) mean(abs(WCT(idx, :)), 1); end plot(t, band_coh); legend(0.5-2Hz,2-8Hz,8-15Hz,15-20Hz);这样能看出哪个频段的耦合在哪个时间段最强比单张相干图信息量大。mean用算术平均就行别用几何平均相干值接近 0 时几何平均会数值不稳。5. wavelet-coherence 实战避坑5 个血泪教训5.1 相干值恒为 1忘了平滑现象跑完相干图整张图全是红色R² 处处接近 1看不出任何结构。原因小波相干定义里如果不做平滑单点的 |W_xy|² / (|W_x|²·|W_y|²) 在数学上恒等于 1Cauchy-Schwarz 等号条件。平滑是让相干值有意义的唯一途径。解决确认wtc内部平滑被正确调用或者自己实现平滑时检查窗长是否非零。如果用的是自己拼的代码检查平滑是不是只作用在幅值上而漏了交叉项。5.2 边界高相干假象没画影响锥现象相干图左右两端出现大片高相干但物理上那段数据不该有耦合。原因小波在时间边界处卷积到零填充区能量虚高相干值跟着虚高。这是边缘效应不是真实耦合。解决始终画 COI 白虚线解读时声明 COI 以外区域不讨论。如果边界区域恰好是你关心的用镜像延拓或反射填充预处理信号但要在文中说明。5.3 频率轴错位尺度当频率用现象相干带的位置和预期频率对不上低频段被拉长、高频段被压缩。原因直接把 CWT 的尺度序列当频率轴画图没做 f fc/(a·Δt) 的映射。解决先生成频率向量反算尺度画图时纵轴用频率。检查方法拿一个已知频率的正弦信号做 CWT看峰值是否落在正确频率上。5.4 显著性线偏松AR1 系数没传现象sig95 圈住一大片区域很多明显是噪声的地方也被判为显著。原因wtc默认假设白噪声背景真实信号多是红色噪声白噪声假设下显著性阈值偏低。解决用aryule或arfit估 AR1 系数传给wtc。样本短于 200 点时 AR1 估计不稳可以用多段平均或直接设一个保守的经验值如 0.5。5.5 相位解读翻车低相干区看相位现象相位箭头图看起来有规律但物理上说不通。原因相干值低的地方相位差是随机噪声没有物理意义。在 R² 0.5 的区域解读相位等于看噪声。解决相位分析前先用相干幅值做掩膜只保留 R² 0.5或 0.6的区域。掩膜阈值没有绝对标准但低于 0.5 的相位基本不可信。6. 把 wavelet-coherence 用稳的几个进阶习惯合成信号验证要成为肌肉记忆。每次换数据、换参数、换母函数先拿一段已知耦合的合成信号跑一遍确认代码没退化。我习惯把合成信号验证写成一个独立脚本改参数后先跑它通过了再上真实数据。这一步花两分钟能省掉后面几小时的这结果到底对不对的纠结。相干图的配色和显示范围也有讲究。imagesc默认把数据范围拉满噪声区的小波动会被放大成显眼颜色。我一般把caxis固定在 [0 1]因为相干值理论上就在 0–1 之间固定色标后不同图之间可以直接对比。如果关心的是弱耦合可以把上限压到 0.8但要在图注里写明。多组数据对比时别只放相干图把相干值的分布也画出来。用histogram(abs(WCT(:)), 50)看整体分布如果大部分值集中在 0.3 以下说明耦合很弱图上的高相干斑点可能只是噪声涨落。这个习惯帮我避过好几次看着有耦合、实际是噪声的误判。最后说个我自己的教训早期做相干分析总想把参数调到最好看结果调出一张漂亮的图但换个数据集就复现不了。后来改成固定一套参数Morlet ω06、对数频率 60 点、时间窗 1.5 倍尺度、尺度窗 0.5 倍尺度、AR1 系数实测只在有明确理由时才动结果的可比性和可复现性一下子上来了。参数不是越调越准是越固定越可信。希望帮到你。本文还有配套的精品资源点击获取