企业官网建设流程全解析

简介：四元数小波是一种将四元数代数与小波分析相结合的高级数学工具，既保留小波在时频域的局部化分析能力，又能利用四元数有效表示信号的旋转与对称特性，在图像处理、信号去噪和通信领域中具有重要应用。资源包以hingfui.m为唯一源码文件，完整实现了四元数信号构建、常用小波基选择、四元数小波分解、多尺度系数分析以及信号逆变换重构等核心流程，可帮助读者快速建立从理论到代码的映射。整个压缩包体积约9KB，仅含1个m文件，代码精简、结构清晰，适合熟悉MATLAB基础、想要了解或验证四元数小波算法的研究人员、工程师及高年级学生参考。目前已有235人学习下载，可作为轻量级参考脚本，用于信道编码误码检测、调制方式识别、信道频率选择性衰落补偿等具体问题的算法验证与二次开发，对深入理解四元数小波的实际应用具有不错的参考价值。

1. 四元数小波：为什么二维信号需要三个相位

普通 DWT 系数是实数，回答的是“这个位置有没有结构、结构有多强”；但如果问题变成“结构相对上一帧往哪个方向移动了多远”，实数系数给不出直接答案。傅里叶变换有相位，相位能测整体位移，却把空间位置摊平。四元数小波（Quaternion Wavelet Transform，QWT）把每个小波系数从实数升级成一个四元数——一个幅值加三个相位——既保留空间定位，又把局部结构的水平、垂直、对角位移信息编码进相位里。做纹理分类、图像配准、光流和质量评价的人，经常会遇到“DWT 不够用、FT 又太全局”的尴尬，QWT 是解决这类问题最常用的工具箱。这篇文章不引入商业化库，用 PyWavelets 和 SciPy 从零搭一套可运行的 QWT 分解，再落到纹理分类和平移验证，整个路径可以直接复现。

2. 四元数小波的结构原理：从希尔伯特变换到四通道滤波

2.1 为什么 DWT 缺相位、傅里叶变换缺局部定位

DWT 的实数系数来自与实小波基的内积，只有幅值信息。对一个正弦条纹，DWT 系数会告诉你这个位置有能量，但无法区分它是往左移了 0.3 像素还是 0.8 像素。傅里叶变换的相位可以测全局平移，但它的基函数无限延伸，任何局部形变都会污染全部相位。QWT 是个折中：基函数依然有紧支撑，但系数是四元数，通过三个相位把“水平位移”、“垂直位移”和“对角协调性”分离出来。

这里还需要回答一个前置问题：为什么不用一维复小波（CWT）直接扩展？复小波把实小波变成 ψ + iψ_h，系数是复数，自带一个相位，适合一维信号的瞬时频率和包络分析。图像是二维的，一个相位无法同时区分水平与垂直两个方向的结构走向。把行、列两个方向的解析滤波结果放到四元数的三个虚部里，才能用一个数同时表达水平位移、垂直位移和对角协调性三个自由度。

2.2 可分离 QWT：四个分量各司其职

QWT 的常见构造方式是“可分离”的：把一维实小波 ψ 和它的希尔伯特变换 ψ_h 分别放到行、列方向组合，得到四个滤波通道。对图像 I 做分解后，每个位置得到一个四元数系数 q = a + bi + cj + dk，其中 a 是实部，b、c、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, N=256): """对 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。先沿 axis=1 分解行方向，再对中间结果沿 axis=0 分解列方向。

def dwt2_rect(img, w_row, w_col, mode='periodization'): """行、列方向使用不同小波滤波器的二维 DWT。""" r_lo, r_hi = pywt.dwt(img, w_row, mode=mode, axis=1) ll = pywt.dwt(r_lo, w_col, mode=mode, axis=0)[0] lh = pywt.dwt(r_lo, w_col, mode=mode, axis=0)[1] hl = pywt.dwt(r_hi, w_col, mode=mode, axis=0)[0] hh = pywt.dwt(r_hi, w_col, mode=mode, axis=0)[1] return ll, (lh, hl, hh) def qwt2(img, wavelet='sym4'): 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, psi

theta 的取值范围是 0 到 π，描述主相移；phi 是 -π/2 到 π/2，描述虚部内部比例；psi 是 -π 到 π，描述 d 相与 c 相的相对转角。三个相位的轴序与 2.2 节表格一致：水平位移对应 phi、垂直位移对应 psi。实际文献中不同作者会对换轴序，跑通流程后应先固定自己的轴序约定，避免后期特征对照时出现相位符号不一致。

3.3 多层 QWT 的金字塔分解结构

单层 QWT 返回 16 个二维数组：4 个通道 × 4 个子带。下一层分解的常见做法是对四个通道的 LL 子带再执行 qwt2，也就是把上一层的近似系数当成图像继续分解：

def qwt2_levels(img, wavelet='sym4', level=3): 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 度相移的精度越高。但支撑域长了，边界伪影的扩散范围也变大，所以要在相位精度与边界之间折中。

sym4 与 db4 长度相同，但 sym4 的滤波器更接近线性相位，相移后支撑中心偏移更小，镜像延拓所需的边沿点也更少。在同样的 128×128 输入、三层分解条件下，sym4 的相位图在子带边缘的暗纹比 db4 少；这个差异在视觉效果上很直观，也是我默认用 sym4 的原因。

4.2 分解层数按最小子带尺寸反推

层数每加 1，最低频子带边长减半。相位数在边长 8 以下的子带上噪声占比过高，特征基本是统计噪声。按这个约束，建议值如下表：

4.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(kind='v', size=128): 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_size=0.25, stratify=y, random_state=0) scaler = StandardScaler().fit(X_tr) clf = SVC(kernel='rbf', C=10).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, shift=1, axis=1) _, (_, _, d_hh) = pywt.dwt2(img, 'db4', mode='periodization') _, (_, _, d_hh_s) = pywt.dwt2(img_shift, 'db4', mode='periodization') 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(f'DWT 幅值变化: {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(f'HH 子带平均相位差: {diff.mean():.3f} rad')

对单一频率的条纹，这个相位差近似等于位移量与空间频率的乘积，频率越高相位差越大。这就是 QWT 相位能用于局部亚像素配准的核心性质。把测试图换成随机噪声时相位差会变得杂乱，因为噪声包含所有频率分量，此时应按子带分别统计更有意义。这类平移测试应当写进 QWT 工具库的冒烟测试函数里，每次改动滤波器参数后先跑一遍再继续调特征。如果拿到的参考实现（比如命名类似 hingfui.zip 的源码包）相位差不满足随频率变化的单调性，优先排查两处：相移滤波器是否做了能量归一化，边界模式是否被改成了 reflect 或 symmetric 而不是 periodization。

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