连续小波变换C语言实现：从cwt.m到嵌入式信号处理与优化

简介：这是一个用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 语言实现：为什么选复数小波

母小波的选择直接决定系数含义。常见做法总结成表：

对于做 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.0，sigma 是高斯包络宽度，一般取 1.0。env 是高斯窗，它保证母小波在 t=0 附近以外快速衰减；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))，这样保证首尾精确落在上下限上，中间等比分布。

尺度数的选择没有绝对标准，但参考经验：

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.dat

gnuplot 里用以下命令生成热图，这条命令能快速检查主循环是否正确，不需要等 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)。当信号长度从几千点涨到几十万点，窗口长度又随尺度变大时，计算时间会很难看。

参考这笔账：

如果你只是拿 cwt.zip 里的代码做离线分析，1 万点以下直接法完全够用。但我一般会在代码里留一个开关：当n * nscales > 1e7时自动切换到 FFT 路径。FFT 在这里不改变结果，只是把每个尺度的卷积从时域线性内积变成频域点乘。

4.2 用卷积定理把每个尺度换成一次 FFT 与 IFFT

CWT 每个尺度的数学本质是信号与一个模板的卷积，卷积定理可以直接套用。流程：

• 对信号 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 influence，coi），落在该区域外的系数标注为不可信。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 默认输出更容易对齐。

输出时若想着重看动态范围，建议先转 dB：

double mag2 = re * re + im * im; double db = 10.0 * log10(mag2 + 1e-12);

参数说明：10 * log10用于功率，20 * log10用于幅度。我一般在存储时保留复数或幅度，只在可视化函数里转 dB。加 1e-12 防止纯零系数取对数时产生 -inf，这个细节在系数稀疏时特别重要。

综合起来的推荐参数组合：

5. 用合成信号校验 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, endpoint=False) 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 的导数，所以在 t=0 附近频率约为 5Hz，在 t=10 附近约为 100Hz。用这个信号跑 CWT，热图上应该看到一条从左下到右上的亮线，而不是水平直线，这能一举验证尺度轴、频率映射和窗口拉伸三件事是否同时正确。

5.2 从 CWT 输出矩阵自动回读峰值频率

CWT 的输出矩阵里，每列是某个时间点在各个尺度上的系数幅度。取某一列的最大值对应的尺度，就能算出该时刻的瞬时频率。命令行一条 awk 就能做峰值回读：

awk '{ for (j=1; j<=NF; j++) if ($j > max[j]) { max[j]=$j; f[j]=NR } } END { for (j=1; j<=NF; 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 实现就可以作为独立模块进入工具链。

本文还有配套的精品资源，点击获取