- # -*- 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
楼主
