简介面向小波变换入门者的MATLAB连续小波变换快捷程序适合正在进行信号处理或图像处理课程学习的学生也适合刚接触小波理论的科研人员快速验证想法。程序聚焦连续小波变换的核心实现压缩包仅由1个.m脚本组成大小仅795B结构单纯便于逐行研读省去了复杂工程配置的干扰。截至目前已有151人学习说明其在同类入门资源中受到一定关注。下载并运行该脚本后读者能直观看到连续小波变换的编程步骤包括小波基选择、尺度与平移参数设置、系数计算与结果展示有助于将教材公式转化为可执行的代码认知。无论是快速验证信号分解效果还是对比不同小波基函数的差异这个脚本都能提供直观的实验载体可作为课程报告、结课设计或课题探索的基线程序。1. 从 wavelet.rar 到连续小波变换先搞懂它在算什么如果你刚下载了一个叫 wavelet.rar 的压缩包解压后大概率会看到一堆 .m 或 .py 文件、几张像频谱图的 png以及一个写着“连续小波变换”的 README。最容易卡住的不是调库而是不知道程序里那个三维图到底在算什么。连续小波变换CWT把一维信号投影到一组由小波基伸缩平移得到的函数上输出随尺度 a 和时间 b 变化的系数矩阵。它和短时傅里叶变换最大的区别是窗口自动随频率变化低频窗宽、高频窗窄正好匹配非平稳信号对时间分辨率的需求。适合做振动故障诊断、脑电分析、语音处理也常被用来给图像增强做多尺度分解打基础。下面沿着“公式→代码→调参→排错”的顺序把连续小波变换从能查到定义到能跑出靠谱谱图的过程讲清楚。回头再看 wavelet.rar 里的代码就不晕了。2. 连续小波变换的公式拆解与变尺度窗口的必然性2.1 CWT 公式里每个符号在代码里对应什么连续小波变换的常见写法是$$X_w(a,b)\frac{1}{\sqrt{|a|}}\int_{-\infty}^{\infty} x(t)\psi^*\left(\frac{t-b}{a}\right)dt$$这里的 $a$ 叫尺度$b$ 叫平移$\psi(t)$ 是小波基函数$*$ 表示复共轭。在 PyWavelets 里一行pywt.cwt(data, scales, wavelet, sampling_period)几乎把这个公式原样搬进了参数data对应 $x(t)$scales数组里的每个值就是一个 $a$wavelet对应 $\psi(t)$sampling_period负责把离散采样点索引换算成物理时间。返回的coefs[i][j]就是第 $i$ 个尺度、第 $j$ 个采样时刻的复系数。看公式时先抓住三个点尺度小对应高频小波基被压缩尺度大对应低频小波基被拉伸。$b$ 只负责平移不改变小波基形状。$\psi(t)$ 的均值必须为零否则积分在无穷远处不收敛这被称为容许条件admissibility condition。Morlet、墨西哥帽、高斯小波都满足这个条件这也是为什么它们能直接用作 CWT 小波基。最小验证代码import numpy as np import pywt fs 1000.0 t np.linspace(0, 1, 1000, endpointFalse) x np.sin(2 * np.pi * 50 * t) 0.5 * np.sin(2 * np.pi * 120 * t) scales np.arange(1, 100, 1) coefs, freqs pywt.cwt(x, scales, morl, sampling_period1/fs) print(coefs.shape) # (99, 1000) print(freqs[:3]) # 最高频的3个频率值coefs.shape的第一维等于len(scales)第二维等于信号长度。freqs的长度也等于len(scales)但顺序是高频到低频所以freqs[0]对应第一个尺度也就是最小时尺度。这里故意用np.arange(1,100,1)做演示实际工程里应该用对数间隔的尺度序列后面第 3 章会说明。2.2 为什么连续小波变换不是“换了个窗”的短时傅里叶短时傅里叶变换STFT用固定长度窗函数截信号。窗长定了时间分辨率和频率分辨率的乘积就受海森堡不确定性原理限制想提高时间分辨率就必须牺牲频率分辨率反之亦然。连续小波变换用 $a$ 去拉伸或压缩 $\psi(t)$分析高频时等效窗很短能定位突变的时刻分析低频时等效窗很长能区分接近的频率成分。这种变尺度窗口正好匹配语音、振动、生物电信号里常见的高频瞬态叠加低频趋势的结构。变换窗口特性时间分辨率频率分辨率典型场景傅里叶变换全局窗无最高平稳信号短时傅里叶变换固定窗恒定恒定缓变信号连续小波变换变尺度窗高频好低频好突变、多尺度信号图像增强里反复提到的“多尺度”也来自这个特性。图像的边缘、纹理在不同尺度下呈现不同粗细单一固定窗滤波很难同时保留细纹理和粗轮廓。小波族天然把尺度维保留下来后续增强可以按尺度分别处理。用 Python 做图像增强时工程上更常用二维离散小波变换wavedec2因为它计算快、可逆、冗余低但选小波基、控制分解层数的逻辑与 CWT 的尺度选择完全同源。2.3 拿到 wavelet.rar 后先别跑主程序先做三件事我拿到这类资料包不会直接双击运行主脚本。里面常有数据也常有画图代码但参数往往是写死的。常见做法是先做三件事确认信号采样率确认数据是否包含 NaN 或直流漂移用小波变换自带的自测函数验证频率轴。这里给一个简单的自测函数def inspect_signal(filename, fs): 读入信号并检查基本质量返回去均值后的数据 x np.loadtxt(filename, delimiter,) # 具体读取按文件格式调整 x np.nan_to_num(x) x x - np.mean(x) # CWT 对直流不敏感去均值可以减少边界伪影 print(flength{len(x)}, fs{fs}, duration{len(x)/fs:.2f}s) return x注意np.loadtxt只适用纯数值文本格式。如果wavelet.rar里是.mat文件就要用scipy.io.loadmat如果是.xlsx用pandas.read_excel。不要盲目照抄读取代码。去均值这步很关键因为连续小波变换本质是通过积分提取信号能量直流分量会造成尺度轴两端出现明显的水平亮带干扰频率峰值的判断。3. 用 Python 和 PyWavelets 跑通连续小波变换的最小闭环3.1 最小可运行版本一条信号、一组尺度、一张图在确认数据质量后最小闭环包含四个动作构造测试信号、选尺度序列、调用pywt.cwt、画出模值谱图。下面是完整代码import numpy as np import pywt import matplotlib.pyplot as plt fs 1000.0 t np.arange(0, 2, 1/fs) # 频率从 100 Hz 开始衰减的振荡信号模拟冲击响应 x np.sin(2 * np.pi * 100 * t) * np.exp(-2 * t) scales np.geomspace(5, 200, 64) # 对数尺度 coefs, freqs pywt.cwt(x, scales, cmor1.5-1.0, sampling_period1/fs) plt.figure(figsize(8, 5)) plt.pcolormesh(t, freqs, np.abs(coefs), shadingauto, cmapjet) plt.yscale(log) plt.xlabel(Time (s)) plt.ylabel(Frequency (Hz)) plt.colorbar(labelMagnitude) plt.show()代码里的np.geomspace(5, 200, 64)是尺度范围的工程默认写法。尺度序列如果太密相邻频率差别不明显图像上会产生横向过度平滑太疏又会漏掉窄带信号。64 个尺度用于预览正式分析用 128。plt.pcolormesh会自动处理不均匀的频率轴比imshow少踩坐标顺序的坑。3.2 scales 和 freqs 的换算写给不懂伪频率的人pywt.cwt在没传sampling_period时freqs返回的不是物理频率而是以采样间隔为单位的“伪频率”。很多教程省略这个参数导致画出来的频率轴数值与信号实际频率对不上。正确的换算关系是物理频率 小波中心频率 / (尺度 × 采样周期)可以用pywt.scale2frequency验证fc 1.0 # cmor1.5-1.0 中的中心频率 freq_expected pywt.scale2frequency(cmor1.5-1.0, scales) / (1/fs) print(np.max(np.abs(freqs - freq_expected)))这个差值正常情况下小于1e-12。如果你用的是morl它的中心频率是 0.8125 而不是 1.0没注意这点就会把所有频率成比例算偏。所以做定量分析时优先选择中心频率明确的小波比如cmorB-F格式其中B是带宽F是中心频率。3.3 “小波变换图像增强python”里为什么更常看到 wavedec2连续小波变换对一维信号的时频分析效果直观但二维图像直接做 CWT 计算量非常大而且连续小波基之间有高冗余。实际搜“小波变换图像增强python”时绝大多数可运行的代码使用的是二维离散小波变换 DWT。它本质上是把 CWT 的尺度按 2 的幂次离散化用滤波器组实现快速分解与重建。import pywt import numpy as np from PIL import Image img np.array(Image.open(lena.png).convert(L)).astype(float) coeffs pywt.wavedec2(img, db4, level3) cA3, (cH3, cV3, cD3), (cH2, cV2, cD2), (cH1, cV1, cD1) coeffs threshold 0.1 * np.max(np.abs(cH1)) for cH in (cH3, cH2, cH1): cH[cH] np.sign(cH[cH]) * np.maximum(np.abs(cH[cH]) - threshold, 0) img_enhanced pywt.waverec2(coeffs, db4) img_enhanced np.clip(img_enhanced, 0, 255).astype(np.uint8)说明wavedec2里level3对应三个尺度。cA3是低频近似cH/cV/cD分别是水平、垂直、对角高频细节。这里的阈值门限用的是全局阈值 0.1 倍最大值实际项目中要按不同层单独设置因为前面尺度的噪声能量分布不一样。db4的消失矩是 4对图像边缘有较好的稀疏表示。如果你手里wavelet.rar里的 CWT 代码无法直接处理图像可以先对图像的每一行做 CWT再按行拼接出二维时频图但更标准的路径还是切到 DWT。3.4 什么时候坚持 CWT什么时候换成 DWT任务建议理由观测瞬时频率随时间变化CWT变分辨率和相位信息图像降噪/增强DWT可逆、速度快、冗余低故障特征频率提取CWT频率轴连续便于峰值搜索CWT 适合分析长度几百到几百万点的一维信号图像这类二维数据即使只做边缘增强也建议先 DWT 粗分解。这样代码性能和结果可复现性都比强行 CWT 好。4. 连续小波变换的 5 个必调参数与实际选参顺序4.1 小波基的参数格式从morl到cmorB-Fpywt.cwt的wavelet参数可以直接传字符串比如morl、cmor1.5-1.0、gaus8。字符串里的数字是有含义的cmor1.5-1.0表示复数 Morlet 小波带宽 1.5中心频率 1.0gaus8表示高斯小波的 8 阶导数。选型没有绝对标准但常见经验是小波参数特点建议场景morl无实数、速度快、中心频率 0.8125快速预览cmorB-FB 带宽F 中心频率复数有相位频率轴直观时频定量分析gausPP 阶导多阶导数适合边缘检测信号突变点定位dbNN 消失矩正交或双正交适合离散图像/压缩带宽 B 和中心频率 F 一旦变化同一个尺度数组对应的物理频率就变了。设置cmor时带宽太小时域支撑长边界效应会占掉很大面积带宽太大则频率分辨率粗糙。我通常从cmor1.5-1.0开始再根据谱图效果微调。4.2 尺度范围的计算目标是覆盖目标频段尺度范围和采样率、小波中心频率的关系scales fc × fs / f其中 fs 是采样率f 是希望分析的物理频率。所以想分析 5 Hz 到 200 Hz、采样率 1000 Hz、中心频率 1.0 的信号尺度下限是 1×1000/2005尺度上限是 1×1000/5200。这个简单换算能避免盲目把尺度取到 1000 以上。用代码落地fc 1.0 fs 1000.0 fmin, fmax 5.0, 200.0 scale_min fc * fs / fmax # 5.0 scale_max fc * fs / fmin # 200.0 scales np.geomspace(scale_min, scale_max, 128) print(scale_min, scale_max)注意fmax不要超过奈奎斯特频率fs/2。超过之后采样数据里已经没有有效能量只会增加高频边界伪影。如果信号有 50 Hz 工频干扰可以按需要把频段下限定在 60 Hz 以上而不是把宽带噪声一起纳入。4.3 尺度点数不是越密越好尺度序列点数是 CWT 唯一的“分辨率旋钮”。点数越多coefs矩阵行数越多计算时间和内存线性上升。对 1 秒采样率 1000 的信号128 个尺度运行时间通常在几十毫秒到几百毫秒之间取 512 个尺度图像会显得平滑但不会产生新信息。经验值预览 64验证参数 128做峰值搜索 256。如果信号本身只有一个窄带频率64 点足够如果存在两个频率接近的振荡尺度点数要足以让它们在频率轴上分开此时可以用解析信号来检验分辨率而不是盲目增加点数。4.4 边界效应和去趋势最容易让谱图骗人的两个点CWT 在信号两端会引入边界伪影因为小波支撑可能超出数据范围。常见的缓解方法是先对信号做对称延拓modeperiodization或symmetric或者分析时把两端各截掉最小尺度对应的小波支撑长度。下面是带边界处理的写法coefs, freqs pywt.cwt( x, scales, cmor1.5-1.0, sampling_period1/fs, methodfft, # FFT 方式计算速度更快 modeperiodization ) # 丢弃首尾各 2% 的数据点避免边界宽带影响 trim int(len(t) * 0.02) coefs coefs[:, trim:-trim]methodfft适合对整段信号做连续分析因为 FFT 假定信号周期延拓所以要求mode也匹配周期延拓。如果信号不是严格周期改用methodconv会稍微慢一些但边界响应更接近直接积分。参数调整后要用下一节的自测信号确认频率峰位置不能让边界亮带影响判断。5. 用“已知频率信号”校验 wavelet.rar 里的 CWT 参数不管从 wavelet.rar 里拿到的是.mat、.csv还是.npy数据落地的第一步应该是构造一个已知频率的仿真信号跑同一套 CWT 参数检查频谱峰值频率是否准确。这个自检能一次性暴露采样率设置、尺度范围、小波中心频率三类错误。def validate_cwt(scales, wavelet, fs1000.0, f050.0): t np.arange(0, 1, 1/fs) x np.sin(2 * np.pi * f0 * t) coefs, freqs pywt.cwt(x, scales, wavelet, sampling_period1/fs) ridge_idx np.argmax(np.abs(coefs), axis0) ridge_freq freqs[ridge_idx] detected np.median(ridge_freq) rel_error abs(detected - f0) / f0 print(fdetected{detected:.2f} Hz, error{rel_error:.2%}) return coefs, freqs如果rel_error超过 5%先检查scales是否覆盖了f0再确认wavelet字符串里有没有把带宽和中心频率写反。实际中常见的错误是cmor1.5-1.0写成cmor1-1.5结果中心频率变成 1.5所有频率读数等比偏移 50%。请对比公式验证fc / (scale * Ts)中使用的 fc。验证通过后再把同一组参数应用到wavelet.rar提供的真实数据。如果真实数据谱图在预期频段之外还有持续亮带先用scipy.signal.detrend去掉线性趋势再看窗函数mode是否用了periodization。如果真实数据和仿真信号的采样率不同不要直接复用 scales先按 4.2 的公式重算。pywt.cwt不会因为你传了sampling_period就自动把 scales 缩放到目标频段它只会原样使用你传入的尺度数组。所以每次换数据先执行scale_min fc * fs / fmax和scale_max fc * fs / fmin这两行代码再生成新的np.geomspace。检查完频率轴再去看相位或包络就不会被彩色图上的边界噪声干扰了。本文还有配套的精品资源点击获取