简介这是一个用C语言实现连续小波变换CWT的源码包适合信号处理初学者、嵌入式开发人员以及需要在C/C工程中集成时频分析功能的工程师。代码通过尺度向量与小波母函数参数对输入信号进行多分辨率分解在保留时间信息的同时提取频域特征弥补了传统傅里叶变换在时频联合分析上的不足。压缩包内共1个文件类型为c源代码整体大小仅2KB结构精简易读便于直接查看函数实现、分析算法流程并移植到其他项目。已有464人学习下载。对于想理解小波变换底层逻辑的读者这份源码提供了一条清晰的入门路径从输入信号sig、尺度数组scales到小波名称wname均可自行调整观察不同尺度和不同小波基下的变换输出实现中涉及的尺度分析和特征提取步骤能帮助理解信号局部特征捕捉方法。无论是用于课堂实验、算法对比还是作为实时信号处理模块的参考实现都具有不错的实用价值。1. CWT 的 C 语言实现把 cwt.m 从验证脚本变成可移植模块手头信号不是越大越好而是要同时看到时间位置和频率成分这件事STM32 上的一个振动采集板经常就能把工程师卡住FFT 给全局频率却看不出轴承故障发生在哪一段短时傅里叶窗口一固定低频分辨率又不够。这时候连续小波变换CWT几乎成了标配选项相关热度词里也总能看到小波变换图像增强、水印鲁棒性攻击、机械故障检测这些场景同时出现。麻烦在于团队里能跑的参考实现大都是 cwt.m一旦要把它移植到 C 语言环境就有两个拦路虎一是 MATLAB 内置的 cwt 既有工具箱又封装了边界处理和归一化代码里看到的只是接口二是网上流传的 cwt.zip 源码包往往写得能跑但不完整频率轴标定、边界裁剪、归一化约定这三件事都藏在细节里直接抄过来十有八九对不上数。这篇文章不是某份源码包的解析而是顺着「cwt.m 逻辑 → C 语言可移植实现」这条最常见路径讲清楚我一般怎么做从连续公式怎么离散化到尺度、母小波、卷积循环怎么写再到 FFT 加速、边界锥和频率校正最后用一段合成信号把 C 输出和 cwt.m 对拍。适合要把 CWT 嵌进 C/C 工具链的人也适合只想搞清楚 cwt.zip 里那些参数到底怎么设的读者。2. CWT 离散化尺度、平移、母小波在 C 语言里的表示2.1 连续 CWT 公式的三次替换以及 cwt.m 的频率标定逻辑连续小波变换的定义式是W(a,b) (1 / sqrt(|a|)) * ∫ x(t) * ψ*((t - b) / a) dt其中 a 是尺度b 是平移量ψ* 是母小波的共轭。要把它写成 C 代码必须做三次离散替换缺一不可。尺度把连续 a 替换成离散序列 a_i常用做法是对数均匀分布也就是让频率轴按倍频程或更细的等比分度。平移把连续 b 替换成采样点索引 j步长就是采样周期 dt。母小波把 ψ((t-b)/a) 替换成有限长离散模板不能对每个点都算一整条无限长的函数否则循环根本跑不动。这里还有一处与 cwt.m 对齐的关键点尺度 a 和物理频率 f 的关系。cwt.m 内部用「伪频率」做输出公式是f fc / (a * dt)其中 fc 是母小波的中心频率dt 是采样周期。所以反过来给定目标频率 f尺度就是a fc / (f * dt)我在 C 实现里不做「尺度到频率」的翻译而是直接让调用者传频率上下限内部用这个公式换算成尺度这样输出可以直接画在频率-时间图上。这就是 5.1 节找峰值能直接和理论频率对上号的原因。2.2 Morlet 母小波的 C 语言实现为什么选复数小波母小波的选择直接决定系数含义。常见做法总结成表母小波实/复数频率选择性实现成本典型用途Morlet复数好中心频率明确低一个复指数乘高斯包络时频分析、特征提取、图像增强中的频谱分量分离墨西哥帽Ricker实数一般旁瓣高低突变检测、边缘提取Haar实数差频域拖尾长最低快速变换、教学演示Morse复数可调参数多中MATLAB 新版 CWT 默认对于做 C 语言实现我建议先用 Morlet参数少频率轴直观而且是 cwt.m 老版本最常用的选择方便对拍验证。复数形式的好处是幅度包络不会因相位对齐问题出现零值做热力图和峰值检测都更稳。Morlet 的连续形式是高斯包络乘复指数写成 C 语言需要先定义一个最简单的复数结构体typedef struct { double re; double im; } cplx; cplx morlet_wavelet(double t, double fc, double sigma) { cplx w; double env exp(-0.5 * (t / sigma) * (t / sigma)); double arg 2.0 * M_PI * fc * t; w.re env * cos(arg); w.im env * sin(arg); return w; }参数说明t 是尺度化后的无量纲时间fc 是中心频率通常取 0.8125 或 1.0sigma 是高斯包络宽度一般取 1.0。env 是高斯窗它保证母小波在 t0 附近以外快速衰减arg 的系数 2π 让母小波在频域峰值刚好落在 fc 上。2.3 模板缓存每个尺度只算一次母小波如果每个平移点都重新调用 morlet_wavelet代价是毫无意义的重复计算。正确做法是先确定窗口宽度再生成一段一次性模板数组。Morlet 的有效支撑大约在 ±3σ 附近拉伸到尺度 a 后窗口半长取half ceil(3.0 * sigma * a)模板长度就是2 * half 1。实际代码里 sigma 和 fc 合并进尺度计算模板生成一次存成数组后续内积只做乘加。这一节解决的是 C 语言实现里最常见的问题算法思路对但每个系数都重新算指数函数导致规模稍大就跑不动。3. C 语言实现 CWT 主循环读文件、尺度扫描与系数输出3.1 单列信号读取文件读写操作里的两个防错点CWT 的输入一般是单列数值文件可能是 .txt、.csv 或采集卡直接导出的数据。读取这一段很多人习惯用 fscanf 循环到底但有两个防错点值得写进代码一是文件可能包含空行或不同分隔符fscanf 的 %lf 会跳过空白所以能用二是必须先数行数分配内存不能假定固定长度。下面是最小可用的读取函数double *read_signal(const char *path, int *n) { FILE *fp fopen(path, r); if (!fp) { perror(open); exit(1); } int cap 4096; double *x malloc(sizeof(double) * cap); int i 0; while (fscanf(fp, %lf, x[i]) 1) { i; if (i cap) { cap * 2; x realloc(x, sizeof(double) * cap); } } fclose(fp); *n i; return x; }逻辑说明第一遍边读边扩容避免预先统计行数造成两遍读取fscanf 返回 1 表示成功读到一个浮点数遇到文件尾返回 EOF 退出。参数 path 是文件路径n 通过指针带回有效采样点数。注意扩容系数取 2这是 realloc 常见的折中——太小会频繁拷贝太大会浪费内存。3.2 对数频率轴生成与尺度换算CWT 的尺度轴不应该线性扫否则低频段分辨率被浪费高频段过密。我一般让调用者指定频率上下限和尺度数内部生成对数均匀序列double *logspace_freq(double fmin, double fmax, int nscales) { double *f malloc(sizeof(double) * nscales); double ratio fmin / fmax; for (int i 0; i nscales; i) { f[i] fmax * pow(ratio, ((double)i) / (nscales - 1)); } return f; }参数说明fmin 是要分析的最低频率fmax 是最高频率nscales 是尺度个数。第 i 个频率的实际公式是fmax * (fmin/fmax)^(i/(nscales-1))这样保证首尾精确落在上下限上中间等比分布。尺度数的选择没有绝对标准但参考经验参数典型值说明fmin信号基频的 0.1 倍太低导致窗口过长边界浪费严重fmaxmin(采样率/2, 目标最高频)不要超过奈奎斯特频率nscales32 到 128分析用 64出出版级图片用 128采样率 fs由硬件决定影响尺度-频率换算必须正确传入3.3 主循环核心按尺度滑窗、按平移点做内积现在进入 CWT 主循环。核心思路是对每个尺度生成一段 Morlet 模板然后在信号上按步长 1 滑动做内积。下面的代码是完整可运行的内核void cwt_core(const double *x, int n, double fs, double fmin, double fmax, int nscales, double *out) { double dt 1.0 / fs; double *freqs logspace_freq(fmin, fmax, nscales); for (int i 0; i nscales; i) { double a 0.8125 / (freqs[i] * dt); // 尺度单位是采样点 int half (int)ceil(3.0 * a); // Morlet 窗口半长 double norm dt / sqrt(a); // L2 归一化因子 for (int j 0; j n; j) { int left j - half; int right j half; if (left 0 || right n) { // 边界先置零 out[i * n j] 0.0; continue; } double re 0.0, im 0.0; for (int k left; k right; k) { double t ((double)(k - j) * dt) / a; cplx psi morlet_wavelet(t, 0.8125, 1.0); re x[k] * psi.re * dt; im - x[k] * psi.im * dt; } out[i * n j] norm * sqrt(re * re im * im); } } free(freqs); }逻辑说明外层 i 循环遍历频率轴内层 j 循环遍历平移点最内层 k 循环完成小波模板与信号的加权重叠积分。窗口随尺度自动变长高频时 half 小、计算快低频时 half 大、做得多这正好对应小波变换「低频看粒度、高频看细节」的特性。边界处先把系数置 0第 4.3 节再处理为更严谨的锥形裁剪。参数说明a 的计算用了前文公式的反推把频率换算成采样点数norm 里的 1/sqrt(a) 是 L2 归一化保证不同尺度下同幅值的正弦分量产生可比的系数这一点很多简版实现会漏掉导致高频分量看起来被放大。复数小波的共轭体现在im - x[k] * psi.im实数信号与小波实部做普通乘加、与虚部做符号翻转积分。3.4 输出与快速可视化用 gnuplot 检查热图输出直接写成一个二维矩阵每行对应一个尺度每列对应一个时间点。最简单的落地方式./cwt_demo signal.txt 1000 2 200 64 coeff.datgnuplot 里用以下命令生成热图这条命令能快速检查主循环是否正确不需要等 Python 环境gnuplot -e set view map; set pm3d; set logscale y; plot coeff.dat matrix using 1:2:3 with image观察步骤先看是否有明显的水平条带再看条带中心是否出现在你期望的信号频率上。如果只有一条水平亮线且横跨全图说明这是一个平稳正弦分量如果亮线随时间是弯曲的说明存在扫频分量。这个初步判断在引入任何优化之前都值得先做一次。4. CWT 的 C 语言优化FFT 加速、边界锥与频率轴标定4.1 直接卷积什么时候不可用复杂度与工程经验第 3 章的滑动内积实现复杂度是 O(nscales × n × half)。当信号长度从几千点涨到几十万点窗口长度又随尺度变大时计算时间会很难看。参考这笔账信号长度 n尺度数平均复杂度工作量直接法表现1k ~ 4k32约 10^6 量级毫秒级直接法足够10k ~ 50k64约 10^7 ~ 10^8 量级秒级勉强可接受100k 以上128约 10^9 量级不可接受必须换 FFT如果你只是拿 cwt.zip 里的代码做离线分析1 万点以下直接法完全够用。但我一般会在代码里留一个开关当n * nscales 1e7时自动切换到 FFT 路径。FFT 在这里不改变结果只是把每个尺度的卷积从时域线性内积变成频域点乘。4.2 用卷积定理把每个尺度换成一次 FFT 与 IFFTCWT 每个尺度的数学本质是信号与一个模板的卷积卷积定理可以直接套用。流程对信号 x 做一次 FFT得到 X。对当前尺度的母小波模板做 FFT得到 PSI。频域点乘 Y X * conj(PSI)再做 IFFT就得到该尺度的全部平移系数。伪代码轮廓// 预处理一次信号 FFT complex *X fft(x, n); // 每个尺度独立处理 for (int i 0; i nscales; i) { complex *psi make_wavelet_template(a, 2 * half 1); // 补零到 n complex *PSI fft(psi, n); for (int k 0; k n; k) { Y[k] X[k] * conj(PSI[k]); } complex *y ifft(Y, n); // 取实部归一化后写入对应行 }逻辑说明一次 FFT 只需要做一次每个尺度的模板 FFT 各做一次整体复杂度从 O(nscales × n × half) 降到 O(nscales × n log n)。这里的归一化因子 dt/sqrt(a) 要在 IFFT 后统一乘避免在频域混入尺度变化引起数值漂移。注意模板长度和信号长度不一致时模板要按自然位置左对齐、右侧补零到 n否则频域乘法会包含循环移位效应。第 4 章不建议提供完整 FFT 实现因为 radix-2 的代码在任何项目里都能找到。实际工程里常见做法是项目原本就有 FFT 库就直接复用没有的话先跑直接法再在确有性能瓶颈时引入避免第一版就被傅里叶变换的边界条件干扰主流程调试。4.3 边界效应与 cwt.m 的 coi 处理第 3.3 节把边界系数直接置 0这会导致一个假象热图在左右边缘突然出现一条暗带。如果不关心边缘这样做能看但如果做特征提取风险就大了——边缘处的小波窗口有一部分落在信号外系数明显偏小与真实的低幅值信号无法区分。cwt.m 的做法是标记锥形影响区cone of influencecoi落在该区域外的系数标注为不可信。C 语言实现里不需要画锥但要在输出矩阵中同步标记double coi a * 3.0; // 该尺度的边界可信半宽 if (j coi || (n - j - 1) coi) { out[i * n j] NAN; // 或传递一个 mask 数组 }参数说明3.0 对应 Morlet 高斯包络的三倍标准差窗口外侧的母小波幅度已经衰减到接近于零所以低于这个距离的位置存在明显能量截断。用 NAN 标记后绘图和统计都要跳过非数否则峰值检测会误判。4.4 幅度归一化与 dB 显示防止峰值读数漂移CWT 系数的幅值不只取决于信号强度还取决于归一化约定。直接法和 FFT 法很容易因为少了某个因子导致结果整体偏大或偏小。我在实现里统一采用 L2 归一化每个系数除以 sqrt(a)离散化时乘 dt。这样有三个好处同幅值不同频率的正弦波在各自中心尺度处峰值一致与 Parseval 定理兼容与 cwt.m 默认输出更容易对齐。输出时若想着重看动态范围建议先转 dBdouble mag2 re * re im * im; double db 10.0 * log10(mag2 1e-12);参数说明10 * log10用于功率20 * log10用于幅度。我一般在存储时保留复数或幅度只在可视化函数里转 dB。加 1e-12 防止纯零系数取对数时产生 -inf这个细节在系数稀疏时特别重要。综合起来的推荐参数组合使用场景信号长度fmin / fmaxnscales边界处理快速验证 4k基频 ~ Nyquist32置 0 即可特征提取10k ~ 50k关注频带上下限64输出 coi mask图像增强/水印大矩阵拉成行按行频谱峰值128FFT 路径 coi5. 用合成信号校验 CWT峰值回读频率并与 cwt.m 对拍5.1 构造一段非平稳信号作为验证基准验证 CWT 实现是否正确最可靠的方式不是看热图像不漂亮而是构造一个频率成分已知的非平稳信号让 CWT 输出反推出频率值。下面用 Python 生成 10 秒、1kHz 采样率、5Hz 到 100Hz 线性扫频的信号写入单列文本import numpy as np fs 1000 t np.linspace(0, 10, 10 * fs, endpointFalse) f0, f1 5.0, 100.0 phase 2 * np.pi * (f0 * t (f1 - f0) * t**2 / (2 * t[-1])) x np.sin(phase) np.savetxt(chirp.txt, x, fmt%.8f)逻辑说明瞬时频率是 phase 对 t 的导数所以在 t0 附近频率约为 5Hz在 t10 附近约为 100Hz。用这个信号跑 CWT热图上应该看到一条从左下到右上的亮线而不是水平直线这能一举验证尺度轴、频率映射和窗口拉伸三件事是否同时正确。5.2 从 CWT 输出矩阵自动回读峰值频率CWT 的输出矩阵里每列是某个时间点在各个尺度上的系数幅度。取某一列的最大值对应的尺度就能算出该时刻的瞬时频率。命令行一条 awk 就能做峰值回读awk { for (j1; jNF; j) if ($j max[j]) { max[j]$j; f[j]NR } } END { for (j1; jNF; j) printf %d %d\n, j, f[j] } coeff.dat因为输出行是按照频率从高到低排列的这里的 NR 实际代表尺度索引再结合生成的频率轴就能换算成物理频率。理论值与回读值的误差来源主要有两个一是对数频率轴的离散间隔二是 Morlet 的有限带宽让峰值出现在最接近真实频率的离散频点上。误差有一个频点以内实现就是可靠的。5.3 与 cwt.m 对拍统一采样、归一化、边界如果手上还有 MATLAB 环境最后一步对照会非常直接。核心步骤把 C 输出和 cwt.m 输出都存成矩阵C 这边每行一个尺度MATLAB 那边用writematrix导出。新版 MATLAB 使用cwt(x)时默认选 Morse 小波不能直接和 Morlet 的 C 实现比数值要求一致需要指定类似cwt(x, morl)的方式或把 C 实现改成与默认小波对齐。比较时不要逐点相比而要比较每列峰值位置的尺度索引。功率归一化约定不同会导致数值整体偏移但峰值位置应一致。我一般先打印两个矩阵在若干时间点上的峰值索引画在一张图上。如果两者在窗口中部完全重合只在边缘出现偏差那就说明边界处理方式不影响主体结果。判据很简单中心频点的偏差在半频点之内系数形状在 coi 内一致这套 CWT 的 C 实现就可以作为独立模块进入工具链。本文还有配套的精品资源点击获取