前言

之前我们介绍了基于STFT时频分析方法来提取转速信息，但也提到了该方法的不足之处：由于采用了固定窗口，难以平衡时间分辨率和频率分辨率，窗口太小，时间分辨率高，能够捕捉信号的快速变化，但是频率分辨率低，无法区分相近的频率成分，窗口太大，频率分辨率高，能够区分相近的频率成分，但是时间分辨率低，无法捕捉信号的快速变化。在实际应用中，选择某个窗函数，希望其时-频窗形状是自适应变化的，对低频信号，其窗口形状自动变得扁平；对高频信号，其窗口形状自动变得瘦长。
为了解决这个问题，我们后面将基于小波分析与脊线提取工具来进行转速计算工作，为了更容易理解，本篇文章先介绍小波分析。

小波分析原理

作为目前首选的时频分析方法，小波分析拥有出色的自适应时频分析能力，我们先来看看其发展历史：

1. 起源：1909年Haar提出了最早的类似小波的基函数（Haar小波），由于该小波不连续，因此没有引起人们的重视，直到70年代，Morlet（地球物理学家）为解决地震信号分析问题，提出了Morlet小波，并建立了分解重构信号的方法，标志着现代小波理论的开端。
2. 理论奠基：1984年Alex Grossman与Morlet合作建立了连续小波变换的数学框架，1986年Meyer构造了一个无穷光滑的小波标准正交基Meyer小波，1989年Mallat提出了多分辨率分析理论，统一了小波与滤波器组设计，奠定离散小波变换（DWT）基础，Ingrid Daubechies也在这几年构造了紧支撑正交小波族（如dbN系列），解决实际计算问题。
3. 应用扩展：1989年Mallat提出了快速小波变换，复杂度仅为O(N)，而后又逐渐在图像压缩、故障诊断等领域崭露头角。

那么为什么它具备如此强大的能力呢？关键就在于基函数的选择。傅里叶变换使用了正弦/余弦函数作为基函数，由于其是一种无限长的信号，不管在时间轴上如何移动，与原始时域信号做内积后得到的系数都不会变（傅里叶变换可以看作原始信号在以正弦/余弦函数为基础构建的向量空间里做投影，变换系数就是投影值），因此无法定位某个时间的频率信息。而小波分析使用了小波基函数，它是一种在时域上能量衰减的信号，我们以Morlet小波为例，它的数学表达式如下：
$$\psi (t)={\pi }^{-\frac{1}{4}}(e^{i{w}_{0}t}-e^{-\frac{w_0^2}{2}})e^{-\frac{t^2}{2}}$$
其中$w_0$是中心频率（通常取$w_0>=5$以确保小波满足零均值条件），${\pi }^{-\frac{1}{4}}$是归一化因子，保证小波的能量为1。而当$w_0>=5$的时候，第二项$e^{-}\frac{w_0^2}{2}$可以忽略，因此Morlet小波可以简化为：
$$\psi (t)={\pi }^{-\frac{1}{4}}{e}^{i{w}_{0}t}{e}^{-\frac{t^2}{2}}$$
这也是Morlet小波的常用形式，其曲线图如下：

绘制的代码如下：

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


plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
plt.rcParams['axes.unicode_minus'] = False   # 步骤二（解决坐标轴负数的负号显示问题）


def morlet_wavelet(t, omega0=6):
    """
    生成Morlet小波函数
    :param t: 时间轴
    :param omega0: 中心频率
    :return: Morlet小波
    """
    return np.pi**(-1/4) * np.exp(1j * omega0 * t) * np.exp(-t**2 / 2)


# 定义时间轴
t = np.linspace(-4, 4, 1000)

# 计算Morlet小波
psi = morlet_wavelet(t, omega0=6)

# 可视化
plt.figure(figsize=(10, 6))
plt.plot(t, np.real(psi))
plt.title('Morlet小波图')
plt.xlabel('时间')
plt.ylabel('幅值')
plt.legend()
plt.show()

当然我们也可以调用pywt库来绘制，代码如下：

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


plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
plt.rcParams['axes.unicode_minus'] = False   # 步骤二（解决坐标轴负数的负号显示问题）


wavelet = pywt.ContinuousWavelet('morl')
[psi, xval] = wavelet.wavefun(length=1000)
plt.plot(xval, psi)
plt.title("Morlet小波图")
plt.xlabel('时间')
plt.ylabel('幅值')
plt.legend()
plt.show()

可以看到，二者的形状基本一致，而且中心频率处能量最大，越向两边走能量越小，直到趋近于0。试想一下，若是使用这样的信号在时间轴上不断移动着与原始时域信号做内积运算，其值越大，则说明该时间处的频率与该小波频率越匹配，如果我们使用多个不同频率的小波函数，就可以实现计算不同时间上不同频率能量的分布了。
我们再来看看连续小波变换的计算公式，就能够更深刻地理解上述表达：
$$W_f(a,b)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{\infty }f(t){\psi }^{\ast }(\frac{t-b}{a})dt$$
其中$f(t)$是原始时域信号，${\psi }^{\ast }$是小波函数的复共轭，a是尺度参数（控制频率，当a较小时，频窗中心自动地调整到较高的频率中的位置，且时频窗形状自动地变为“瘦窄”状，当a较大时，时-频窗形状自动地变为“扁平”状），b是位移参数（控制时间），通过调整a和b，可以让小波函数具有不同的频率和时域分辨率，整个公式就代表信号与小波函数的复共轭在整个时间轴上的内积运算。
如果小波函数是复值，那么复共轭可以定义为：
$${\psi }^{\ast }(t)=Re(\psi (t))-i\ast Im(\psi (t))$$
如果小波函数是实值，那么复共轭就是其本身。
在小波变换中，尺度参数和频率的转换公式如下：
$$scale\ast f=Fs\ast wcf$$
其中$Fs$是信号的采样频率，$wcf$是小波的中心频率。

Python实现

有了上面的理论基础，我们可以用一个简单的例子看看，如何基于Morlet小波实现连续小波变换。假定原始时域信号为两个不同频率的正弦函数叠加，整体代码如下所示:

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


plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
plt.rcParams['axes.unicode_minus'] = False   # 步骤二（解决坐标轴负数的负号显示问题）


# 生成离散信号
N = 1000
t = np.linspace(0, 1, N)
f = np.sin(2 * np.pi * 10 * t) + np.sin(2 * np.pi * 20 * t)

# 计算CWT
scales = np.arange(1, 128)

# Morlet小波
coef, freqs = pywt.cwt(f, scales, 'morl')  

# 绘制时频图
plt.figure(figsize=(12, 6))
plt.contourf(t, freqs, abs(coef))
plt.colorbar(label='幅值')
plt.title('连续小波变换时频分析')
plt.xlabel('时间')
plt.ylabel('频率')
plt.show()

这里可以看到在0.01和0.02两个频率上的能量最为凸显,和我们生成的正弦信号频率一致(由于采样频率为0.001,10正弦Hz信号对应的最终频率为0.01,20Hz正弦信号对应的最终频率为0.02)。

总结

本篇文章介绍了小波变换的基本原理和Python实现，其本质上就是通过一个衰减的小波基替代了无限长的正弦/余弦基，从而可以定位频率的时间信息。当然除了连续小波分析外，还有离散小波变换、小波包分解等诸多分支，留待后续介绍。