python连续小波分析CWT

转自：python连续小波分析

0 引言

我们学过内积，内积的物理含义：两个图形的相似性，若两个图形完全正交，则内积为0，若两个图形完全一样，则系数为1（相对值）。

小波变换的实质是：原信号与小波基函数的相似性。小波系数就是小波基函数与原信号相似的系数。（英文文献中是这样解释：The definition of wavelet transform shows that the wavelet analysis is a measure of similarity the basis functions(wavelets)and the original function .The coefficients caculated indicate how close the function is to the danghter wavelet at that particular scale )

连续小波变换：小波函数与原信号对应点相乘，再相加，得到对应点的小波变换系数，平移小波基函数，再计算小波函数与原信号对应点相乘，再相加，这样就得到一系列的小波系数。

对于离散小波变换（由于很多小波函数不是正交函数，因此需要一个尺度函数）所以，原信号函数可以分解成尺度函数和小波函数的线性组合，在这个函数中，尺度函数产生低频部分，小波函数产生高频部分。
（本段摘自CSDN博主「不负韶华T」：https://blog.csdn.net/heifan2014/article/details/72561758）

1 实例

1.1 结果图

上面为原始信号，下面为连续小波变换得到的时频数据

涉及到代码不会的知识点均可在下链接查找：
小波变换介绍

小波时频(尺度图)的绘制原理与实现

pywt模块使用说明

Python Matplotlib add_subplot 和 subplots_adjust详解及代码详细说明 配图片说明

matplotlib.pyplot contourf()函数的使用

1.2 代码

import numpy as np
import pywt
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pylab as plt

sampling_rate = 1024  # 采样频率
t = np.arange(0, 1.0, 1.0/sampling_rate)  # 0-1.0之间的数，步长为1.0/sampling_rate
f1 = 100  # 频率
f2 = 200
f3 = 300
data = np.piecewise(t, [t < 1, t < 0.8, t < 0.3], 
                    [lambda t: np.sin(2 * np.pi * f1 * t), 
                     lambda t: np.sin(2 * np.pi * f2 * t), 
                     lambda t: np.sin(2 * np.pi * f3 * t)]
                    )
wavename = "cgau8"  # 小波函数
totalscal = 256     # totalscal是对信号进行小波变换时所用尺度序列的长度(通常需要预先设定好)
fc = pywt.central_frequency(wavename)  # 计算小波函数的中心频率
cparam = 2 * fc * totalscal  # 常数c
scales = cparam/np.arange(totalscal, 1, -1)  # 为使转换后的频率序列是一等差序列，尺度序列必须取为这一形式（也即小波尺度）
[cwtmatr, frequencies] = pywt.cwt(data, scales, wavename, 1.0/sampling_rate)  # 连续小波变换模块

plt.figure(figsize=(8, 4))
plt.subplot(211)  # 第一整行
plt.plot(t, data)
plt.xlabel(u"time(s)")
plt.title(u"300Hz 200Hz 100Hz Time spectrum")
plt.subplot(212)  # 第二整行

plt.contourf(t, frequencies, abs(cwtmatr))
plt.ylabel(u"freq(Hz)")
plt.xlabel(u"time(s)")
plt.subplots_adjust(hspace=0.4) #调整边距和子图的间距 hspace为子图之间的空间保留的高度，平均轴高度的一部分
plt.show()

1.3 结果分析

• 小波变换需要提前设定 totalscal（对信号进行小波变换时所用尺度序列的长度），本文为 256，小波变换计算结果返回了 cwtamtr和frequencies，其结构分别为 (255, 1024)，(255,)，数据一共 1024 个

• cwtmatr ：小波变换的返回结果之一，包括实部和虚部，结构大小为 (255, 1024)

  • 在连续小波变换时，当选用小波母函数是复数而非实数时，则该小波变换就定义为连续复小波变换，变换后的系数为复数，即WTf=WTr+jWTi，其中WTr和WTi为复小波系数的实部和虚部。
  • 复小波变换就是用实部和虚部两个实小波同时对信号进行变换。复小波的实部和虚部通常是正交的，所以信号的复小波变换实际上是将该信号在两个互相正交的空间同时做实小波变换，相比实小波变换只能得到原信号一个空间中的信息，复小波同时提供了两个正交空间的信息。
  • 这样，复小波变换除了提供类似于实小波变换的幅值信息外，还可以提供两个正交空间中信号分量间特有的相位信息。可以得到的实部、虚部、幅值和相位四种基本信息进行适当的组合构成复小波变换的复合信息。这些特点有助于更准确地捕捉畸变信号和提供更丰富的信号细微特征差异，克服实小波无相位信息和平移敏感性的缺点。
    

• frequencies：小波变换返回另一结果，结构为 {ndarray: (255, )}，为计算出来的时频数据
  

2 cwt 使用介绍

pywt.cwt(data, scales, wavelet)
连续小波变换

• Parameters:
  • data : array_like， 信号数组
  • scales : 要使用的小波尺度（s）。
    可以用 *f = scale2frequency(wavelet, scale)/sampling_period 来确定物理频率大小。f的单位是赫兹，采样周期的单位为秒。
  • wavelet : Wavelet 对象或名字
  • sampling_period : float
    频率输出的采样周期。coefs的计算值与sampling_period的选择无关。scales不按抽样周期进行缩放。
  • axis: int, optional
    计算CWT的轴。如果不给出，则使用最后一个轴。
• Returns:
  • coefs : array_like，
    给定尺度和小波的输入信号的连续小波变换。第一维是scales，其余为data
  • frequencies : array_like。
    如果采样周期的单位是秒，那么频率的单位是赫兹。否则，假设采样周期为1。

3. 参考链接

[1] matplotlib.pyplot contourf()函数的使用 2018.11；