使用Python 进行信号小波分析

李彩龙

（中国人民大学信息学院 北京市 100872）

小波最先被地球物理学家用来分析通过爆炸方式产生的人造地震数据，得到地表下的岩层“图像”，以便于发现石油和进行探矿等。地震法探矿产生了许多的二维图像数据，层叠后形成一个地表下岩层的三维图像。传感器排列成线阵，爆炸发生在线阵的另一端，产生地壳中的地震波。传感器实时的收集数据，记录下地震数据。

首先被传感器记录到的是沿着地表传播的直达波。这个一般不太重要。接下来的是地表下岩层的反射波，这个比较重要。反射波到达传感器的时间信息表达了反射岩的位置，其频率反应了该岩层的结构信息。用地震法探矿的关键是如何选用适合的信号分析方法。傅立叶变换（Fourier Transform，FT）仅能提供信号频率信息，并不能给出每个时刻的频率。另外一个方法短时傅立叶变换（STFT）会好一些。这个时间域被分割成小的等时间间隔，然后分别进行傅立叶变换。分析结果包含了时间和频率信息。由于时间间隔固定，某些持续时间非常短、频率很高的脉冲信号发生的时刻很难被检测到。

小波是指“一个持续时间很短的震荡波形”，如果用ψ(t)表示一个母小波，ψ(t)的伸缩和平移也是小波，由一个母小波及其伸缩平移形成的函数集来表示一个信号，这就是小波变换。

小波变换（Wavelet Transform，WT）可以跟踪时间和频率信息，即可以检测出短时脉冲，也可以检测出长的变慢波，即具有可变的时间和频率分辨率，在低频时具有良好的频率分辨率，在高频时有良好的时间分辨率。

随着小波变换理论研究的深入和成熟，已经被广泛的应用于信号处理、图像处理、数据压缩、奇异值估计、非线性分析和雷达信号处理等方面。

Python，是一种广泛使用的解释型、高级和通用的编程语言，众多的开源科学计算软件包都提供了Python 接口，如计算机视觉库OpenCV、可视化工具库VTK 等。Python专用计算扩展库，如NumPy、SciPy、matplotlab、Pandas、scikit-learn 等。

开发工具上可用Spyder，在安装Anaconda 时已经安装好，可以从开始菜单或者使用命令行spyder 启动。Spyder是一个免费、开源的集成开发环境，用Python 编写，用于Python 程序的开发，由科学家、工程师和数据分析师设计。它综合开发工具的高级编辑、分析、调试和分析功能与科学软件包的数据探索、交互执行、深入检查和漂亮的可视化功能独特地结合在一起。可以单步执行，查看运行后的各变量的值，查看序列数据线图等，非常方便。其他的集成开发工具（IDE）还包括：

Visual Studio Code（简称 VS Code）是一款由微软开发且跨平台的免费源代码编辑器。该软件支持语法高亮、代码自动补全（又称 IntelliSense）、代码重构功能，并且内置了命令行工具和 Git 版本控制系统。用户可以更改主题和键盘快捷方式实现个性化设置，也可以通过内置的扩展程序商店安装扩展以拓展软件功能。Visual Studio Code 默认支持非常多的编程语言，要使用Python 进行开发，需要先安装相应的Python 扩展（Python extension for Visual Studio Code）。

PyCharm，主要用于Python 语言开发，由捷克公司JetBrains 开发，提供代码分析、图形化调试器，集成测试器、集成版本控制系统，并支持使用Django 进行网页开发。PyCharm 是一个跨平台开发环境，拥有Microsoft Windows、macOS 和Linux 版本。社区版在Apache 许可证下发布，另外还有专业版在专用许可证下发布，其拥有许多额外功能。

1 PyWavelets介绍

PyWavelets 是Python 的、遵循MIT 许可协议的免费开源软件，可以通过网址 https://pywavelets.readthedocs.io/ 查看相关信息。

PyWavelets 的主要特性包括：

一维、二维和 nD 前向和逆向离散小波变换（DWT 和IDWT）

一维、二维和 nD 多级 DWT 和 IDWT

一维、二维和 nD 平稳小波变换（Undecimated Wavelet Transform）

一维和二维小波包分解与重构

一维连续小波变换

计算小波和标度函数的近似值

超过100 个内置小波滤波器并支持自定义小波

单精度和双精度计算

真实复杂的计算

结果与 Matlab Wavelet Toolbox (TM) 兼容

2 PyWavelets安装方式

在命令行中使用pip install PyWavelets

使用Anaconda Python 分发版的用户可以使用如下命令安装

conda install pywavelets

3 PyWavelets的使用

PyWavelets 模块导入语句为：import pywt

4 小波族及小波名称

小波分析是会涉及到不同的小波族，在PyWavelets 中，支持的小波族名称可以通过如下的方法调用获得：

pywt.families(short=True)

返回可用内置小波族的列表。

参数short：布尔类型，可选，表示是否使用简写，默认为True

返回：列表类型，支持的小波组名字列表

目前内置的小波族有：

获得小波族下具体的小波名称的方法如下，在这个方法中可以获得具体可用的小波名称，后边在进行小波分析的时候都需要传递小波名称的参数或者小波对象。

pywt.wavelist(family=None, kind='all')

返回给定族名称的可用子波名称列表。

参数 family : 字符串类型, 可选。表示简写的小波族名称。如果族名称为None（默认），则返回所有内置小波的名称。否则，该函数返回属于给定族的小波的名称。有效名称为：'haar', 'db', 'sym', 'coif', 'bior', 'rbio', 'dmey', 'gaus', 'mexh','morl', 'cgau', 'shan', 'fbsp', 'cmor'。

参数 kind : 字符串类型，可选，表示小波族的类型，可选项为{‘all’, ‘continuous’, ‘discrete’}，默认为all。如果指定了小波族名字，则忽略该参数。否则指的是连续小波（'continuous'）或者离散小波（'discrete'）。

返回字符串列表类型，可用小波名称列表。

例如：pywt.wavelist('coif') 返回：['coif1', 'coif2', 'coif3','coif4', 'coif5', 'coif6', 'coif7', 'coif8', 'coif9', 'coif10', 'coif11','coif12', 'coif13', 'coif14', 'coif15', 'coif16', 'coif17']。

pywt.wavelist(kind='continuous') 返回：['cgau1', 'cgau2','cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7', 'cgau8', 'cmor', 'fbsp','gaus1', 'gaus2', 'gaus3', 'gaus4', 'gaus5', 'gaus6', 'gaus7', 'gaus8','mexh', 'morl', 'shan']。

5 一维连续小波变换（Continuous Wavelet Transform）

连续小波转换可以建构一个具有良好时域和频域局部化的时频讯号，能看到每个时刻的频率分布情况。在实际使用中，结合matplotlib.pyplot.pcolormesh（用于绘制频谱图）和matplotlib.pyplot.contour（用于绘制等高线）可以呈现出丰富的表达。在工作中发现，pcolormesh 函数的cmap 参数也至关重要，不同的值有不同的表现力，需要多实验，寻找到适合的取值。

cwt(data, scales, wavelet, sampling_period=1.,method='conv', axis=-1)

参数：data：输入的数据，类数组类型

scales：类数组，所使用的小波尺度。可以使用f =scale2frequency(wavelet, scale)/sampling_period 来确定什么物理频率f。这里，当sampling_period 以秒为单位时，f 的单位是赫兹。

wavelet：所用的小波对象或者名称

sampling_period：浮点数，频率输出的采样周期（可选）。coefs 的计算的值与sampling_period 的选择无关（即，scales不按采样周期缩放）。

method：取值为{‘conv’, ‘fft’}, 可选，用于计算CWT 的方法。可以是以下任何一种：

conv：使用numpy.convolve。

fft：使用频域卷积。

auto：使用基于对每个尺度的计算复杂性的估计进行自动选择。

conv 方法的复杂度为O(len(scale) * len(data))。fft 方法是O(N * log2(N))，N = len(scale) + len(data) - 1。它非常适合大尺寸信号，但在小尺寸信号上略慢于conv。

axis：整型，可选。用于计算CWT 的轴。如果没有给出，则使用最后一个轴。

返回值：

coefs：类数组类型，对于给定的尺度和小波的输入信号的连续小波变换。coefs 的第一轴对应于尺度。剩余的轴与data 的形状相匹配。

frequencies：类数组类型，如果采样周期单位为秒且给定，则频率单位为赫兹。否则，假设采样周期为1。

在实际的使用中，可以使用如下的方式计算：

长期以来，离散小波变换在数字信号处理、地震预报、石油勘探、医学断层处理、编码理论、量子物理及概率论等领域中都有广泛的应用。使用离散小波变换，先得到每个频段的系数，然后由系数重建波形，这样就可以看到每个频段的数据波形，方便观察和进行特征提取。

6 单级离散小波变换

pywt.dwt(data, wavelet, mode='symmetric', axis=-1)

其中：data 是输入待分析信号，类数组类型；

wavelet 是小波对象或名字， 可以看通过pywt.wavelist(kind='discrete') 的返回值，查看所有可用于离散变换的小波名字。

mode 是信号边界扩展模式，在pywt.Modes.modes 中列出了所有支持的模式，['zero', 'constant', 'symmetric', 'periodic','smooth', 'periodization', 'reflect', 'antisymmetric', 'antireflect']；

axis：整数类型，可选，用于计算DWT 的轴。如果没有给出，则使用最后一个轴。

返回值：(cA, cD)，元组，近似（低频信息）和细节系数（高频信息）

例：

7 多级一维离散小波变换

pywt.wavedec(data, wavelet, mode='symmetric',level=None, axis=-1)

data、wavelet、mode、axis 参数同上边的dwt 函数。

level：整型，可选，分解级别（必须大于等于0）。如果级别为“无”（默认），则将使用dwt_max_level 函数计算该级别。

返回值：列表，[cA_n, cD_n, cD_n-1, …, cD2, cD1]，系数数组的有序列表，其中n 表示分解级别。结果的第一个元素（cA_n）是近似系数数组，下面的元素（cD_n-cD_1）是细节系数数组。

低频部分代表了趋势，也叫近似信号；高频部分代表了噪声，也叫细节信号。

通常小波分解的频率与采样率有关，若n 层分解，则各个频率段的小小为fs/2/2n。如采样率为1000Hz，则根据奈奎斯特采样定理，最高频率为500Hz。做3 层分解，则：

一阶细节的频率为250-500Hz，一阶近似的频率小于250Hz；

二阶阶细节的频率为125-250Hz，二阶近似的频率小于125Hz；

三阶细节的频率为62.5-125Hz，三阶近似的频率小于62.5Hz；

例：

返回：ndarray，从系数中重建数据的一维数组。可用于显示分解各层的波形。

例：

8 小波包分解

小波包分析为信号提供了一种更加精细的分析方法，它将频带进行多层次分解，对多级分析中没有细分的高频部分进一步分解。并能根据被分析信号的特征，自适应的选择相应的频带，使之与信号频谱相匹配，即：既对低频部分进行分解，也对高频部分进行分解。

针对一维小波包分解，主要用的到是类pywt.Wavelet Packet，其构造方法为：

__init__(data, wavelet[, mode='symmetric'[,maxlevel=None[, axis=-1]]])

参数：data：与节点相关联的数据，或者说是待分析数据

wavelet：要使用的小波名字

mode：信号边界延拓模式，同前，默认为'symmetric'

maxlevel：最大允许的分解级别。如果未指定，则使用pywt.dwt_max_level（）基于小波和数据长度来计算。如未指定，可以调用maxlevel 属性获得。

axis：要转换的阵列的轴， 默认为-1

该类的常用方法和属性

data：与节点相关联的数据

get_leaf_nodes([decompose=False])

遍历分解树并返回叶子节点（没有任何子节点的节点）。

参数：decompose：如果分解为True，则该方法将尝试将树分解为最大级别。

get_level(level[, order="natural"[, decompose=True]])

获得给定的分解级别的节点。

参数：level：指定分解的级别数

order：指定节点的顺序，自然顺序（'natural'），频率顺序 （'freq'），默认为自然顺序。

decompose：如果设置为True，则该方法将尝试将数据分解为指定的级别。

例子：

ecg = pywt.data.ecg() #pywt 提供的例子数据

wp_ecg = pywt.WaveletPacket(ecg, "db5", maxlevel=3) #

使用'db5'小波分解，最大层数为3

wp_ecg.maxlevel #该值为3

leaf_nodes = wp_ecg.get_leaf_nodes(True) #所有的叶子节点

print([node.path for node in leaf_nodes]) #打印叶子节点的内容，显示为 ['aaa', 'aad', 'ada', 'add', 'daa', 'dad', 'dda', 'ddd']

aaa = wp_ecg['aaa'].data #获取aaa 节点的数据

示例数据集

在pywt 中提供了一些事例数据集，位于pywt.data 模块中，针对一维数据主要有：

pywt.data.ecg() ：心电图（ECG）波形（1024 个样本）

pywt.data.nino()：海面温度（264 个样本）

pywt.data.demo_signal(name='Bumps', n=None)：各种合成一维测试信号，

参数：name：要生成的测试信号类型（名称不区分大小写）。如果名称设置为“list”，则返回可访问测试函数的列表。

n：整型或者None，测试信号的长度。除具有固定长度的“Gabor”和“sineoneoverx”外，应为所有测试信号提供该参数。

返回：np.ndarray

9 NumPy简介

在前文中提到的类数组类型，既包括Python的数组类型，也包括NumPy 库的ndarray 类型。NumPy 是Python 语言的一个扩展程序库。支持高阶大量的维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库。

NumPy 提供了与MATLAB 相似的功能与操作方式，因为两者皆为解释型语言，并且都可以让用户在针对数组或矩阵运算时提供较标量运算更快的性能。两者相较之下，MATLAB 提供了大量的扩展工具箱（例如Simulink）；而NumPy 则是根基于Python 语言之上。此外NumPy 也可以结合其它的Python 扩展库，例如SciPy 以及Matplotlib 等

NumPy 的核心功能是ndarray（即n-dimensional array，多维数组）数据结构。这是一个表示多维度、同质并且固定大小的数组对象。