四元数小波QWT：Python从零实现与纹理分类实战

简介四元数小波是一种将四元数代数与小波分析相结合的高级数学工具既保留小波在时频域的局部化分析能力又能利用四元数有效表示信号的旋转与对称特性在图像处理、信号去噪和通信领域中具有重要应用。资源包以hingfui.m为唯一源码文件完整实现了四元数信号构建、常用小波基选择、四元数小波分解、多尺度系数分析以及信号逆变换重构等核心流程可帮助读者快速建立从理论到代码的映射。整个压缩包体积约9KB仅含1个m文件代码精简、结构清晰适合熟悉MATLAB基础、想要了解或验证四元数小波算法的研究人员、工程师及高年级学生参考。目前已有235人学习下载可作为轻量级参考脚本用于信道编码误码检测、调制方式识别、信道频率选择性衰落补偿等具体问题的算法验证与二次开发对深入理解四元数小波的实际应用具有不错的参考价值。1. 四元数小波为什么二维信号需要三个相位普通 DWT 系数是实数回答的是“这个位置有没有结构、结构有多强”但如果问题变成“结构相对上一帧往哪个方向移动了多远”实数系数给不出直接答案。傅里叶变换有相位相位能测整体位移却把空间位置摊平。四元数小波Quaternion Wavelet TransformQWT把每个小波系数从实数升级成一个四元数——一个幅值加三个相位——既保留空间定位又把局部结构的水平、垂直、对角位移信息编码进相位里。做纹理分类、图像配准、光流和质量评价的人经常会遇到“DWT 不够用、FT 又太全局”的尴尬QWT 是解决这类问题最常用的工具箱。这篇文章不引入商业化库用 PyWavelets 和 SciPy 从零搭一套可运行的 QWT 分解再落到纹理分类和平移验证整个路径可以直接复现。2. 四元数小波的结构原理从希尔伯特变换到四通道滤波2.1 为什么 DWT 缺相位、傅里叶变换缺局部定位DWT 的实数系数来自与实小波基的内积只有幅值信息。对一个正弦条纹DWT 系数会告诉你这个位置有能量但无法区分它是往左移了 0.3 像素还是 0.8 像素。傅里叶变换的相位可以测全局平移但它的基函数无限延伸任何局部形变都会污染全部相位。QWT 是个折中基函数依然有紧支撑但系数是四元数通过三个相位把“水平位移”、“垂直位移”和“对角协调性”分离出来。变换系数类型相位信息空间定位典型任务DWT实数无好去噪、压缩DFT复数全局一个相位无全局配准QWT四元数三个局部相位好纹理、局部配准这里还需要回答一个前置问题为什么不用一维复小波CWT直接扩展复小波把实小波变成 ψ iψ_h系数是复数自带一个相位适合一维信号的瞬时频率和包络分析。图像是二维的一个相位无法同时区分水平与垂直两个方向的结构走向。把行、列两个方向的解析滤波结果放到四元数的三个虚部里才能用一个数同时表达水平位移、垂直位移和对角协调性三个自由度。2.2 可分离 QWT四个分量各司其职QWT 的常见构造方式是“可分离”的把一维实小波 ψ 和它的希尔伯特变换 ψ_h 分别放到行、列方向组合得到四个滤波通道。对图像 I 做分解后每个位置得到一个四元数系数 q a bi cj dk其中 a 是实部b、c、d 是三个虚部。分量行方向滤波器列方向滤波器主要编码aψ_hψ_h近似能量与整体相移bψψ_h水平方向的结构位移cψ_hψ垂直方向的结构位移dψψ对角/斜向纹理细节这里的关键是 ψ_h 不能随便挑一个高通滤波器它必须是 ψ 的 90 度相移版本只有相位差严格为 90 度b、c、d 三个通道之间才保持正交后面提取的相位才有几何意义。另一个容易忽略的点是a 通道的滤波器组合不是通常意义的 LL低通-低通而是两个相移低通的组合。这样做的理由是四个通道必须共享同一套解析包络a 的实部才能与 b、c、d 的虚部构成一致的四元数如果 a 直接取 DWT 的 LL幅值和相位之间的几何关系就不自洽后面提取的相位会带有系统性误差。2.3 用 scipy.signal.hilbert 构造 90 度相移滤波器scipy 的 hilbert 返回解析信号它的虚部就是原信号的希尔伯特变换也就是 90 度相移。直接对 4 个系数的 db2 滤波器做 hilbert 会得到严重拖尾的结果所以常见做法是补零到足够长度再做变换。import numpy as np from scipy.signal import hilbert import pywt def phase_shift_filter(f, N256): 对 FIR 滤波器做 90 度相移返回相同长度的实系数滤波器。 n len(f) f_pad np.zeros(N) start (N - n) // 2 f_pad[start:start n] f f_h np.imag(hilbert(f_pad)) out f_h[start:start n].copy() out * np.sqrt(np.sum(f * f) / np.sum(out * out)) # 保持能量量级 return out def make_analytical_wavelet(base): 给定实小波对象返回滤波器的 90 度相移版本构成的新小波。 dec_lo, dec_hi, rec_lo, rec_hi base.filter_bank return pywt.Wavelet(qwt_custom, filter_bank(phase_shift_filter(dec_lo), phase_shift_filter(dec_hi), rec_lo, rec_hi))phase_shift_filter处理的是补零后的 256 点序列取回中心 n 个系数后滤波器保持了原来的支撑位置。最后的能量归一化保证相移前后的分解系数处在同一量纲这直接影响后续多尺度分解时幅值特征的可比性。实际使用中我一般先在信号上对比 ψ 和 ψ_h 的小波系数包络确认相移滤波器的幅值响应没有明显失真再进入正式的 QWT 分解流程。3. 用 PyWavelets 从零实现 QWT分解、幅值与相位3.1 行列独立滤波器的四通道分解pywt.dwt2 默认行列使用相同滤波器而 QWT 需要行列不同所以要把二维分解拆成两步一维 DWT。先沿 axis1 分解行方向再对中间结果沿 axis0 分解列方向。def dwt2_rect(img, w_row, w_col, modeperiodization): 行、列方向使用不同小波滤波器的二维 DWT。 r_lo, r_hi pywt.dwt(img, w_row, modemode, axis1) ll pywt.dwt(r_lo, w_col, modemode, axis0)[0] lh pywt.dwt(r_lo, w_col, modemode, axis0)[1] hl pywt.dwt(r_hi, w_col, modemode, axis0)[0] hh pywt.dwt(r_hi, w_col, modemode, axis0)[1] return ll, (lh, hl, hh) def qwt2(img, waveletsym4): base pywt.Wavelet(wavelet) w_h make_analytical_wavelet(base) a dwt2_rect(img, w_h, w_h) # 实部分量 b dwt2_rect(img, base, w_h) # i 分量 c dwt2_rect(img, w_h, base) # j 分量 d dwt2_rect(img, base, base) # k 分量 return a, b, c, d先做行分解再列分解ll 是低频列与低频行的组合对应四元数实部的近似频带lh、hl、hh 分别对应水平、垂直、对角三个方向细节。默认用 sym4是因为它近似线性相位相移滤波器边界效应比 db 系列更小这个取舍在第四章展开。3.2 幅值与三个相位的计算四元数系数的极坐标表示是 q |q|(cos(θ/2) μ sin(θ/2))其中 μ 是单位纯四元数轴。把轴展开后得到三个欧拉角这里采用下面代码中的定义def qwt_band_arrays(a, b, c, d): 把四个通道的各子带组织成 shape(H, W, 4) 的四元数数组。 bands {} for idx, name in enumerate([LL, LH, HL, HH]): if idx 0: q np.stack([a[0], b[0], c[0], d[0]], axis-1) else: q np.stack([a[1][idx - 1], b[1][idx - 1], c[1][idx - 1], d[1][idx - 1]], axis-1) bands[name] q return bands def qwt_mag_phase(q): 输入 shape(H, W, 4) 的四元数系数返回幅值和三个相位。 a_, b_, c_, d_ q[..., 0], q[..., 1], q[..., 2], q[..., 3] mag np.sqrt(a_**2 b_**2 c_**2 d_**2) theta np.arctan2(np.sqrt(b_**2 c_**2 d_**2), a_) phi np.arctan2(np.sqrt(c_**2 d_**2), b_) psi np.arctan2(d_, c_) return mag, theta, phi, psitheta 的取值范围是 0 到 π描述主相移phi 是 -π/2 到 π/2描述虚部内部比例psi 是 -π 到 π描述 d 相与 c 相的相对转角。三个相位的轴序与 2.2 节表格一致水平位移对应 phi、垂直位移对应 psi。实际文献中不同作者会对换轴序跑通流程后应先固定自己的轴序约定避免后期特征对照时出现相位符号不一致。3.3 多层 QWT 的金字塔分解结构单层 QWT 返回 16 个二维数组4 个通道 × 4 个子带。下一层分解的常见做法是对四个通道的 LL 子带再执行 qwt2也就是把上一层的近似系数当成图像继续分解def qwt2_levels(img, waveletsym4, level3): coeffs [] for _ in range(level): a, b, c, d qwt2(img, wavelet) coeffs.append((a, b, c, d)) img a[0] # 以实部的 LL 作为下一层输入 return coeffs这个递归方式与 DWT 的多分辨率分析一致只是每层多了三个通道的细节。优点是每一层子带尺寸减半统计特征自然形成金字塔缺点是特征维数随层数线性增长层数超过 4 后统计特征里会混入大量边界相位噪声。输入尺寸如果不是 2 的整数次幂分解前要先修剪或填充到 2^level 的整数倍否则 periodization 模式的输出长度会不符合预期。4. QWT 参数整定滤波器、层数与边界模式4.1 小波滤波器选择长度与相移精度的权衡相移滤波器由原滤波器经 FFT 域希尔伯特变换得到滤波器越长频域采样点越多90 度相移的精度越高。但支撑域长了边界伪影的扩散范围也变大所以要在相位精度与边界之间折中。滤波器长度相移精度边界伪影建议haar/db12差轻不推荐db24一般较重快速原型db4 / sym48较好中等默认db8 / sym816好较重大图或精度优先sym4 与 db4 长度相同但 sym4 的滤波器更接近线性相位相移后支撑中心偏移更小镜像延拓所需的边沿点也更少。在同样的 128×128 输入、三层分解条件下sym4 的相位图在子带边缘的暗纹比 db4 少这个差异在视觉效果上很直观也是我默认用 sym4 的原因。4.2 分解层数按最小子带尺寸反推层数每加 1最低频子带边长减半。相位数在边长 8 以下的子带上噪声占比过高特征基本是统计噪声。按这个约束建议值如下表输入短边建议层数最低子带边长642161283162563325124324.3 边界模式必须选 periodization这一条很容易被忽略。symmetric 延拓在图像边缘制造镜像对称对普通 DWT 无所谓但对 QWT 影响明显相移滤波器在镜像边界处产生的响应与内部真实响应相位相反会把相位特征从边缘开始污染第一层全部子带。periodization 假设信号周期延拓与希尔伯特变换的周期假设一致。代价是要求输入边长能被 2^level 整除所以代码里通常先做一次预处理h, w img.shape[:2] level QWT_CONFIG[level] img img[:h - h % (2 ** level), :w - w % (2 ** level)]这行裁剪保证 periodization 不报长度错误也避免了分解完子带尺寸不齐的问题。还要注意periodization 模式下 pywt 对无法整除的输入会先复制到最近的可分解长度导致输出尺寸和预期不一致提前裁剪可以避免这种隐式行为也让不同实验的系数形状完全对齐。4.4 一个可复现的参数配置把可调参数集中到一个字典里便于实验对照QWT_CONFIG { wavelet: sym4, # 8 tap兼顾相位精度与边缘伪影 level: 3, # 128x128 输入时最低子带 16x16 mode: periodization, hilbert_pad: 256, # 相移滤波器的 FFT 补零长度 }需要换输入尺寸时只改 level 与裁剪逻辑不动分解主函数。hilbert_pad 只在构造滤波器时生效对分解速度几乎没有影响。5. 实战QWT 幅值 相位特征做纹理分类5.1 构造四类合成纹理数据集为了验证 QWT 特征的有效性先用合成数据做验证四类纹理容易分离且可控——竖直条纹、水平条纹、45 度斜纹和随机噪声。每类生成 32 张 128×128 图频率随机取 2 到 6再叠加高斯噪声def make_texture(kindv, size128): y, x np.mgrid[:size, :size] freq np.random.choice([2, 3, 4, 5, 6]) angle {v: 0, h: np.pi / 2, d: np.pi / 4}.get(kind) if kind n: img np.random.randn(size, size) else: img np.sin(2 * np.pi * freq * (x * np.cos(angle) y * np.sin(angle))) return img 0.05 * np.random.randn(size, size)合成数据最大的好处是类别标签绝对干净也方便在消融实验时判断分类器到底在用什么线索下结论。5.2 特征定义幅值统计 相位 sin/cos 编码相位是角量直接对角度求均值会遇到 -π 和 π 附近跳变的问题所以把每个相位编码成 sin 与 cos 两个实数再统计。每层 4 个子带每个子带特征为幅值 mean/std 三个相位的 sin 均值、cos 均值、sin 标准差、cos 标准差共 8 维。三层 QWT 得到一个 96 维特征向量。def qwt_feature_vec(img, config): feats [] for a, b, c, d in qwt2_levels(img, config[wavelet], config[level]): for name in [LL, LH, HL, HH]: q qwt_band_arrays(a, b, c, d)[name] mag, th, ph, ps qwt_mag_phase(q) feats.append([mag.mean(), mag.std()]) for ang in (th, ph, ps): s, c np.sin(ang), np.cos(ang) feats [s.mean(), c.mean(), s.std(), c.std()] return np.concatenate(feats)特征里没有直接用相位均值是避免角度环绕造成均值漂移sin/cos 编码相当于把角度投影到单位圆上求出来的均值方向才是真正的主相位方向。5.3 SVM 分类与消融对照数据生成、分解、特征提取后训练一个 RBF 核 SVM固定随机种子保证可复现from sklearn.svm import SVC from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split X, y [], [] for kind, label in [(v, 0), (h, 1), (d, 2), (n, 3)]: for _ in range(32): X.append(qwt_feature_vec(make_texture(kind), QWT_CONFIG)) y.append(label) X np.asarray(X) X_tr, X_te, y_tr, y_te train_test_split( X, y, test_size0.25, stratifyy, random_state0) scaler StandardScaler().fit(X_tr) clf SVC(kernelrbf, C10).fit(scaler.transform(X_tr), y_tr) print(clf.score(scaler.transform(X_te), y_te))在这个配置下幅值相位联合特征分类准确率通常在 0.95 以上单独用幅值统计时分类器容易把 45 度斜纹与竖直条纹混淆因为两者在幅值分布上接近而相位分布差别明显。单独用相位时随机噪声类容易被误判成某一种定向纹理因为噪声系数幅值小相位对小幅值噪声系数敏感而幅值统计能识别“这是无方向的能量”。6. 平移验证实验QWT 幅值的稳定性与相位-位移耦合6.1 DWT 与 QWT 对 1 像素平移的幅值稳定性拿一张 45 度条纹图对图像做 1 像素的水平循环平移分别用 db4 DWT 和 QWT 提取 HH 子带幅值比较归一化变化量img make_texture(d) img_shift np.roll(img, shift1, axis1) _, (_, _, d_hh) pywt.dwt2(img, db4, modeperiodization) _, (_, _, d_hh_s) pywt.dwt2(img_shift, db4, modeperiodization) dwt_delta np.abs(np.linalg.norm(d_hh) - np.linalg.norm(d_hh_s)) / np.linalg.norm(d_hh) b1 qwt_band_arrays(*qwt2(img))[HH] b2 qwt_band_arrays(*qwt2(img_shift))[HH] qmag, qmag_s qwt_mag_phase(b1)[0], qwt_mag_phase(b2)[0] qwt_delta np.abs(np.linalg.norm(qmag) - np.linalg.norm(qmag_s)) / np.linalg.norm(qmag) print(fDWT 幅值变化: {dwt_delta:.2%}, QWT 幅值变化: {qwt_delta:.2%})注意这里的对比必须使用同一组图像唯一区别是平移一个像素否则幅值差异主要来自纹理本身的空间不均匀性测试就失去意义。DWT 的 HH 系数对平移非常敏感幅值变化通常在 20% 以上QWT 把位移信息放进了相位幅值变化通常能压到 5% 以内。这个对照是验证 QWT 实现是否正确的第一道检查如果 QWT 幅值变化也很大大概率是相移滤波器构造或边界模式的问题。6.2 相位差与位移的线性关系继续用上一节的 QWT 分解结果计算 HH 子带相位 theta 的差异。theta 是周期角先做 mod 到 [-π, π] 再取平均th, th_s qwt_mag_phase(b1)[1], qwt_mag_phase(b2)[1] diff (th - th_s np.pi) % (2 * np.pi) - np.pi print(fHH 子带平均相位差: {diff.mean():.3f} rad)对单一频率的条纹这个相位差近似等于位移量与空间频率的乘积频率越高相位差越大。这就是 QWT 相位能用于局部亚像素配准的核心性质。把测试图换成随机噪声时相位差会变得杂乱因为噪声包含所有频率分量此时应按子带分别统计更有意义。这类平移测试应当写进 QWT 工具库的冒烟测试函数里每次改动滤波器参数后先跑一遍再继续调特征。如果拿到的参考实现比如命名类似 hingfui.zip 的源码包相位差不满足随频率变化的单调性优先排查两处相移滤波器是否做了能量归一化边界模式是否被改成了 reflect 或 symmetric 而不是 periodization。本文还有配套的精品资源点击获取