简介
在本实验中，我们将学习如何使用稀疏编码方法将信号转换为雷克子波（Ricker wavelet）的稀疏组合。雷克子波（也称为墨西哥帽小波或高斯函数的二阶导数）并不是表示此类分段常数信号的特别好的核函数。因此，可以看出添加不同宽度的原子有多重要，这也促使我们学习最适合信号类型的字典。
我们将使用 SparseCoder 估计器直观地比较不同的稀疏编码方法。右边更丰富的字典在大小上并没有更大，为了保持在相同的数量级，我们进行了更重的下采样。
虚拟机使用提示
虚拟机启动完成后，点击左上角切换到 笔记本 标签页，以访问 Jupyter Notebook 进行练习。
有时，你可能需要等待几秒钟让 Jupyter Notebook 完成加载。由于 Jupyter Notebook 的限制，操作的验证无法自动化。
如果你在学习过程中遇到问题，可以随时向 Labby 提问。课程结束后提供反馈，我们会及时为你解决问题。
导入库
我们将首先导入必要的库。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import SparseCoder
定义雷克子波函数
我们将定义生成雷克子波及雷克子波字典的函数。
def ricker_function(resolution, center, width):
"""离散下采样的雷克子波（墨西哥帽小波）"""
x = np.linspace(0, resolution - 1, resolution)
x = (
(2 / (np.sqrt(3 * width) * np.pi**0.25))
* (1 - (x - center) ** 2 / width**2)
* np.exp(-((x - center) ** 2) / (2 * width**2))
)
return x
def ricker_matrix(width, resolution, n_components):
"""雷克子波（墨西哥帽小波）字典"""
centers = np.linspace(0, resolution - 1, n_components)
D = np.empty((n_components, resolution))
for i, center in enumerate(centers):
D[i] = ricker_function(resolution, center, width)
D /= np.sqrt(np.sum(D**2, axis=1))[:, np.newaxis]
return D
生成信号
我们将生成一个信号，并使用 Matplotlib 对其进行可视化。
resolution = 1024
subsampling = 3 ## 下采样因子
width = 100
n_components = resolution // subsampling
## 生成一个信号
y = np.linspace(0, resolution - 1, resolution)
first_quarter = y < resolution / 4
y[first_quarter] = 3.0
y[np.logical_not(first_quarter)] = -1.0
## 可视化信号
plt.figure()
plt.plot(y)
plt.title("原始信号")
plt.show()
计算小波字典
我们将计算一个小波字典，并使用 Matplotlib 对其进行可视化。
## 计算一个小波字典
D_fixed = ricker_matrix(width=width, resolution=resolution, n_components=n_components)
D_multi = np.r_[
tuple(
ricker_matrix(width=w, resolution=resolution, n_components=n_components // 5)
for w in (10, 50, 100, 500, 1000)
)
]
## 可视化小波字典
plt.figure(figsize=(10, 5))
for i, D in enumerate((D_fixed, D_multi)):
plt.subplot(1, 2, i + 1)
plt.imshow(D, cmap=plt.cm.gray, interpolation="nearest")
plt.title("小波字典 (%s)" % ("固定宽度" if i == 0 else "多种宽度"))
plt.axis("off")
plt.show()
稀疏编码
我们将使用不同方法对信号执行稀疏编码，并可视化结果。
## 按以下格式列出不同的稀疏编码方法：
## (标题，变换算法，变换α,
## 变换非零系数数量，颜色)
estimators = [
("正交匹配追踪（OMP）", "omp", None, 15, "海军蓝"),
("套索（Lasso）", "lasso_lars", 2, None, "青绿色"),
]
lw = 2
plt.figure(figsize=(13, 6))
for subplot, (D, title) in enumerate(
zip((D_fixed, D_multi), ("固定宽度", "多种宽度"))
):
plt.subplot(1, 2, subplot + 1)
plt.title("针对 %s 字典的稀疏编码" % 标题)
plt.plot(y, lw=lw, linestyle="--", label="原始信号")
## 进行小波逼近
for title, algo, alpha, n_nonzero, color in estimators:
coder = SparseCoder(
dictionary=D,
transform_n_nonzero_coefs=n_nonzero,
transform_alpha=alpha,
transform_algorithm=algo,
)
x = coder.transform(y.reshape(1, -1))
density = len(np.flatnonzero(x))
x = np.ravel(np.dot(x, D))
squared_error = np.sum((y - x) ** 2)
plt.plot(
x,
color=color,
lw=lw,
label="%s: %s 个非零系数，\n%.2f 误差" % (标题, density, squared_error),
)
## 软阈值去偏
coder = SparseCoder(
dictionary=D, transform_algorithm="threshold", transform_alpha=20
)
x = coder.transform(y.reshape(1, -1))
_, idx = np.where(x!= 0)
x[0, idx], _, _, _ = np.linalg.lstsq(D[idx, :].T, y, rcond=None)
x = np.ravel(np.dot(x, D))
squared_error = np.sum((y - x) ** 2)
plt.plot(
x,
color="深橙色",
lw=lw,
label="带去偏的阈值处理:\n%d 个非零系数，%.2f 误差"
% (len(idx), squared_error),
)
plt.axis("tight")
plt.legend(shadow=False, loc="最佳")
plt.subplots_adjust(0.04, 0.07, 0.97, 0.90, 0.09, 0.2)
plt.show()
总结
在本实验中，我们学习了如何使用稀疏编码方法和 SparseCoder 估计器将信号变换为雷克子波的稀疏组合。我们还比较了不同的稀疏编码方法并可视化了结果。