傅里叶变换到小波变换

傅立叶变换是最入门的，也是最先了解的，通过傅立叶变换，了解缺点，改进，慢慢的就成了小波变换。从傅立叶变换、短时傅立叶变换，小波变换等等，还有EMD变换。当然，其中会看到很多的名词，例如，内积，基，归一化正交，投影，Hilbert空间，多分辨率，父小波，母小波，这些不同的名词也是学习小波路上的标志牌。

傅里叶变换

信号处理最常用的工具就是傅里叶变换，但我们大多只是作为一个工具，会用就行了。但是深入理解其原理，才能更好的应用该工具。

背景

傅立叶是一位法国数学家和物理学家的名字，英语原名是Jean Baptiste Joseph Fourier(1768-1830)。 Fourier对热传递很感兴趣，于1807年在法国科学学会上发表了一篇论文，运用正弦曲线来描述温度分布，论文里有个在当时颇具争议性的命题：任何连续周期信号可以由一组适当的正弦曲线组合而成。当时审查这个论文的人，其中有两位是历史上著名的数学家拉格朗日(Joseph Louis Lagrange, 1736-1813)和拉普拉斯(Pierre Simon deLaplace, 1749-1827)，当拉普拉斯和其他审查者投票通过并要发表这个论文时，拉格朗日坚决反对，在近50年的时间里，拉格朗日坚持认为傅立叶的方法无法表示带有棱角的信号，如在方波中出现非连续变化斜率。法国科学学会屈服于拉格朗日的权威，拒绝了傅立叶的工作，幸运的是，傅立叶还有其它事情可忙，他参加了政治运动，随拿破仑远征埃及，法国大革命后因为怕被推上断头台而一直在逃难。直到拉格朗日死后15年这个论文才被发表出来。

那么，谁是对的呢？
拉格朗日是对的：正弦曲线无法组合成一个带有棱角的信号。但是，我们可以用正弦曲线来非常逼近地表示它，逼近到两种表示方法不存在能量差别，这个角度来说傅立叶是对的。

意义

傅里叶变换：任何连续周期信号都可以由一组适当的正弦曲线组合而成。
这个变换有什么意义呢？为什么要用正弦曲线的组合来代替原来的曲线呢？用方波或三角波来代替可以吗？
用方波或三角波来代替，分解信号的方法是无穷的，但分解信号的目的是为了更加简单地处理原来的信号。用正余弦来表示原信号会更加简单，因为正余弦拥有其他信号所不具备的性质：正弦曲线保真度。一个正弦曲线信号输入后，输出的仍是正弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的，且只有正弦曲线才拥有这样的性质，正因如此我们才不用方波或三角波来表示。

数学角度：
任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。

物理角度：
它其实是帮助我们改变传统的时间域分析信号的方法转到从频率域分析问题的思维。

下面的一幅立体图形可以帮助我们更好得理解这种角度的转换：

有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。

快速傅里叶变换FFT

首先，按照被变换的输入信号类型不同，傅立叶变换可以分为 4种类型：

• 非周期性连续信号傅立叶变换（Fourier Transform）
• 周期性连续信号傅立叶级数(Fourier Series)
• 非周期性离散信号离散时域傅立叶变换（Discrete Time Fourier Transform）
• 周期性离散信号离散傅立叶变换(Discrete Fourier Transform)

下面是四种原信号图例：

傅里叶级数，在时域是一个周期且连续的函数，而在频域是一个非周期离散的函数。
傅里叶变换，则是将一个时域非周期的连续信号，转换为一个在频域非周期的连续信号。

这里我们要讨论是离散信号，对于连续信号我们不作讨论，因为计算机只能处理离散的数值信号，我们的最终目的是运用计算机来处理信号的。所以对于离散信号的变换只有离散傅立叶变换（DFT）才能被适用，对于计算机来说只有离散的和有限长度的数据才能被处理，对于其它的变换类型只有在数学演算中才能用到，在计算机面前我们只能用DFT方法，我们要讨论的FFT也只不过是DFT的一种快速的算法。

DFT的运算过程：

FFT的过程大大简化了在计算机中进行DFT的过程，简单来说，如果原来计算DFT的复杂度是NN次运算（N代表输入采样点的数量），进行FFT的运算复杂度是Nlg10(N)，因此，计算一个1,000采样点的DFT，使用FFT算法只需要计算3,000次，而常规的DFT算法需要计算1,000,000次！
典型的时域2分裂算法图示如下：

做N点FFT之后，可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。FFT运算量：Nlog2N（2为对数的底）

如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，相应的采样点也为原来2倍，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即延长采样时间，所以频率分辨率和采样时间是倒数关系，就是说，要想分辨出频率间隔越小的频率（频率分辨率越高），采样时间越长越好。

存在问题

第一个就是傅立叶变换是整个时域，所以没有局部特征，这个也是他的基函数决定的看图，同时如果在时域张有了突变，那么在频域就需要大量的三角波去拟合，这也是傅立叶变换性质决定的。

第二个就是面对非平稳信号，傅立叶变换可以看到由哪些频域组成，但是不知道各成分对应的时刻是什么，也就是没有时频分析，看不出来信号频域随着时间变换的情况，反过来说就是，一个的频图对应好几个时域图，不知道是哪个，这个在实际应用中就不好了，看图：

做FFT后，我们发现这三个时域上有巨大差异的信号，频谱（幅值谱）却非常一致。尤其是下边两个非平稳信号，我们从频谱上无法区分它们，因为它们包含的四个频率的信号的成分确实是一样的，只是出现的先后顺序不同。
它只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。

总结：
1.不可处理非平稳信号；
2.基函数固定，基是正交的，分析重建用的基一样；
3.时频分辨率固定，不能刻画时间域上信号的局部特性；

短时傅里叶变换STFT

为了解决普通傅里叶的缺陷，短时间傅里叶变换出现了。在使用傅里叶变换前对原始信号通过滑动窗口的形式将信号分成几个相同长度。

小波变换

首先说说什么是波？就是在时间域或者空间域的震荡方程，比如傅里叶用的正弦波， 正弦波能量无穷，同样的幅度在整个无穷大区间里面震荡：

而小波的“小”就是一种能量在时域非常集中的波，即能量有限且集中在某一点附近：

小波还有个特点，可长可短可胖可瘦（伸缩平移）。小波有什么好处呢？它对于瞬时时变的信号很有用。
小波基函数会伸缩、会平移（其实是两个正交基的分解）。缩得窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。

某一个尺度（宽窄）下乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。那么我们就知道信号包含该频率的成分的多少。

分析动态频谱最好的方法就是使用小波变换。小波变换在频域和时域上都有很好的表现。它不仅能够告诉我们在信号种有哪些频率出现，而且能够告诉我们在信号的什么时候出现。这是通过不同的伸缩变换实现的。首先，我们看一个大尺度/窗口的信号，并分析“大”特征，然后我们看一个小尺度的信号，以便分析更小的特征。图2显示了不同方法的时间分辨率和频率分辨率。

基（basis）

傅立叶变换和小波变换，都会听到分解和重构，其中这个就是根本，因为他们的变化都是将信号看成由若干个东西组成的，而且这些东西能够处理还原成比原来更好的信号。那怎么分解呢？那就需要一个分解的量，也就是常说的基。

基（basis）：在线性代数中我们都曾学过，它表示空间中的一系列线性独立的向量，空间中的其他向量都可以通过它的线性组合表示。

对于傅立叶变换的基是不同频率的正弦曲线，所以傅立叶变换是把信号波分解成不同频率的正弦波的叠加和，而对于小波变换就是把一个信号分解成一系列的小波，那小波是什么呢，定义中告诉我们小波，因为这个小波实在是太多，一个是种类多，还有就是同一种小波还可以尺度变换，但是小波在整个时间范围的幅度*均值是0，具有有限的持续时间和突变的频率和振幅，可以是不规则，也可以是不对称,很明显正弦波就不是小波。

基本小波（母小波）：
小波函数需满足如下的的条件：
• 有限的能量
• 均值为零
• 可以正交也可以不正交
• 可以是双正交的，也可以不是
• 可以是对称的，也可以不是
• 可以是负数，也可以是实数
• 如果它是复数，通常分为表示振幅的实部和表示相位的虚部
• 被归一化为单位能量

第一个条件意味着它在时间和频率上具有局部性的;它是可积的，小波与信号的内积始终存在。第二个条件意味着小波在时域中均值为零，在时域中频率为零时均值为零。这对于保证小波变换的可积性和计算小波变换的逆是必要的。

常见为小波有：

import pywt
from matplotlib import pyplot as plt
#print(pywt.families(short=False))
discrete_wavelets = ['db5', 'sym5', 'coif5', 'bior2.4']
continuous_wavelets = ['mexh', 'morl', 'cgau5', 'gaus5']
​
list_list_wavelets = [discrete_wavelets, continuous_wavelets]
list_funcs = [pywt.Wavelet, pywt.ContinuousWavelet]
​
fig, axarr = plt.subplots(nrows=2, ncols=4, figsize=(16,8))
for ii, list_wavelets in enumerate(list_list_wavelets):
    func = list_funcs[ii]
    row_no = ii
    for col_no, waveletname in enumerate(list_wavelets):
        wavelet = func(waveletname)
        family_name = wavelet.family_name
        biorthogonal = wavelet.biorthogonal
        orthogonal = wavelet.orthogonal
        symmetry = wavelet.symmetry
        if ii == 0:
            _ = wavelet.wavefun()
            wavelet_function = _[0]
            x_values = _[-1]
        else:
            wavelet_function, x_values = wavelet.wavefun()
        if col_no == 0 and ii == 0:
            axarr[row_no, col_no].set_ylabel("Discrete Wavelets", fontsize=16)
        if col_no == 0 and ii == 1:
            axarr[row_no, col_no].set_ylabel("Continuous Wavelets", fontsize=16)
        axarr[row_no, col_no].set_title("{}".format(family_name), fontsize=16)
        axarr[row_no, col_no].plot(x_values, wavelet_function)
        axarr[row_no, col_no].set_yticks([])
        axarr[row_no, col_no].set_yticklabels([])
​
plt.tight_layout()
plt.show()

第一行包含四个离散小波，第二行包含四个连续小波。可在 Wavelet Properties Browser （http://wavelets.pybytes.com/）查看不同的小波的情况。

在每个小波族中，可以有许多不同的小波子类别属于这个族。你可以通过系数的数量(消失矩的数量)和分解的级别来区分小波的不同子类别。
这是在下面为一个族的小波称为“Daubechies”。

db_wavelets = pywt.wavelist('db')[:5]
print(db_wavelets)
*** ['db1', 'db2', 'db3', 'db4', 'db5']
​
fig, axarr = plt.subplots(ncols=5, nrows=5, figsize=(20,16))
fig.suptitle('Daubechies family of wavelets', fontsize=16)
for col_no, waveletname in enumerate(db_wavelets):
    wavelet = pywt.Wavelet(waveletname)
    no_moments = wavelet.vanishing_moments_psi
    family_name = wavelet.family_name
    for row_no, level in enumerate(range(1,6)):
        wavelet_function, scaling_function, x_values = wavelet.wavefun(level = level)
        axarr[row_no, col_no].set_title("{} - level {}\n{} vanishing moments\n{} samples".format(
            waveletname, level, no_moments, len(x_values)), loc='left')
        axarr[row_no, col_no].plot(x_values, wavelet_function, 'bD--')
        axarr[row_no, col_no].set_yticks([])
        axarr[row_no, col_no].set_yticklabels([])
plt.tight_layout()
plt.subplots_adjust(top=0.9)
plt.show()

Daubechies，一般音译为“多贝西”。上图中我们看到多贝西小波族，第一列中，我们可以看到一阶(db1)的Daubechies小波，在二阶(db2)的第二列中，一直到第五列中的第五阶。小波包含多达20阶(db20)的Daubechies小波。
阶数表示消失力矩的个数。所以db3有3个消失力矩db5有5个消失力矩。消去矩的个数与小波的逼近阶数和平滑度有关。如果一个小波有p个消失矩，它可以近似p - 1次多项式。
在选择小波时，我们还可以指出分解的级别。默认情况下，PyWavelets为输入信号选择可能的最大分解级别。分解的最大级别(参见pywt.dwt_max_level())取决于输入信号长度和小波的长度。
可以看出，随着消失矩数的增加，小波的多项式阶数增加，小波变得更加平滑。随着分解程度的增加，这个小波的样本数以增加的形式表示。

母小波经过平移和伸缩，可得到一组函数，称之为小波基函数，简称小波基。

以下列出的15种小波基是Matlab中支持的15种：

小波基的选择：https://blog.csdn.net/heifan2014/article/details/72530858

内积

内积用来刻画两个向量的夹角。在Hilbert空间，当内积为0时，两个向量正交，若g为Hilbert空间里的正交基的时候，内积为f向基上的正交投影。
假设有两个函数f和g，两个向量求内积，就是把它们相同的位置上对应的点的乘积做一个累加，就是x的每一个点，对应的f(x)和g(x)做乘积，再累加。但f和g是无限函数，x又是一个连续值，所以表示成积分：

小波变换

如前边所说，小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。效果如下图：
现在来看一下小波公式：

从公式可以看出，不同于傅里叶变换，变量只有频率ω，小波变换有两个变量：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函数的伸缩，移量 τ控制小波函数的移。尺度就对应于频率（反比），*移量 τ就对应于时间。如下图

当伸缩、移到这么一种重合情况时，也会相乘得到一个大的值。这时候和傅里叶变换不同的是，这不仅可以知道信号有这样频率的成分，而且知道它在时域上存在的具体位置。
而当我们在每个尺度下都平移着和信号乘过一遍后，我们就知道信号在每个位置都包含哪些频率成分。
（1） 解决了局部性

（2）解决时频分析
做傅里叶变换只能得到一个频谱，做小波变换却可以得到一个时频谱！

既然小波是在时间上局域的，我们可以叠加不同时间段的小波。我们从信号开始的时间点慢慢想信号结尾移动小波。这个方法也是所说的卷积。在对原始（母）小波做完这些，我们可以改变小波尺寸，使得它变大然后继续上面的操作。这个过程如下图所示。

一个直观例子（和傅里叶对比）：
小波变换对不管一维还是高维的大部分信号都能很好表示。这个和傅里叶级数有很大区别。傅里叶级数擅长把一维的、类三角波连续信号映射到一维系数序列，但对于突变信号或任何高维的非三角波则无能为力。
假设有一个直流信号：

用傅里叶将其展开：会发现形式非常简单，只有一个级数的系数不是0，其他所有级数系数都是0。

再看另一个信号：直流信号上增加一个突变

我们再用傅里叶将其展开：此时所有傅里叶级数都是非0，为甚么？因为傅里叶必须用三角波来展开信号，对于这种突变信号，傅里叶就会用大量三角波去拟合，就像这样：

但再怎么拟合，都只是近似，这也就是吉伯斯Gibbs现象。Gibbs是由于展开式再间断点邻域不能均匀收敛所引起的，即使N趋于无穷大，这一现象也依然存在。
那我们用小波展开这个突变信号：

可以看到，小波基在突变部分有值，其他部分对于的级数都为0。比如用三层小波展开，那么只有3个系数不为0，其他都为0，这是由于小波基的特殊性，任何小波和常量函数的内积都趋近于0。所以，这就要求小波基在一个周期的积分为0。

代码实践

Python中有几个库提供了小波变换的功能，主要包括PyWavelets（pywt）、pycwt、和scipy的小波模块。每个库都有其特点和适用场景，以下是它们的主要区别和用途：

1.PyWavelets（pywt）
特点： PyWavelets是一个开源的小波变换库，提供了广泛的离散小波变换（DWT）、逆小波变换（IDWT）、以及连续小波变换（CWT）功能。它支持多种小波和多级变换。
适用场景： 非常适合于需要进行信号分析、信号压缩、去噪或特征提取的应用。由于其提供了丰富的小波选项和灵活的多级分析功能，它被广泛用于科学研究和工业应用。

import numpy as np
import pywt

# 创建一个模拟信号
# 这里使用简单的正弦波叠加随机噪声作为示例
signal_length = 2**11  # 信号长度，确保能够进行10层分解
signal = np.sin(np.linspace(0, 8 * np.pi, signal_length)) + np.random.normal(0, 0.5, signal_length)

# 选择小波类型并执行10层小波变换
# 'db1'是Daubechies小波家族中最简单的小波
coeffs = pywt.wavedec(signal, 'db1', level=10)

# 输出每层的细节和近似系数的长度，确认变换正确执行
for i, coeff in enumerate(coeffs):
    print(f"层 {i}: 长度 {len(coeff)}")

# 注意：这里的输出长度展示了每一层分解的系数长度
# 第一层是最低频率的近似，之后每层的细节系数代表了不同的频率范围

2.pycwt
特点： pycwt专注于提供连续小波变换（CWT）的实现。它允许用户详细地分析信号的局部特征，以及信号频率随时间的变化。
适用场景： pycwt特别适合于分析和处理具有复杂频率行为的信号，如地震信号、生物医学信号等。由于连续小波变换提供了信号的详细时频表示，这个库在需要精细时频分析的领域中特别有用。

import numpy as np
import pycwt as wavelet
import matplotlib.pyplot as plt

# 创建一个模拟信号，例如简单的正弦波加上一些噪声
N = 1000
t = np.linspace(0, 1, N)
signal = np.cos(2 * np.pi * 7 * t) + np.random.normal(0, 0.5, N)

# 定义时间序列的采样间隔
dt = t[1] - t[0]

# 使用Morlet小波进行连续小波变换
wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(signal, dt, dj=1/12, s0=-1, J=-1, wavelet='morlet')

# 绘制小波变换的功率谱
power = np.abs(wave) ** 2
plt.contourf(t, freqs, power, 100, cmap='jet')
plt.colorbar(label='Power')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (sec)')
plt.title('Wavelet Power Spectrum')
plt.show()

要注意的是，pycwt专注于连续小波变换，这与离散小波变换（DWT）是有区别的。连续小波变换提供了一种不同的方式来分析信号，在某些应用中可能更加有用，比如在处理非平稳信号时。

3.Scipy
特点： scipy是一个广泛用于科学计算的库，其中包括了一个小波模块。scipy的小波部分提供了一些基础的小波变换功能，虽然不如PyWavelets那样全面，但足以满足一些基本应用需求。
适用场景： 适合于需要在科学计算中使用小波变换，但对小波变换的高级功能需求不大的用户。由于scipy是一个综合性的科学计算库，使用它的小波模块可以方便地与scipy提供的其他数值计算和信号处理功能结合使用。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import cwt, ricker
from scipy import signal

# 创建信号：一个简单的正弦波加上噪声
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 7 * t) + np.random.normal(0, 0.5, t.shape)

# 定义用于CWT的小波函数和尺度
widths = np.arange(1, 31)
wavelet = ricker  # 使用Ricker小波（也称为Mexican hat小波）

# 执行连续小波变换
cwtmatr = cwt(signal, wavelet, widths)

# 绘制结果
plt.imshow(cwtmatr, extent=[0, 1, 1, 31], cmap='PRGn', aspect='auto',
           vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
plt.colorbar(label='Magnitude')
plt.xlabel('Time')
plt.ylabel('Scale')
plt.title('Continuous Wavelet Transform (CWT) of the signal')
plt.show()

CWT和DWT的区别

在计算机实际处理信号时，连续小波变换（CWT）和离散小波变换（DWT）的主要区别在于它们的应用目标、计算方法、以及结果的解释方式。

1. 应用目标
   CWT：CWT提供了一种在连续尺度上分析信号的方法，允许用户检查信号在不同频率（或尺度）和时间点的行为。它特别适合于分析局部特征和非平稳信号，例如瞬时频率随时间变化的信号。
   DWT：DWT是对信号进行多尺度分析的有效工具，通常用于信号压缩和去噪。与CWT相比，DWT在尺度和位置上都是离散的，这使得它在处理大量数据时更为高效和简洁。
2. 计算方法
   CWT：在计算上，CWT通过将信号与一系列缩放和平移后的母小波进行卷积来实现。小波的伸缩和平移是连续的，因此可以产生非常精细的尺度分辨率。结果是一个二维表示，显示了信号的时间-频率特征。
   DWT：相反，DWT使用一套离散的小波基函数来逼近信号，在伸缩和平移上是离散的，通常按照2的幂次来缩放（称为二进制尺度）和平移小波函数。DWT将信号分解成一系列近似和细节系数，每一次分解都增加了分析的尺度（频率范围变低）。
3. 结果解释
   CWT结果：CWT的结果是一个时间-尺度（或时间-频率）平面上的连续分布，可以用于直观地分析信号的局部变化。它特别适合于识别和分析信号中的瞬态特征，如突发事件、边缘或其他非周期性变化。
   DWT结果：DWT的结果是一组逐步下采样的近似系数和一系列细节系数，每一级分解都对应于信号在不同频率层次上的表现。这种多层次的结构使得DWT非常适合于图像压缩、去噪和数据压缩。
4. 计算效率和实现复杂性
   CWT：由于CWT的连续性质，它在计算上通常比DWT更为复杂和计算量更大。尽管现代算法和优化手段已经大大提高了CWT的计算效率，但它仍然需要更多的计算资源。
   DWT：DWT的计算过程更简洁高效，特别是对于大规模数据处理。这得益于其离散的特性和能够使用快速算法（如快速小波变换，FWT）进行优化。

总之，CWT和DWT在信号处理中各有优势，选择哪一种方法取决于具体的应用需求、计算资源以及预期的结果类型。