今天这篇写小波变换，这也是这个系列**最后一篇，他解决的问题同希尔伯特-黄变换类似，都是观察随时间变化的频谱。但貌似小波变换应用的更为普遍，比如图片压缩、地震波检测之类。像希尔伯特-黄变换一样，我只知道小波变换的基本原理，但对其推导过程以及更深层的原理并没有系统的学习过，但这并不妨碍我们使用他人写好的小波变换去分析信号，必要时再去学习也不迟。
为了便于理解代码，先粗略的讲讲小波变换的基本原理。由于没有系统学习过小波变换，只从网上翻阅了些粗浅的资料，有疏漏之处还请谅解与指出。
小波变换首先离不开小波，所谓小波，指的是有限长度，有一定尺度的的波。长度好理解，而尺度对应于频率，当频率高时，尺度就小，当频率低时，尺度就大。如下图所示，不同尺度下的小波形状。
我们知道，傅里叶变换是把一个时域信号展开为若干正弦函数的组合，而小波变换就是把一个时域信号展开为若干不同尺度、不同位置的小波组合[1]。
现在，假设我们有了一个小波，怎么通过小波发现时域信号频率与时间的关系呢？可分为以下几个步骤：
1、用这个小波与时域信号相乘[2]，然后看其结果，结果越大，表明相似性越高。而相似度的高低又可以表征这一段的时域信号是否具有某个频率。相似度为0，表示时域信号没有该频率。相似度大，表示时域信号有该频率。在小波变换中，衡量小波与时域信号相似度的参数称为小波系数（这里描述不严谨，只是做为一个门外汉来说便于理解而已）
2、移动这个小波，与下一段时域信号做类似1步骤的动作，求解小波系数，直到整个时域信号分析完成，这样就得到了在相同尺度下，不同位置的一组小波系数
3、使用不同尺度的小波，重复1-2步骤
由1-3步骤，我们就能得到一个小波系数矩阵，这个矩阵中，同一行的每个数据表示一个小波与每一段时域信号做运算得到的相同尺度下，不同位置的小波系数。不同行则表示不同尺度的小波与时域信号做运算得到一批小波系数。于是，根据小波系数矩阵，我们就知道在哪个时间点出现了哪些频率。
有了上述的准备，我们就可以用Python来做小波变换了，见下面代码。
# -*- coding: utf-8 -*-
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
fs = 200
t = np.arange(fs) / fs
#这里使用的是两个简单信号叠加，且频率不随时间变换，这么做是为了便于观察频谱图
s = np.cos(2*np.pi*5*t) + np.sin(2*np.pi*60*t)
#生成1-100Hz的频率，一共100个点(分辨率1Hz)
freq = np.linspace(1, fs/2, 100)
w = 30
#计算不同频率下的小波尺度，对于morlet2小波，他的尺度计算公式a=w*fs/(2*f*np.pi)，其中
#f为时域频率，fs为采样率，参考[4]
a = w*fs / (2*freq*np.pi)
#signal.cwt()函数对信号做连续小波分析，原型参考[3]
#signal.morlet2是一个小波函数，原型参考[4]
#a、w是小波函数signal.morlet2需要的参数，这个函数会按照a、w生成一系列不同尺度的小波
#signal.cwt()函数返回值cwtm就是小波系数矩阵
#前面说过，系数矩阵中，每一列是相同尺度下、不同位置的小波系数，每一行是不同尺度下的小波系数
#第一行对应最高频率（这里是100Hz）小波系数，矩阵最后一行为最低频率(这里是1Hz)的小波系数
cwtm = signal.cwt(s, signal.morlet2, a, w=w)
#画出小波系数矩阵坐标图，横轴为时间，纵轴为频率
#np.abs(cwtm)表示求cwtm的模（他是复数），然后把不同大小的模值与cmap指定的颜色集一一映射
plot1 = plt.pcolormesh(t, freq, np.abs(cwtm), cmap='viridis', shading='gouraud')
#plt.colorbar（）画出颜色对应的模值，便于观察
plt.colorbar(plot1)
plt.ylabel('Frequency [Hz]')
plt.xlabel('t(s)')
plt.text(0,5,'5Hz',color='r')
plt.text(0,60,'60Hz',color='r')
plt.show()
频谱图如下：
可以看到，黄、绿颜色值表示的值是最大的（3-5，参考上图右边的颜色值指示器），对应的小波系数也是最大的，而频谱图中，恰好是5Hz、60Hz处这两个颜色最深。因此，利用cwt（）函数，我们就可以分析出这个时域信号频率与时间的关系了。
在scipy的signal模块中，小波分析的功能并不强大，这里仅仅是了解一下小波变换到底是在做什么，要进使用小波变换去解决实际问题，可以看看[5]，这里面有完备的小波变换库。如果你研究了，也希望你能写写关于这个库的使用方法并告诉我，让我们也学习学习。
[1]https://www.zhihu.com/question/22864189/answer/40772083，来自“bellaxia”的回答
[2]https://www.zhihu.com/question/22864189/answer/40772083，来自“咚懂咚懂咚”的回答
[3]https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.signal.cwt.html#scipy.signal.cwt
[4]https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.signal.morlet2.html#scipy.signal.morlet2
[5]https://github.com/PyWavelets/pywt
|
楼主
标题置顶
标题高亮
