2026a

# modwtmra

基于 MODWT 的多分辨率分析

函数库: TyWavelet

# 语法

mra = modwtmra(w)
mra = modwtmra(w, wname)
mra = modwtmra(w, Lo, Hi)
mra = modwtmra(___; boundary = "reflection")

# 说明

mra = modwtmra(w) 返回最大重叠离散小波变换（MODWT）矩阵 w 的多分辨率分析（MRA）。 MODWT 矩阵 w 是 modwt 函数的输出。默认情况下，modwtmra 假设你使用带有周期性边界处理的 "sym4" 小波获得 w。示例

mra = modwtmra(w, wname) 使用与 wname 对应的小波构建 MRA。wname 小波必须是用于获得 MODWT 的同一小波。示例

mra = modwtmra(w, Lo, Hi) 使用尺度滤波器 Lo 和小波滤波器 Hi 构建 MRA。Lo 和 Hi 滤波器必须是用于获得 MODWT 的相同滤波器。示例

mra = modwtmra(___; boundary = "reflection") 在构建 MRA 的过程中使用反射边界条件，使用前面语法中的任何参数。如果你指定 "reflection"，modwtmra 假定 w 的列维是偶数，等于原始信号长度的两倍。示例

如果你用小波管理器添加了一个名为 "reflection" 的小波，你必须在使用这个选项之前重新命名该小波。默认情况下，modwtmra 在边界使用周期性扩展。

# 示例

使用 MODWTMRA 进行完美重建

获取简单时间序列信号的 MODWTMRA 并演示完美重构。

创建时间序列信号

using TyWavelet
t = 1:10
x = sin.(2 * pi * 2 * t)

获取 MODWT 和 MODWTMRA，并对 MODWTMRA 行求和。

m = modwt(x)
mra = modwtmra(m)
xrec = sum(mra; dims=1)

用绝对值的最大值来说明原始信号和重建信号之间的差异极小。最大绝对值约为，这表明重建是完美的。

maximum(abs.(x .- xrec))

4.4087284769323405e-15

使用非默认小波的 MRA

使用 db2 小波构建心电信号四级以下的 MRA。数据取自 Percival & Walden (2000)，第 125 页（数据最初由华盛顿大学的 William Constantine 和 Per Reinhall 提供）。心电信号的采样频率为 180 赫兹。

using TyWavelet
using TyBase
using TyPlot
pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/wecg.mat"
y = load(source_path)
wecg = y["wecg"]
lev = 4
wtecg = modwt(wecg, "db2", lev)
mra = modwtmra(wtecg, "db2")

绘制心电图波形和 MRA 图。

t = (0:(length(wecg) - 1)) / 180
subplot(6, 1, 1)
plot(t, wecg)
for kk in 2:(lev + 2)
    subplot(6, 1, kk)
    plot(t, mra[kk - 1, :])
end
xlabel("Time (s)")

使用默认小波的 MRA

对南方涛动指数数据进行多分辨率分析。取样周期为一天。绘制与天的比例相对应的8级细节。这一尺度的细节捕捉到了大约一年的涛动。

using TyWavelet
using TyBase
using TyPlot
pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/soi.mat"
y = load(source_path)
soi = y["soi"]

wtsoi = modwt(soi)
mrasoi = modwtmra(wtsoi)
plot(mrasoi[8, :])
title("Level 8 Details")

使用最小带宽缩放和小波滤波器的 MRA

使用最小带宽缩放和具有四个系数的小波滤波器，获取Deutsch Mark - U.S. Dollar 汇率数据的 MRA。

using TyWavelet
using TyBase
pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/DM_USD.mat"
y = load(source_path)
DM_USD = y["DM_USD"]
Lo = [0.4801755; 0.8372545; 0.2269312; -0.1301477]
Hi = qmf(Lo)
wdm = modwt(DM_USD, Lo, Hi)
mra = modwtmra(wdm, Lo, Hi)

使用反射边界的 MRA

使用 "反射 "边界处理获取心电信号的 MRA。数据摘自 Percival & Walden (2000)，第 125 页（数据最初由华盛顿大学的 William Constantine 和 Per Reinhall 提供）。

using TyWavelet
using TyBase
pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/wecg.mat"
y = load(source_path)
wecg = y["wecg"]

wtecg = modwt(wecg; boundary="reflection")
mra = modwtmra(wtecg; boundary="reflection")

证明 MRA 的列数等于原始信号的元素数。

isequal(size(mra, 2), length(wecg))

true

多信号 MRA

加载 Espiga3 [3] 的 23 个脑电图通道数据。通道按列排列。数据采样频率为 200 Hz。

using TyWavelet
using TyBase
pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/Espiga3.mat"
y = load(source_path)
Espiga3 = y["Espiga3"]

获取多信号的 MRA。

w = modwt(Espiga3)
mra = modwtmra(w)

MODWT 与 MODWTMRA 的比较

本例演示了函数 MODWT 和 MODWTMRA 之间的区别。MODWT 将信号能量划分为细节系数和缩放系数。MODWTMRA 将信号投影到小波子空间和缩放子空间。

选择 sym6 小波。加载并绘制心电图（ECG）信号。心电图信号的采样频率为 180 赫兹。数据摘自 Percival 和 Walden (2000)，第 125 页（数据最初由华盛顿大学的 William Constantine 和 Per Reinhall 提供）。

using TyWavelet
using TyBase
using TyPlot
pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/wecg.mat"
y = load(source_path)
wecg = y["wecg"]
t = (0:(length(wecg) - 1)) / 180
wv = "sym6"
plot(t, wecg)
grid("on")
title("Signal Length = $(length(wecg))")
xlabel("Time (s)")
ylabel("Amplitude")

求信号的 MODWT。

wtecg = modwt(wecg, wv)

输入数据是函数在多个时间点上的评估样本。该函数可表示为不同尺度和平移的缩放函数\phi (x)和小波的线性组合：，其中，和是小波分解的级数。第一个求和是信号的粗比例近似值，而是连续尺度的细节。MODWT 返回个系数和个细节系数的扩展。wtecg 中的每一行都包含不同比例的系数。

在对长度为的信号进行 MODWT 处理时，有多级分解（默认）。每一级都会生成详细系数。仅返回最后一级的缩放系数。在本例中，由于，，wtecg 中的行数为。

MODWT 将能量划分为不同的尺度和缩放系数：，其中，为输入数据，是尺度上的细节系数，而是最终缩放系数。

计算每个尺度的能量，并评估它们的总和。

energy_by_scales = vec(sum(wtecg .^ 2; dims=2))
Levels = [
    "D1"
    "D2"
    "D3"
    "D4"
    "D5"
    "D6"
    "D7"
    "D8"
    "D9"
    "D10"
    "D11"
    "A11"
]
energy_table = DataFrame(; Levels, energy_by_scales)

12×2 DataFrame
 Row │ Levels  energy_by_scales 
     │ String  Float64
─────┼──────────────────────────
   1 │ D1              14.0635
   2 │ D2              20.6121
   3 │ D3              37.7156
   4 │ D4              25.1233
   5 │ D5              17.4368
   6 │ D6               8.98517
   7 │ D7               1.29062
   8 │ D8               4.7278
   9 │ D9              12.2052
  10 │ D10             76.4283
  11 │ D11             76.2684
  12 │ A11              3.41918

energy_total = sum(energy_by_scales)

298.2758926235762

通过计算信号的能量，并将其与所有尺度的能量之和进行比较，确认 MODWT 是保能量的。

energy_ecg = sum(wecg .^ 2)
maximum(abs.(energy_total .- energy_ecg))

7.442508831445593e-10

提取信号的 MODWTMRA。

mraecg = modwtmra(wtecg, wv)

MODWTMRA 返回函数在不同小波子空间和最终缩放空间上的投影。也就是说，MODWTMRA 返回和在个时间点求值的个 $\left { f_{j}(x) \right } $。mraecg 中的每一行都是在不同子空间上的投影。这意味着将所有投影相加即可恢复原始信号。而 MODWT 的情况并非如此。添加 wtecg 中的系数并不能恢复原始信号。

选择一个时间点，添加在该时间点评估的的投影，并与原始信号进行比较。

time_point = 1000
abs.(sum(mraecg[:, time_point]) .- wecg[time_point])

3.0876690093606385e-13

确认，与 MODWT 不同，MODWTMRA 不是一种能量守恒变换。

energy_ecg = sum(wecg .^ 2)
energy_mra = sum(mraecg .^ 2)
abs(energy_mra - energy_ecg)

115.7053080368704

MODWTMRA 对信号进行零相位滤波。特征将进行时间对齐。通过绘制原始信号及其投影图来证明这一点。为了更好地说明对齐情况，请放大。

plot(t, wecg, "b")
hold("on")
plot(t, mraecg[4, :], "-")
hold("off")
grid("on")
xlim([4 8])
legend("Signal", "Projection", "Location", "northwest")
xlabel("Time (s)")
ylabel("Amplitude")

使用相同比例的 MODWT 系数绘制类似的曲线图。请注意，特征将不会进行时间对齐。MODWT 并非对输入进行零相滤波。

plot(t, wecg, "b")
hold("on")
plot(t, wtecg[4, :], "-")
hold("off")
grid("on")
xlim([4 8])
legend("Signal", "Coefficients", "Location", "northwest")
xlabel("Time (s)")
ylabel("Amplitude")

# 输入参数

w - MODWT 变换

矩阵

将信号或多信号的 MODWT 变换降至 LEV 级别，分别指定为矩阵或三维数组。

对于 N 点信号的 MODWT，w 是一个 LEV+1 × N 的矩阵，对于 N × NC 的多信号的 MODWT，是一个 LEV+1 × N × NC 的数组。默认情况下，imodwt 假定你使用 "sym4" 小波获得 MODWT，并进行周期性边界处理。

wname - 合成小波

"sym4" (默认) | 字符串标量

合成小波，指定为一个字符串标量。合成小波必须是用于用 modwt 函数获得 MODWT 的小波。

Lo, Hi - 滤波器

偶数长度的实值向量

滤波器，指定为一对偶数长度的实值向量。Lo 是尺度滤波器，Hi 是小波滤波器。这些滤波器必须满足正交小波的条件。Lo 和 Hi 的长度必须相等。参见 wfilters 获取更多信息。你不能同时指定 wname 和一个滤波器对 Lo, Hi。

TIP

默认情况下，wfilters 函数会返回两对与你指定的正交或双正交小波相关的滤波器。为了与数值包中 MODWT 实现的惯例一致，当你指定一个正交小波 wname 时，modwt 函数内部使用 wfilters 返回的第二对滤波器。比如说

wts = modwt(x, "db2")

相当于

_, _, Lo, Hi = wfilters("db2"); wts = modwt(x, Lo, Hi)

这个惯例与大多数小波工具箱离散小波变换函数在分解信号时遵循的惯例不同。大多数函数在内部使用第一对滤波器。

# 输出参数

mra - 多分辨率分析

矩阵 | 三维数组

mra 是一个 LEV+1 × N 的矩阵或 LEV+1 × N × NC 的数组，其中 LEV 是 MODWT 的层级，N 是被分析信号的长度。mra 的第 k 行包含第 k 级的细节。mra 的第 (LEV+1) 行包含第 LEV 级的平滑。

默认情况下，mra 的大小与输入 w 相同。如果你指定了反射边界处理，那么 mra 的列维大小是输入 w 的一半。

# 参考文献

[1] Percival, Donald B., and Andrew T. Walden. Wavelet Methods for Time Series Analysis. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge ; New York: Cambridge University Press, 2000.

[2] Whitcher, Brandon, Peter Guttorp, and Donald B. Percival. “Wavelet Analysis of Covariance with Application to Atmospheric Time Series.” Journal of Geophysical Research: Atmospheres 105, no. D11 (June 16, 2000): 14941–62. https://doi.org/10.1029/2000JD900110.

[3] Mesa, Hector. “Adapted Wavelets for Pattern Detection.” In Progress in Pattern Recognition, Image Analysis and Applications, edited by Alberto Sanfeliu and Manuel Lazo Cortés, 3773:933–44. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005. https://doi.org/10.1007/11578079_96.

# 另请参阅

modwt | imodwt