简介一个基于C语言实现的小波变换程序包面向信号处理、图像分析领域的工程师与科研人员以多分辨率分析方式同时提取时间与频率特征弥补傅里叶变换丢失时域信息的不足。压缩包仅1个c文件大小约2KB核心代码cwt.c中提供输入信号sig、尺度向量scales、小波母函数wname等关键参数入口支持自定义尺度范围与小波类型便于移植到实时或离线任务中。已有464人学习下载。通过学习这份代码用户可以掌握连续小波变换的C语言实现思路理解尺度变换与卷积计算流程并以此为模板扩展更多小波基函数或集成至自己的信号处理项目。整体精简、无冗余文件适合具备一定C语言基础和小波理论知识的读者快速阅读与二次开发。1. 拿到 cwt.zip 却要在 C 语言里做小波变换先把 cwt.m 的思路放下手边有 cwt.zip、cwt.m甚至已经跑出来几张时频图却被要求“用 C 语言实现一遍 cwt”——这个需求在嵌入式信号处理里很常见。原因不外乎两种目标平台没有 MATLAB 运行时或者算法要进实时采集链路不能再靠脚本算完再画图。这里有个反直觉的结论cwt 的 C 实现本身不复杂复杂的是把尺度、采样率、小波中心频率这三个参数在代码里对齐否则同样的信号在 cwt.m 和 C 程序里出来的峰值位置、幅值比例会差很多。MATLAB 里 cwt 函数返回的是一个现成的复数矩阵换到 C 语言你得自己决定缓冲区怎么开、核函数怎么截断、边界怎么处理。这篇文章按“离散公式 - 最小可运行代码 - 参数调优 - 和 cwt.m 对齐”的顺序把 C 语言场景下实现小波变换的路线讲完整。2. cwt 的离散化把尺度 a、平移 b、采样率 fs 凑到同一个数组里2.1 从连续积分到数组求和三个替换必须同时做连续小波变换的定义式是W(a,b) 1 / sqrt(a) * ∫ x(t) * ψ*((t - b) / a) dt其中 ψ* 是小波母函数的共轭a 是尺度b 是平移。这个式子在 C 语言里不能直接求因为 t 是连续时间而信号 x 是一串按采样周期排列的离散点。把它变成程序能跑的求和要同时做三个替换第一积分换成求和dt 变成差分项第二小波函数的自变量 (t - b) / a 变成 (k - b) * T_s / a其中 T_s 1 / fs第三对每个平移 b只对落在小波核有效范围内的数据点累加超过核截断范围的部分不做贡献。三者缺一个结果就不是 cwt而是某种带通滤波或者卷积变体。C 语言里实现时我一般把尺度 a 的单位直接定成“秒”核函数的横轴就用采样点序号换算成秒。这样最小可运行版本里不需要引入复杂的单位换算逻辑。double dt 1.0 / fs; // fs 是采样率单位 Hz // 尺度 a 的单位是秒核自变量 t 的单位是秒 double t (double)(i - half) * dt / a;这里 scale 用秒而不是用“采样点数”原因是后续和频率换算时更直观。对于一个中心角频率为 w0 的 Morlet 小波它对应的中心频率是 w0 / (2 * PI * a)这里 a 用了秒频率直接就是 Hz。如果 a 换成采样点数公式里还得再乘一个 fs容易在移植 cwt.m 时对不上数。2.2 尺度序列怎么定对数间隔加最小尺度约束cwt 的尺度不是随便给几个数一般要按对数均匀分布铺开。原因是小波的尺度在低频端变化慢在高频端变化快线性间隔会让低频端分辨率不够高频端又冗余。常用做法是void log_spaced_scales(double *scales, int nscales, double s_min, double s_max) { double ratio pow(s_max / s_min, 1.0 / (nscales - 1)); scales[0] s_min; for (int i 1; i nscales; i) { scales[i] scales[i - 1] * ratio; } }这个函数做的事情很直接从最小尺度 s_min 开始每次乘一个固定比率直到 s_max。比率按 nscales 的个数均匀拆分。这样在双对数坐标下尺度点是一串等间隔的点跟 cwt.m 里大部分实现默认的排列方式一致。s_min 不能小于信号最高频率对应的半个周期。假设信号最高有效频率是 150 HzMorlet 用 w06那么最小尺度大约是 6 / (2 * PI * 150)约 0.0064 秒。如果再小小波核的振荡频率超过奈奎斯特或者超过信号自身带宽算出来的只有核的高频尾巴没有实际物理意义。2.3 为什么 cwt 的核用复数而不是实数用 Morlet 小波做 cwt核函数天然是复数的因为 Morlet 本质上是一个复指数乘高斯窗。实数核也可以做小波变换但它们通常是差分型小波比如 Mexican Hat得到的系数是实数频率选择性比复数核差。Morlet 核的表达式是ψ(τ) π^(-1/4) * exp(i * w0 * τ) * exp(-τ^2 / 2)其中 τ 是核的时间轴。在 C 语言里这个表达式对应两个组成部分复指数 e^(iw0τ) 决定振荡频率高斯窗 exp(-τ²/2) 决定时间局域性。w0 越大振荡越密集频率分辨率越好但时间分辨率越差。cwt 里对这个参数的处理直接决定了变换结果的质量这也是后文参数调优里要重点展开的内容。复数核带来的一个直接后果是C 语言实现时必须用complex double或手写复数结构体。如果只把核存成实部就丢了相位信息幅值计算也不准确。2.4 核截断高斯窗决定了核的有限长度理论上 Morlet 核在时间轴上无限延伸但高斯窗 exp(-τ²/2) 衰减很快。cwt 在实现时一般截断到 ±4 个标准差范围。τ 的标准差是 1所以截断到 τ ∈ [-4, 4] 时窗函数尾部只有 exp(-8) ≈ 0.000335完全可以忽略。每个尺度的核长度不是固定的它跟尺度成正比。尺度越大核在时间上拉得越宽覆盖的采样点数越多。这个特性让某些优化变得不直观固定 N 的情况下小尺度卷积很快大尺度卷积可能要遍历数倍于 N 的核长度。边界问题也因此在大尺度下更严重。下一节实现代码时核长度按4 * scale / dt来取再用奇数保证关于中心对称。3. 用 C 语言实现 cwt 的最小可用版本核缓存加卷积循环3.1 数据存储结果矩阵用一维数组核用 complex doubleC 语言没有 MATLAB 那种二维复数矩阵所以要自己定义内存布局。结果矩阵我一般开成一个一维数组长度为nscales * n第 s 个尺度、第 b 个时间点存放在out[s * n b]。这样不管之后是写文件、传参还是在 GUI 里显示都是一段连续内存方便用memcpy或直接fwrite导出也方便在 OpenMP 里按尺度拆分并行。double complex *out calloc((size_t)nscales * n, sizeof(double complex)); if (!out) { /* 处理分配失败 */ }分配时用calloc而不是malloc的好处是如果第一个尺度的卷积因为边界问题少算了几项未初始化区域直接是零不会出现随机复数噪声把调试过程带偏。3.2 生成 Morlet 共轭核一个尺度的完整核模板cwt 公式里用的是 ψ*即共轭。所以在生成核时直接生成共轭形式能省掉卷积循环里每一点都取共轭的操作。下面这个函数把某一尺度下的 Morlet 共轭核写入预分配缓冲区。static void morlet_conj_kernel(double complex *psi, int npoints, double scale, double w0, double dt) { int mid npoints / 2; double norm 1.0 / sqrt(scale) * pow(M_PI, -0.25); for (int i 0; i npoints; i) { double tau (double)(i - mid) * dt / scale; double env exp(-0.5 * tau * tau); double re env * cos(w0 * tau) * norm; double im -env * sin(w0 * tau) * norm; // 取共轭后的虚部 psi[i] re I * im; } }逻辑说明i - mid把缓冲区下标换算成相对中心的偏移量tau是偏移量对应的连续时间除以scale后代入母函数。env是高斯窗决定核的衰减cos是实部-sin是共轭后的虚部。norm是尺度归一化系数保证不同尺度下的能量一致。参数说明npoints必须是奇数否则核不对称卷积结果会产生一个采样点的相位偏移w0取 6 时Morlet 小波近似满足允许条件直接用于滤波不会引入明显的均值漂移scale是秒不是样本数这样tau永远是无量纲时间比较稳妥。3.3 直接卷积核心遍历尺度和采样点内层扫核void cwt_direct(const double *x, int n, double fs, const double *scales, int nscales, double w0, double complex *out) { double dt 1.0 / fs; size_t stride (size_t)n; for (int s 0; s nscales; s) { int len (int)(4.0 * scales[s] / dt) 1; if (len 5) len 5; if (len % 2 0) len; if (len 4 * n) len 4 * n; double complex *psi malloc((size_t)len * sizeof(double complex)); if (!psi) continue; morlet_conj_kernel(psi, len, scales[s], w0, dt); int half len / 2; double complex *row out (size_t)s * stride; for (int b 0; b n; b) { double complex acc 0.0 0.0 * I; for (int k 0; k len; k) { int idx b k - half; if (idx 0 || idx n) continue; acc x[idx] * psi[k]; } row[b] acc * dt; } free(psi); } }这段代码是最直接的 cwt C 语言实现没有任何优化但可读性最高适合作为移植 cwt.m 时核对结果的基准。外层循环按尺度展开中层循环按时间平移点展开内层循环把核和信号逐点乘加。边界处直接跳过越界索引等价于零填充延拓这是最简单但误差最大的方式。参数说明len取4 * scale / dt对应 τ ∈ [-4, 4]尾部误差在千分之一左右stride用n保证每个尺度在输出里是连续一块最后乘dt是因为离散积分必须带上采样间隔否则结果量纲不对核长度超过4*n时强制截断避免尺度太大导致核比信号还长内存和计算量失控。3.4 复杂度评估什么时候这个实现能直接上线对比项直接卷积FFT 卷积每个尺度复杂度O(N * L_k)O(N log N L_k log L_k)核远短于信号时最快无额外内存固定开销高不划算核与信号同量级时明显变慢优势明显额外内存只需一组核需要 FFT 缓冲区和频域掩模嵌入式定点平台舍入误差容易控制动态范围难控制当 N1024、尺度 30 个、平均核长 200 点时直接卷积大约做 614 万次复数乘加普通桌面 CPU 在毫秒级完成完全可以直接用。只有当 N 到几万、或嵌入式时钟只有几十 MHz 时才值得换成 FFT 卷积。第一版先跑通直接卷积再用已知信号验证后面再考虑优化这个顺序最不容易出错。4. 从 cwt.m 搬到 C 语言常见的五处差异和四个调试检查点4.1 cwt.m 和 C 实现的返回格式差异多数 cwt.m 脚本返回的矩阵是nsacles × N按行存尺度、按列存时间点。如果你的脚本最终要画 spectrogram那 C 实现导出的数据要按同样的行列顺序布局否则读回来图全反了。cwt.zip 这类压缩包里通常是一组相关脚本它们的矩阵方向、尺度排列方式未必相同。移植前先确认一件事原始 cwt.m 用的是wcoher风格接口还是自写卷积循环。前者按列存时间后者可能已经对核做了翻转。另一个常见差异是幅值归一化。MATLAB 的 cwt 函数默认有 L1 归一化某些旧版脚本和工具箱实现用 L2 归一化。C 实现里你要明确写一个归一化系数并在测试信号上和 cwt.m 对齐。最简单的方法是构造一个已知频率的余弦信号分别在 cwt.m 和 C 代码里跑比较同一位置系数实部和虚部的比例关系而不是直接比较绝对值。4.2 边界延拓的选择零填充、对称延拓、周期延拓直接卷积里if (idx 0 || idx n) continue;是零填充。零填充的缺点是信号两端会突然截断到零产生高频假成分这些假成分在时频图上表现为和时间轴平行的亮线。对称延拓在语音、振动信号里更常用它把边界之外的样本按镜面反射补齐int idx b k - half; if (idx 0) idx -idx - 1; if (idx n) idx 2 * n - idx - 1;这段替换的是内层循环里的越界判断。对称延拓假设信号在边界处是连续的但不是周期性的对非整周期截断的信号更友好。周期延拓适用于信号本身是周期的场景比如旋转机械的整周采集数据。移植 cwt.m 时先查原脚本用的延拓方法。MATLABcwt函数默认在内部做对称延拓如果你在 C 实现里用零填充开头和结尾各一段系数的幅值会系统性偏小而中间部分基本一致。调试时只看中间 60% 的数据大概率能快速定位是否只是边界策略差异。4.3 检查点一把核导出成文本对比调试 cwt C 实现时最有效的工具是文本导出。把某个固定尺度的核打印成两列表格一列是 tau一列是复数核的实部虚部然后在 MATLAB 里用相同的 w0、scale、fs 手动生成核对比差值。这一步能排除 90% 的公式写错问题。for (int i 0; i len; i) { fprintf(fp, %d %.9f %.9f\n, i, creal(psi[i]), cimag(psi[i])); }这里打印用定点十进制而不是%g是为了避免短浮点格式把细微的虚部舍入掉。对比时关注两个指标峰值位置应该处在len/2处幅值应该等于1/sqrt(scale) * pi^(-0.25)这两个都不满足就说明 tau 的计算或者归一化系数有问题。4.4 检查点二输出矩阵的行列方向和尺度排列C 语言里输出缓冲区只有一个维度一旦写成了out[b * nscales s]后面所有按列取数据的步骤都会出错。建议在 cwt_direct 函数入口处加一个断言式的自检// 用最小尺度跑一遍单点信号 double x[1] {1.0}; double complex y[1]; cwt_direct(x, 1, 1000.0, scales[0], 1, w0, y); double expect creal(y[0]);单点信号的 cwt 实际等价于该尺度下小波核在零点处共轭值的积分近似它的实部虚部可以直接手算验证。把这个测试写进回归用例比每次都重新对照 cwt.m 高效得多。5. cwt C 语言实现的参数调优ω0、尺度数量、浮点类型的选择5.1 中心频率 ω0 是分辨率的开关Morlet 小波的 ω0 取值直接决定时频平面的格子形状。ω0 越小比如 2小波振荡次数少时间分辨率好但频率分辨率差相邻频率成分容易糊在一起ω0 越大比如 10振荡次数多频率分辨率好但时间分辨率差瞬态跳变会被拉宽。ω0振荡周期数适用场景2 ~ 41 ~ 2 个周期故障脉冲、瞬态检测定位突变时刻6约 3 个周期通用分析平衡时频分辨率8 ~ 124 个以上周期谐波、窄带成分分离频率精度优先如果你只是把 cwt.m 的默认参数搬到 C 代码先固定 ω0 6不要一上来就微调。只有在时频图上发现相邻两个频率峰分不开或者峰值位置被拉宽到不能接受时再调大 ω0。反过来如果信号里关心的是突变时刻就调小 ω0。5.2 尺度数量的取舍太少锯齿太多冗余尺度数量 nscales 决定时频图的纵抽分辨率。太小时频率轴呈锯齿状峰值之间有明显棱角太多时计算量和内存线性增长但频率分辨率提升有限因为相邻尺度的核在频域上的重叠度已经很高。经验值是变频范围 2 Hz 到 150 Hz取 48 到 64 个尺度已经能在谱图上看到平滑的峰脊。如果需要做自动峰值提取而不是人眼观察32 个尺度足够。int nscales 48; double s_min w0 / (2.0 * M_PI * 150.0); double s_max w0 / (2.0 * M_PI * 2.0); double *scales malloc((size_t)nscales * sizeof(double)); log_spaced_scales(scales, nscales, s_min, s_max);这里 s_min 和 s_max 由频率范围反推。如果你要分析的信号频率范围是 10 到 400 Hz把两个 150 和 2 换成对应值即可。5.3 float、double、还是定点按平台选桌面平台直接用complex double最稳妥。嵌入式平台有时为了性能换成complex float但要注意Morlet 核里exp(-0.5 * tau * tau)在 tau 绝对值较大时接近零float 的尾数精度不足会造成核尾部出现不连续的微小台阶表现为频谱上有微弱的梳状噪声。在 Cortex-M 系列也没有硬件浮点单元时更常见的做法是把核预计算成整数定点数。这个过程通常固定尺度数量在 PC 上把每个尺度的核写入头文件目标端只做乘加累加不再运行时计算 exp 和 cos。代价是修改 ω0 或采样率后要重新生成头文件迭代变慢。// 定点表示示例Q15 格式乘累加用 int32 累加器 int16_t psi_q15[len]; // 从 double 转换得到 int32_t acc 0; for (int k 0; k len; k) { int32_t xq (int32_t)x_q15[idx]; acc xq * psi_q15[k]; }这个片段只画出了定点卷积的核心骨架实际使用时还要处理累加器的饱和、移位回缩以及边界延拓在定点下的映射。第一版先保证和 double 版本比对时误差在 1% 以内再谈吞吐。5.4 直接卷积和 FFT 卷积的真实切换阈值前面的复杂度表格只是理论值实际取舍还要看平台。在 PC 上复乘是单周期指令循环开销占比大直接卷积在 N4096、核长小于 512 时并不慢。在嵌入式处理器上如果编译器没有做循环展开内层循环的边界判断和索引计算可能占一半时间。我个人的经验阈值是当核长大于信号长度的一半时换 FFT 卷积。FFT 卷积把每个尺度转换到频域核和信号的卷积变成频域逐点相乘再逆变换回时域。边界效应依然存在但可以通过延拓后在频域相乘实现线性卷积用 FFTW 这类库实现时要注意输出长度是 N len - 1而不是 N。// FFT 卷积的伪代码骨架实际需要调用具体 FFT 库 // L next_pow2(n len - 1) // fft_signal FFT(x, L); 零填充到 L 长度 // fft_kernel FFT(psi, L) // y IFFT(fft_signal .* fft_kernel, L) // row[b] y[b half] * dt伪代码里y[b half]取的是线性卷积结果中有效区域的中段。用 FFT 卷积时核不随平移变化所以每个尺度只需要做一次 FFT一共 nscales 次。这比直接卷积在每个平移点都做内层累加要高效得多尤其在大尺度核很长时。6. 验证方法把 C 语言 cwt 的结果和 cwt.m 对齐6.1 用合成信号做基准写一个固定频率的合成信号比如 50 Hz 正弦叠加 120 Hz 正弦采样率 1000 Hz时长 1 秒。这个信号简单到根本没有争议cwt 的幅值峰值应该落在频率对应尺度的位置。在 C 代码里构造信号时要避免使用 rand用确定性的正弦函数保证每次运行完全一致。double x[1000]; for (int n 0; n 1000; n) { double t n / 1000.0; x[n] sin(2.0 * M_PI * 50.0 * t) 0.5 * sin(2.0 * M_PI * 120.0 * t); }然后调用 cwt_direct把结果写进二进制文件。关键是要把尺度序列、fs、实际核长度一起写进文件头否则无法确定哪一列对应哪个频率。6.2 从尺度反推频率并比对峰位置CWT 的尺度到频率的换算关系是 f w0 / (2πa)。以 ω06 为例尺度 0.0191 秒对应 50 Hz尺度 0.00796 秒对应 120 Hz。验证时扫一遍输出矩阵找每个尺度下幅值最大的时间点再把最大的尺度反推频率。理论上两个峰应该出现在这两行附近。如果峰值偏了一个尺度多半是尺度序列的起始点不对如果峰值偏了几个时间点多半是核的对称中心没有对齐len/2如果峰值频率对了但幅值差了一倍多半是归一化系数 L1/L2 的问题。6.3 把验证过程固化成脚本最后一件事也是我最推荐做的一步把验证过程固化成脚本用命令行的方式运行 C 程序输出结果再和参考值做比对。参考值不一定是 MATLAB 重新算一遍可以直接用上一次确认无误的 C 输出存成静态文件。对这个合成信号峰值位置和幅值都已经确定断言写成接近自然语言的形式也没有问题。./cwt_demo input.bin output.bin ./check_peak output.bin --freq 50 --tol 2这里--tol 2表示允许峰值位置偏差两个时间点。在 CI 或夜间构建里挂上这个断言以后改动核长度、尺度数量或延拓方式时回归不通过会先于你发现 cwt C 实现哪里被改坏了。本文还有配套的精品资源点击获取 C语言实现CWT小波变换：从离散公式到参数调优与MATLAB对标
简介一个基于C语言实现的小波变换程序包面向信号处理、图像分析领域的工程师与科研人员以多分辨率分析方式同时提取时间与频率特征弥补傅里叶变换丢失时域信息的不足。压缩包仅1个c文件大小约2KB核心代码cwt.c中提供输入信号sig、尺度向量scales、小波母函数wname等关键参数入口支持自定义尺度范围与小波类型便于移植到实时或离线任务中。已有464人学习下载。通过学习这份代码用户可以掌握连续小波变换的C语言实现思路理解尺度变换与卷积计算流程并以此为模板扩展更多小波基函数或集成至自己的信号处理项目。整体精简、无冗余文件适合具备一定C语言基础和小波理论知识的读者快速阅读与二次开发。1. 拿到 cwt.zip 却要在 C 语言里做小波变换先把 cwt.m 的思路放下手边有 cwt.zip、cwt.m甚至已经跑出来几张时频图却被要求“用 C 语言实现一遍 cwt”——这个需求在嵌入式信号处理里很常见。原因不外乎两种目标平台没有 MATLAB 运行时或者算法要进实时采集链路不能再靠脚本算完再画图。这里有个反直觉的结论cwt 的 C 实现本身不复杂复杂的是把尺度、采样率、小波中心频率这三个参数在代码里对齐否则同样的信号在 cwt.m 和 C 程序里出来的峰值位置、幅值比例会差很多。MATLAB 里 cwt 函数返回的是一个现成的复数矩阵换到 C 语言你得自己决定缓冲区怎么开、核函数怎么截断、边界怎么处理。这篇文章按“离散公式 - 最小可运行代码 - 参数调优 - 和 cwt.m 对齐”的顺序把 C 语言场景下实现小波变换的路线讲完整。2. cwt 的离散化把尺度 a、平移 b、采样率 fs 凑到同一个数组里2.1 从连续积分到数组求和三个替换必须同时做连续小波变换的定义式是W(a,b) 1 / sqrt(a) * ∫ x(t) * ψ*((t - b) / a) dt其中 ψ* 是小波母函数的共轭a 是尺度b 是平移。这个式子在 C 语言里不能直接求因为 t 是连续时间而信号 x 是一串按采样周期排列的离散点。把它变成程序能跑的求和要同时做三个替换第一积分换成求和dt 变成差分项第二小波函数的自变量 (t - b) / a 变成 (k - b) * T_s / a其中 T_s 1 / fs第三对每个平移 b只对落在小波核有效范围内的数据点累加超过核截断范围的部分不做贡献。三者缺一个结果就不是 cwt而是某种带通滤波或者卷积变体。C 语言里实现时我一般把尺度 a 的单位直接定成“秒”核函数的横轴就用采样点序号换算成秒。这样最小可运行版本里不需要引入复杂的单位换算逻辑。double dt 1.0 / fs; // fs 是采样率单位 Hz // 尺度 a 的单位是秒核自变量 t 的单位是秒 double t (double)(i - half) * dt / a;这里 scale 用秒而不是用“采样点数”原因是后续和频率换算时更直观。对于一个中心角频率为 w0 的 Morlet 小波它对应的中心频率是 w0 / (2 * PI * a)这里 a 用了秒频率直接就是 Hz。如果 a 换成采样点数公式里还得再乘一个 fs容易在移植 cwt.m 时对不上数。2.2 尺度序列怎么定对数间隔加最小尺度约束cwt 的尺度不是随便给几个数一般要按对数均匀分布铺开。原因是小波的尺度在低频端变化慢在高频端变化快线性间隔会让低频端分辨率不够高频端又冗余。常用做法是void log_spaced_scales(double *scales, int nscales, double s_min, double s_max) { double ratio pow(s_max / s_min, 1.0 / (nscales - 1)); scales[0] s_min; for (int i 1; i nscales; i) { scales[i] scales[i - 1] * ratio; } }这个函数做的事情很直接从最小尺度 s_min 开始每次乘一个固定比率直到 s_max。比率按 nscales 的个数均匀拆分。这样在双对数坐标下尺度点是一串等间隔的点跟 cwt.m 里大部分实现默认的排列方式一致。s_min 不能小于信号最高频率对应的半个周期。假设信号最高有效频率是 150 HzMorlet 用 w06那么最小尺度大约是 6 / (2 * PI * 150)约 0.0064 秒。如果再小小波核的振荡频率超过奈奎斯特或者超过信号自身带宽算出来的只有核的高频尾巴没有实际物理意义。2.3 为什么 cwt 的核用复数而不是实数用 Morlet 小波做 cwt核函数天然是复数的因为 Morlet 本质上是一个复指数乘高斯窗。实数核也可以做小波变换但它们通常是差分型小波比如 Mexican Hat得到的系数是实数频率选择性比复数核差。Morlet 核的表达式是ψ(τ) π^(-1/4) * exp(i * w0 * τ) * exp(-τ^2 / 2)其中 τ 是核的时间轴。在 C 语言里这个表达式对应两个组成部分复指数 e^(iw0τ) 决定振荡频率高斯窗 exp(-τ²/2) 决定时间局域性。w0 越大振荡越密集频率分辨率越好但时间分辨率越差。cwt 里对这个参数的处理直接决定了变换结果的质量这也是后文参数调优里要重点展开的内容。复数核带来的一个直接后果是C 语言实现时必须用complex double或手写复数结构体。如果只把核存成实部就丢了相位信息幅值计算也不准确。2.4 核截断高斯窗决定了核的有限长度理论上 Morlet 核在时间轴上无限延伸但高斯窗 exp(-τ²/2) 衰减很快。cwt 在实现时一般截断到 ±4 个标准差范围。τ 的标准差是 1所以截断到 τ ∈ [-4, 4] 时窗函数尾部只有 exp(-8) ≈ 0.000335完全可以忽略。每个尺度的核长度不是固定的它跟尺度成正比。尺度越大核在时间上拉得越宽覆盖的采样点数越多。这个特性让某些优化变得不直观固定 N 的情况下小尺度卷积很快大尺度卷积可能要遍历数倍于 N 的核长度。边界问题也因此在大尺度下更严重。下一节实现代码时核长度按4 * scale / dt来取再用奇数保证关于中心对称。3. 用 C 语言实现 cwt 的最小可用版本核缓存加卷积循环3.1 数据存储结果矩阵用一维数组核用 complex doubleC 语言没有 MATLAB 那种二维复数矩阵所以要自己定义内存布局。结果矩阵我一般开成一个一维数组长度为nscales * n第 s 个尺度、第 b 个时间点存放在out[s * n b]。这样不管之后是写文件、传参还是在 GUI 里显示都是一段连续内存方便用memcpy或直接fwrite导出也方便在 OpenMP 里按尺度拆分并行。double complex *out calloc((size_t)nscales * n, sizeof(double complex)); if (!out) { /* 处理分配失败 */ }分配时用calloc而不是malloc的好处是如果第一个尺度的卷积因为边界问题少算了几项未初始化区域直接是零不会出现随机复数噪声把调试过程带偏。3.2 生成 Morlet 共轭核一个尺度的完整核模板cwt 公式里用的是 ψ*即共轭。所以在生成核时直接生成共轭形式能省掉卷积循环里每一点都取共轭的操作。下面这个函数把某一尺度下的 Morlet 共轭核写入预分配缓冲区。static void morlet_conj_kernel(double complex *psi, int npoints, double scale, double w0, double dt) { int mid npoints / 2; double norm 1.0 / sqrt(scale) * pow(M_PI, -0.25); for (int i 0; i npoints; i) { double tau (double)(i - mid) * dt / scale; double env exp(-0.5 * tau * tau); double re env * cos(w0 * tau) * norm; double im -env * sin(w0 * tau) * norm; // 取共轭后的虚部 psi[i] re I * im; } }逻辑说明i - mid把缓冲区下标换算成相对中心的偏移量tau是偏移量对应的连续时间除以scale后代入母函数。env是高斯窗决定核的衰减cos是实部-sin是共轭后的虚部。norm是尺度归一化系数保证不同尺度下的能量一致。参数说明npoints必须是奇数否则核不对称卷积结果会产生一个采样点的相位偏移w0取 6 时Morlet 小波近似满足允许条件直接用于滤波不会引入明显的均值漂移scale是秒不是样本数这样tau永远是无量纲时间比较稳妥。3.3 直接卷积核心遍历尺度和采样点内层扫核void cwt_direct(const double *x, int n, double fs, const double *scales, int nscales, double w0, double complex *out) { double dt 1.0 / fs; size_t stride (size_t)n; for (int s 0; s nscales; s) { int len (int)(4.0 * scales[s] / dt) 1; if (len 5) len 5; if (len % 2 0) len; if (len 4 * n) len 4 * n; double complex *psi malloc((size_t)len * sizeof(double complex)); if (!psi) continue; morlet_conj_kernel(psi, len, scales[s], w0, dt); int half len / 2; double complex *row out (size_t)s * stride; for (int b 0; b n; b) { double complex acc 0.0 0.0 * I; for (int k 0; k len; k) { int idx b k - half; if (idx 0 || idx n) continue; acc x[idx] * psi[k]; } row[b] acc * dt; } free(psi); } }这段代码是最直接的 cwt C 语言实现没有任何优化但可读性最高适合作为移植 cwt.m 时核对结果的基准。外层循环按尺度展开中层循环按时间平移点展开内层循环把核和信号逐点乘加。边界处直接跳过越界索引等价于零填充延拓这是最简单但误差最大的方式。参数说明len取4 * scale / dt对应 τ ∈ [-4, 4]尾部误差在千分之一左右stride用n保证每个尺度在输出里是连续一块最后乘dt是因为离散积分必须带上采样间隔否则结果量纲不对核长度超过4*n时强制截断避免尺度太大导致核比信号还长内存和计算量失控。3.4 复杂度评估什么时候这个实现能直接上线对比项直接卷积FFT 卷积每个尺度复杂度O(N * L_k)O(N log N L_k log L_k)核远短于信号时最快无额外内存固定开销高不划算核与信号同量级时明显变慢优势明显额外内存只需一组核需要 FFT 缓冲区和频域掩模嵌入式定点平台舍入误差容易控制动态范围难控制当 N1024、尺度 30 个、平均核长 200 点时直接卷积大约做 614 万次复数乘加普通桌面 CPU 在毫秒级完成完全可以直接用。只有当 N 到几万、或嵌入式时钟只有几十 MHz 时才值得换成 FFT 卷积。第一版先跑通直接卷积再用已知信号验证后面再考虑优化这个顺序最不容易出错。4. 从 cwt.m 搬到 C 语言常见的五处差异和四个调试检查点4.1 cwt.m 和 C 实现的返回格式差异多数 cwt.m 脚本返回的矩阵是nsacles × N按行存尺度、按列存时间点。如果你的脚本最终要画 spectrogram那 C 实现导出的数据要按同样的行列顺序布局否则读回来图全反了。cwt.zip 这类压缩包里通常是一组相关脚本它们的矩阵方向、尺度排列方式未必相同。移植前先确认一件事原始 cwt.m 用的是wcoher风格接口还是自写卷积循环。前者按列存时间后者可能已经对核做了翻转。另一个常见差异是幅值归一化。MATLAB 的 cwt 函数默认有 L1 归一化某些旧版脚本和工具箱实现用 L2 归一化。C 实现里你要明确写一个归一化系数并在测试信号上和 cwt.m 对齐。最简单的方法是构造一个已知频率的余弦信号分别在 cwt.m 和 C 代码里跑比较同一位置系数实部和虚部的比例关系而不是直接比较绝对值。4.2 边界延拓的选择零填充、对称延拓、周期延拓直接卷积里if (idx 0 || idx n) continue;是零填充。零填充的缺点是信号两端会突然截断到零产生高频假成分这些假成分在时频图上表现为和时间轴平行的亮线。对称延拓在语音、振动信号里更常用它把边界之外的样本按镜面反射补齐int idx b k - half; if (idx 0) idx -idx - 1; if (idx n) idx 2 * n - idx - 1;这段替换的是内层循环里的越界判断。对称延拓假设信号在边界处是连续的但不是周期性的对非整周期截断的信号更友好。周期延拓适用于信号本身是周期的场景比如旋转机械的整周采集数据。移植 cwt.m 时先查原脚本用的延拓方法。MATLABcwt函数默认在内部做对称延拓如果你在 C 实现里用零填充开头和结尾各一段系数的幅值会系统性偏小而中间部分基本一致。调试时只看中间 60% 的数据大概率能快速定位是否只是边界策略差异。4.3 检查点一把核导出成文本对比调试 cwt C 实现时最有效的工具是文本导出。把某个固定尺度的核打印成两列表格一列是 tau一列是复数核的实部虚部然后在 MATLAB 里用相同的 w0、scale、fs 手动生成核对比差值。这一步能排除 90% 的公式写错问题。for (int i 0; i len; i) { fprintf(fp, %d %.9f %.9f\n, i, creal(psi[i]), cimag(psi[i])); }这里打印用定点十进制而不是%g是为了避免短浮点格式把细微的虚部舍入掉。对比时关注两个指标峰值位置应该处在len/2处幅值应该等于1/sqrt(scale) * pi^(-0.25)这两个都不满足就说明 tau 的计算或者归一化系数有问题。4.4 检查点二输出矩阵的行列方向和尺度排列C 语言里输出缓冲区只有一个维度一旦写成了out[b * nscales s]后面所有按列取数据的步骤都会出错。建议在 cwt_direct 函数入口处加一个断言式的自检// 用最小尺度跑一遍单点信号 double x[1] {1.0}; double complex y[1]; cwt_direct(x, 1, 1000.0, scales[0], 1, w0, y); double expect creal(y[0]);单点信号的 cwt 实际等价于该尺度下小波核在零点处共轭值的积分近似它的实部虚部可以直接手算验证。把这个测试写进回归用例比每次都重新对照 cwt.m 高效得多。5. cwt C 语言实现的参数调优ω0、尺度数量、浮点类型的选择5.1 中心频率 ω0 是分辨率的开关Morlet 小波的 ω0 取值直接决定时频平面的格子形状。ω0 越小比如 2小波振荡次数少时间分辨率好但频率分辨率差相邻频率成分容易糊在一起ω0 越大比如 10振荡次数多频率分辨率好但时间分辨率差瞬态跳变会被拉宽。ω0振荡周期数适用场景2 ~ 41 ~ 2 个周期故障脉冲、瞬态检测定位突变时刻6约 3 个周期通用分析平衡时频分辨率8 ~ 124 个以上周期谐波、窄带成分分离频率精度优先如果你只是把 cwt.m 的默认参数搬到 C 代码先固定 ω0 6不要一上来就微调。只有在时频图上发现相邻两个频率峰分不开或者峰值位置被拉宽到不能接受时再调大 ω0。反过来如果信号里关心的是突变时刻就调小 ω0。5.2 尺度数量的取舍太少锯齿太多冗余尺度数量 nscales 决定时频图的纵抽分辨率。太小时频率轴呈锯齿状峰值之间有明显棱角太多时计算量和内存线性增长但频率分辨率提升有限因为相邻尺度的核在频域上的重叠度已经很高。经验值是变频范围 2 Hz 到 150 Hz取 48 到 64 个尺度已经能在谱图上看到平滑的峰脊。如果需要做自动峰值提取而不是人眼观察32 个尺度足够。int nscales 48; double s_min w0 / (2.0 * M_PI * 150.0); double s_max w0 / (2.0 * M_PI * 2.0); double *scales malloc((size_t)nscales * sizeof(double)); log_spaced_scales(scales, nscales, s_min, s_max);这里 s_min 和 s_max 由频率范围反推。如果你要分析的信号频率范围是 10 到 400 Hz把两个 150 和 2 换成对应值即可。5.3 float、double、还是定点按平台选桌面平台直接用complex double最稳妥。嵌入式平台有时为了性能换成complex float但要注意Morlet 核里exp(-0.5 * tau * tau)在 tau 绝对值较大时接近零float 的尾数精度不足会造成核尾部出现不连续的微小台阶表现为频谱上有微弱的梳状噪声。在 Cortex-M 系列也没有硬件浮点单元时更常见的做法是把核预计算成整数定点数。这个过程通常固定尺度数量在 PC 上把每个尺度的核写入头文件目标端只做乘加累加不再运行时计算 exp 和 cos。代价是修改 ω0 或采样率后要重新生成头文件迭代变慢。// 定点表示示例Q15 格式乘累加用 int32 累加器 int16_t psi_q15[len]; // 从 double 转换得到 int32_t acc 0; for (int k 0; k len; k) { int32_t xq (int32_t)x_q15[idx]; acc xq * psi_q15[k]; }这个片段只画出了定点卷积的核心骨架实际使用时还要处理累加器的饱和、移位回缩以及边界延拓在定点下的映射。第一版先保证和 double 版本比对时误差在 1% 以内再谈吞吐。5.4 直接卷积和 FFT 卷积的真实切换阈值前面的复杂度表格只是理论值实际取舍还要看平台。在 PC 上复乘是单周期指令循环开销占比大直接卷积在 N4096、核长小于 512 时并不慢。在嵌入式处理器上如果编译器没有做循环展开内层循环的边界判断和索引计算可能占一半时间。我个人的经验阈值是当核长大于信号长度的一半时换 FFT 卷积。FFT 卷积把每个尺度转换到频域核和信号的卷积变成频域逐点相乘再逆变换回时域。边界效应依然存在但可以通过延拓后在频域相乘实现线性卷积用 FFTW 这类库实现时要注意输出长度是 N len - 1而不是 N。// FFT 卷积的伪代码骨架实际需要调用具体 FFT 库 // L next_pow2(n len - 1) // fft_signal FFT(x, L); 零填充到 L 长度 // fft_kernel FFT(psi, L) // y IFFT(fft_signal .* fft_kernel, L) // row[b] y[b half] * dt伪代码里y[b half]取的是线性卷积结果中有效区域的中段。用 FFT 卷积时核不随平移变化所以每个尺度只需要做一次 FFT一共 nscales 次。这比直接卷积在每个平移点都做内层累加要高效得多尤其在大尺度核很长时。6. 验证方法把 C 语言 cwt 的结果和 cwt.m 对齐6.1 用合成信号做基准写一个固定频率的合成信号比如 50 Hz 正弦叠加 120 Hz 正弦采样率 1000 Hz时长 1 秒。这个信号简单到根本没有争议cwt 的幅值峰值应该落在频率对应尺度的位置。在 C 代码里构造信号时要避免使用 rand用确定性的正弦函数保证每次运行完全一致。double x[1000]; for (int n 0; n 1000; n) { double t n / 1000.0; x[n] sin(2.0 * M_PI * 50.0 * t) 0.5 * sin(2.0 * M_PI * 120.0 * t); }然后调用 cwt_direct把结果写进二进制文件。关键是要把尺度序列、fs、实际核长度一起写进文件头否则无法确定哪一列对应哪个频率。6.2 从尺度反推频率并比对峰位置CWT 的尺度到频率的换算关系是 f w0 / (2πa)。以 ω06 为例尺度 0.0191 秒对应 50 Hz尺度 0.00796 秒对应 120 Hz。验证时扫一遍输出矩阵找每个尺度下幅值最大的时间点再把最大的尺度反推频率。理论上两个峰应该出现在这两行附近。如果峰值偏了一个尺度多半是尺度序列的起始点不对如果峰值偏了几个时间点多半是核的对称中心没有对齐len/2如果峰值频率对了但幅值差了一倍多半是归一化系数 L1/L2 的问题。6.3 把验证过程固化成脚本最后一件事也是我最推荐做的一步把验证过程固化成脚本用命令行的方式运行 C 程序输出结果再和参考值做比对。参考值不一定是 MATLAB 重新算一遍可以直接用上一次确认无误的 C 输出存成静态文件。对这个合成信号峰值位置和幅值都已经确定断言写成接近自然语言的形式也没有问题。./cwt_demo input.bin output.bin ./check_peak output.bin --freq 50 --tol 2这里--tol 2表示允许峰值位置偏差两个时间点。在 CI 或夜间构建里挂上这个断言以后改动核长度、尺度数量或延拓方式时回归不通过会先于你发现 cwt C 实现哪里被改坏了。本文还有配套的精品资源点击获取