python装小波变换包 python 小波变换工具

 

信号分析是为了时间和频率之间的相互关系。傅里叶变换提供了有关频域的信息，但有关时间的局部化信息却基本丢失。傅里叶变换只能提供信号在整个时域上的频率，不能提供信号在某个局部时间段上的频率信息。与傅里叶变换不同，小波变换是通过缩放母小波（Mother Wavelet）的宽度获得信号的频域特征，通过平移母小波获得信号的时间信息。对母小波的缩放和平移是为了计算小波系数，这些小波系数反映了小波和局部信号之间的相关程度。小波分析就是把一个信号分解为将母小波经过缩放和平移之后的一系列小波，因此小波是小波变换的基函数。小波变换可以理解为用经过缩放和平移的一系列小波函数，代替傅里叶变换的正弦波和余弦波进行傅里叶变换的结果。与傅里叶变换相比，小波变换是空间（时间）和频率的局部变换，通过伸缩平移运算，对信号逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节。

• 尺度函数 ： scaling function （在一些文档中又称为父函数 father wavelet ）
• 小波函数 ： wavelet function（在一些文档中又称为母函数 mother wavelet）
• 小波变换是将原始图像与小波基函数以及尺度函数进行内积运算,所以一个尺度函数和一个小波基函数就可以确定一个小波变换

从时域上看，相差很大的两个信号，在频域上却非常相近。一个很自然的方法是加窗（短时距傅里叶变换），将长时间信号分成数个较短的等长信号，然后再分别对每个窗进行傅里叶变换，从而得到频率随时间的变化，这就是短时距傅里叶变换

 

 不仅可以看到信号中有哪些频率，还可以看到不同的频率成分在什么时间出现。傅里叶变换类似棱镜，可以将不同的信号分解。小波变换类似于显微镜，不仅知道信号中有哪些成分，还可以知道各种成分在什么位置出现。

 

 

小波分析优于傅里叶分析之处在于：小波分析在时域和频域同时具有良好的局部化性质，因为小波函数是紧支集，而正弦、余弦的区间是无穷区间，所以小波变换可以对高频成分采用逐渐精细的时域或空域取代步长，从而可以聚焦到对象的任意细节。

 通过小波变换得到N个小波系数，而二维离散小波变换输入的是二维矩阵，每个步骤输出的是近似图像，水平细节，垂直细节和对角细节

 

以haar小波变换过程为例来理解小波变换。

例：求只有4个像素[9 7 3 5]的图像的哈尔小波变换系数。 计算步骤如下：

步骤1：求均值(averaging)，也叫Approximation。计算相邻像素对的平均值，得到一幅分辨率比较低的新图像，新的图像的分辨率是原来的1/2，相应的像素值为：[8 4]

步骤2：求差值(differencing)。很明显，用2个像素表示这幅图像时，图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像素组成的原始图像，就需要存储一些图像的细节系数(detail coefficient)，以便在重构时找回丢失的信息。方法是使用这个像素对的差值除以2，结果为[8 4 1 -1]
步骤3：重复第1，2步，把由第一步分解得到的图像进一步分解成分辨率更低的图像和细节系数。[6 2 1 -1]

经过小波变换后图像会生成低频信息和高频信息。低频信息对应于求均值，高频信息对应于求差值。
均值是局部的平均值，变化缓慢，属于低频信息，存储图片的轮廓信息，近似信息差值是局部的波动值，变化较快，属于高频信息，存储图片的细节信息，局部信息，另外含有噪音

 

 二维离散小波变换

 

 

 

1、二维图像单级变换：dwt2() 

haar小波单级变换之后：
低频信息的取值范围为：[0,510]
高频信息的取值范围为：[-255,255]

import numpy as np
import pywt
import cv2
import matplotlib.pyplot as plt

img = cv2.imread("cat.jpg")
img = cv2.resize(img, (448, 448))
# 将多通道图像变为单通道图像
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY).astype(np.float32)

plt.figure('二维小波一级变换')
coeffs = pywt.dwt2(img, 'haar')
cA, (cH, cV, cD) = coeffs

# 将各个子图进行拼接，最后得到一张图
AH = np.concatenate([cA, cH], axis=1)
VD = np.concatenate([cV, cD], axis=1)
img = np.concatenate([AH, VD], axis=0)

# 显示为灰度图
plt.imshow(img,'gray')
plt.title('result')
plt.show()

在拼接子图之前，应该先对各个子图进行处理。未处理的情况下，因为高频部分的像素值极小甚至小于0，所以高频区域呈黑色。

最简单的处理方式为：将高频信息均加255，得到如下结果：

 另一种显示方式：

import numpy as np
import pywt
import cv2
import matplotlib.pyplot as plt

img = cv2.imread("cat.jpg")
img = cv2.resize(img, (448, 448))
# 将多通道图像变为单通道图像
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY).astype(np.float32)

plt.figure('二维小波一级变换')
coeffs = pywt.dwt2(img, 'haar')
cA, (cH, cV, cD) = coeffs

plt.subplot(221), plt.imshow(cA, 'gray'), plt.title("A")
plt.subplot(222), plt.imshow(cH, 'gray'), plt.title("H")
plt.subplot(223), plt.imshow(cV, 'gray'), plt.title("V")
plt.subplot(224), plt.imshow(cD, 'gray'), plt.title("D")
plt.show()

 2、二维图像多级分解：wavedec2()

import numpy as np
import pywt
import cv2
import matplotlib.pyplot as plt

img = cv2.imread("cat.jpg")
img = cv2.resize(img, (448, 448))
# img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB).astype(np.float32)
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY).astype(np.float32)

plt.figure('二维图像多级分解')
coeffs = pywt.wavedec2(img, 'haar', level=2)
cA2, (cH2, cV2, cD2), (cH1, cV1, cD1) = coeffs

# 将每个子图的像素范围都归一化到与CA2一致  CA2 [0,255* 2**level]
AH2 = np.concatenate([cA2, cH2+510], axis=1)
VD2 = np.concatenate([cV2+510, cD2+510], axis=1)
cA1 = np.concatenate([AH2, VD2], axis=0)

AH = np.concatenate([cA1, (cH1+255)*2], axis=1)
VD = np.concatenate([(cV1+255)*2, (cD1+255)*2], axis=1)
img = np.concatenate([AH, VD], axis=0)

plt.imshow(img,'gray')
plt.title('2D WT')
plt.show()

更多例子：

 

import numpy as np
import pywt
import cv2
import matplotlib.pyplot as plt


def haar_img():
    img_u8 = cv2.imread("len_std.jpg")
    img_f32 = cv2.cvtColor(img_u8, cv2.COLOR_BGR2GRAY).astype(np.float32)

    plt.figure('二维小波一级变换')
    coeffs = pywt.dwt2(img_f32, 'haar')
    cA, (cH, cV, cD) = coeffs

    # 将各个子图进行拼接，最后得到一张图
    AH = np.concatenate([cA, cH], axis=1)
    VD = np.concatenate([cV, cD], axis=1)
    img = np.concatenate([AH, VD], axis=0)
    return img

if __name__ == '__main__':
    img = haar_img()

    plt.imshow(img, 'gray')
    plt.title('img')
    plt.show()

 Daubechies （db2）演示

import numpy as np
import pywt
from matplotlib import pyplot as plt
from pywt._doc_utils import wavedec2_keys, draw_2d_wp_basis

x = pywt.data.camera().astype(np.float32)
shape = x.shape
max_lev = 3       # 要绘制多少级分解
label_levels = 3  # 图上要显式标注多少层

fig, axes = plt.subplots(2, 4, figsize=[14, 8])
for level in range(0, max_lev + 1):
    if level == 0:
        # 显示原始图像
        axes[0, 0].set_axis_off()
        axes[1, 0].imshow(x, cmap=plt.cm.gray)
        axes[1, 0].set_title('Image')
        axes[1, 0].set_axis_off()
        continue

    # 绘制标准DWT基的子带边界
    # draw_2d_wp_basis(shape, wavedec2_keys(level), ax=axes[0, level], label_levels=label_levels)
    # axes[0, level].set_title('{} level\ndecomposition'.format(level))

    # 计算二维 DWT
    c = pywt.wavedec2(x, 'db2', mode='periodization', level=level)

    # 独立规范化每个系数数组以获得更好的可见性
    c[0] /= np.abs(c[0]).max()
    for detail_level in range(level):
        c[detail_level + 1] = [d / np.abs(d).max() for d in c[detail_level + 1]]

    # 显示归一化系数 (normalized coefficients)
    arr, slices = pywt.coeffs_to_array(c)
    axes[1, level].imshow(arr, cmap=plt.cm.gray)
    axes[1, level].set_title('Coefficients\n({} level)'.format(level))
    axes[1, level].set_axis_off()
plt.tight_layout()
plt.show()

 Biorthogonal （bior）演示

 

import numpy as np
import matplotlib.pyplot as plt
import pywt
import pywt.data


def test_pywt():
    original = pywt.data.camera()

    # 图像的小波变换，图形的逼近和细节
    titles = ['Approximation', ' Horizontal detail', 'Vertical detail', 'Diagonal detail']
    coeffs2 = pywt.dwt2(original, 'bior1.3')
    LL, (LH, HL, HH) = coeffs2
    plt.imshow(original)
    plt.colorbar(shrink=0.8)
    fig = plt.figure(figsize=(12, 3))
    for i, a in enumerate([LL, LH, HL, HH]):
        ax = fig.add_subplot(1, 4, i + 1)
        ax.imshow(a, interpolation="nearest", cmap=plt.cm.gray)
        ax.set_title(titles[i], fontsize=10)
        ax.set_xticks([])
        ax.set_yticks([])

    fig.tight_layout()
    plt.show()


if __name__ == '__main__':
    test_pywt()