从傅里叶变换到小波变换，并不是一个完全抽象的东西，其实也可以讲得很形象的。小波变换有着明确的物理意义，如果我们从它被提出时所面对的问题来看，思路就会清晰不少。本文按照傅里叶->短时傅里叶变换->小波变换的顺序，讲解一下为什么会出现小波以及小波究竟是怎样的思路？
一、傅里叶变换
傅里叶变换可谓有史以来最伟大、最深刻的数学发现之一。但不幸的是，初次见它的公式似乎一般很难理解其中的内涵。例如，对满足狄里赫利(Dirichlet)条件的周期信号做傅里叶变换可以得到一组傅里叶级数，其可以表示为：
说到数学中所谓的“变换”，其实质是变换看待问题的角度，而并非变换问题自身。（同理，有兴趣也可以了解下另外一个变换：拉普拉斯变换中的S是个什么鬼？。）
为了更好地理解“变换”的概念，现在给出一个Better Explained上关于“变换”的简单比喻：
假如有一种饮料叫做：橘子香蕉牛奶冰沙，你非常爱喝，于是想要知道如何自制橘子香蕉牛奶冰沙，以满足日常的饮用需求！
你知道要想做出和市场上口感一样的冰沙，最重要的是要知道原材料的种类和比例。那么，这个从一杯橘子香蕉牛奶冰沙中获取原材料和比例的过程也就可以类比于傅里叶变换的过程了。这里姑且称之为“橘子香蕉牛奶冰沙的傅里叶变换”。
以下有几个问题先帮你整理下思路：
1、 “橘子香蕉牛奶冰沙的傅里叶变换”是做什么的？
答：得到一杯水果冰沙后，找出其中包含的各种成分及其性质。
2、怎么做？
答：让冰沙通过某种“过滤器”，以便提取出其中的每种成分。
3、为什么这样做？
答：各个单一的组成成分比冰沙本身更容易分析，比较和修改。
4、如何自制橘子香蕉牛奶冰沙？
答：将过滤得到的各个组成成分按分析所得的比例混合即可。
示意图如上所示，如果倒入了9个单位量的水果冰沙：
经过“香蕉过滤器”，得到1个单位量的香蕉；
经过“橘子过滤器”，得到2个单位量的橘子；
经过“牛奶过滤器”，得到3个单位量的牛奶；
经过“冰沙过滤器”，得到3个单位量的冰沙。
这样，我们就获得了制造冰沙的“配方”，只需将各个成分按过滤得到的比例加以混合，就可以得到和市场上一模一样的橘子香蕉牛奶冰沙了！
从消费者的角度，看到的是包含有“香蕉”、“橘子”、“牛奶”和“冰沙”的水果冰沙；而从冰沙的制造者角度来说，关心的是制作冰沙的配方，即成分和比例！
“橘子香蕉牛奶冰沙的傅里叶变换”将我们的视角从消费者转向生产者，从“我有水果冰沙“转向“水果冰沙是怎么制作的？”，这就是两种不同的角度，而实现这两种角度切换的就是上图中的“过滤器”（也即“变换”）。
傅里叶级数的直观表达
傅里叶认为“任何”周期信号都可以表示为一系列成“谐波关系”的正弦信号的叠加。
正弦函数无需多言，大家都清楚，但为了直观化表达，这里将正弦(或余弦)信号同时以两种运动形式来表示，分别是：
(1)一种是以时间为横轴、位移为纵轴，某一点的往复运动，也就是通常所说的正弦波，或者说是振荡信号；
(2)某一点绕另一点的匀速圆周运动。两种情况综合起来为下图所示。正弦波就是一个圆周运动在一条直线上的投影。
这两种表示方法之间并没有什么本质上的区别，就如同描绘转角大小，一圈可以用角度表示为360°，也可以用弧度表示为2π弧度一样。只是采用了两种不一样的表达形式而已，见：一圈为何是360°？。
当我们描述不论是上述的往复运动还是匀速圆周运动，必须且只需三个量即可唯一确定该运动的状态。若要描述匀速圆周运动，需要知道圆周运动轨迹的大小（即半径或幅度）；圆周运动的快慢（即角速度或频率）；以及运动的起始位置（即初始相位角），两信号起始位置之间的角度差又称为相位差。
公式表示为：
其中， 即为正弦波的幅值，
为角速度或角频率，
为初始相位角，
为时间，
为转角。
那么所有圆周运动（或振荡信号）组合起来得到的位置随时间的变化情况也就是我们最终的信号。这和从原材料得到最终的“橘子香蕉牛奶冰沙”过程类似。
同样，如果反过来，傅里叶级数能够将任何周期信号分解成一个（甚至是由无穷多个元素组成的）简单振荡信号的集合。
下面给出三种比较常见且简单的信号的分解与合成过程，这三种信号分别是方波、锯齿波和三角波。这三种信号有一个共同点是：它们都是由无数个无相位差（假设初试相位角均为零）、且成谐波关系的正弦波构成的周期性信号。
（1）方波
方波也称为矩形波，但是这种“方方正正”的信号的确可以分解为无限多个正弦信号的组合。下图展示了方波的傅里叶级数的前50项的叠加过程，如果项数继续增加，则最终趋近方波。
虽然组成方波的这些信号都是正弦信号，但是这些正弦信号之间还需要满足一定的条件。考虑组成方波的正弦信号，方波可由以下公式表示，其中 为奇数：
这里， 称为基波频率，而
等均为ω的整数倍。这些大于基波频率，且是基波频率整数倍的各次分量称为谐波。对于方波，基波的各偶数次谐波的幅值为零。这些谐波成分也就是组成方波的原材料。
这里，引入“频域”的概念，如下图：
最前面黑色的线就是所有正弦波叠加而成的总和，也就是越来越接近矩形波的那个图形。而后面依不同颜色排列而成的正弦波就是组合为矩形波的各个分量。这些正弦波按照频率从低到高从前向后排列开来，而每一个波的振幅满足上面的合成公式。每两个正弦波之间都还有一条直线，那并不是分割线，而是振幅为 0 的偶数谐波。
如果不关心相位或假设所有正弦波之间的相位差为零，按照图示方向看去，时域的方波信号就被投影到了频域。因为前面的方波信号的横轴为时间轴，而在频域，横轴为频率。这样，一组随时间变化的时域正弦信号被表示为了频域的一组离散点。频域每个离散点的横坐标代表一个谐波频率，而其纵坐标则代表该频率的谐波所对应的振动幅度。
让我们从另一个角度去看待合成的过程。前面说到，正弦波就是一个圆周运动在一条直线上的投影。而频域各个离散点也可以理解为一个始终在绕一个圆旋转的点。这些点转动的速度就对应各个谐波频率，而转动的半径就对应各个谐波的幅值。只要将这些谐波叠加起来，就能得到最终的信号。
（2）锯齿波
考虑组成锯齿波的正弦信号，锯齿波可由以下公式表示， 这里 为正整数：
下图展示了锯齿波的傅里叶级数的前50项的叠加过程，如果项数继续增加，则最终趋近锯齿波。
从圆周运动的角度看叠加过程如下图所示：
（3）三角波
对于三角波，与上面的两种类似，下图展示了三角波的傅里叶级数的前25项的叠加过程，如果项数继续增加，则最终趋近三角波。
从圆周运动的角度看叠加过程如下图所示：
总结
通过上面的例子可以看到，对于满足狄里赫利(Dirichlet)条件的周期信号，可以分解为一组成谐波关系的正弦信号，或者说该周期信号做傅里叶变换可以得到一组傅里叶级数。
对于周期信号，既然知道了其中的各个成分是成谐波关系的，那么频率成分就确定了。所以在不考虑相位差的情况下，问题关键是如何得到这些成谐波关系的正弦信号前的系数（或者说，谐波的幅值，也即是各个成分的大小）。而傅里叶变换的公式恰恰就给了我们解决该问题途径。也就是本文最开始那个公式了。由待分析的周期信号 ，可以积分得到其中所包含的谐波成分的幅值
，而将这些频率成分全部相加则可以重构出原周期信号。
为了方便理解，将傅里叶级数的求解用下式表达：
注意，这里用复指数信号来表示正弦信号，这也是两种可以互相转化的表达方法，见下图，或见：“上帝公式”真的神到无法触碰？。
傅里叶变换的局限性
我们知道傅里叶变化可以分析信号的频谱，那么为什么还要提出小波变换？答案就是“对非平稳过程，傅里叶变换有局限性”。看如下一个简单的信号：
做完FFT（快速傅里叶变换）后，可以在频谱上看到清晰的四条线，信号包含四个频率成分。一切没有问题。但是，如果是非平稳信号呢？
如上图，最上边的是频率始终不变的平稳信号。而下边两个则是频率随着时间改变的非平稳信号，它们同样包含和最上信号相同频率的四个成分。
做FFT后，我们发现这三个时域上有巨大差异的信号，频谱却非常一致。尤其是下边两个非平稳信号，我们从频域上无法区分它们，因为它们包含的四个频率的信号的成分确实是一样的，只是出现的先后顺序不同。
可见，傅里叶变换处理非平稳信号有天生缺陷。它只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。
然而平稳信号大多是人为制造出来的，自然界的大量信号几乎都是非平稳的，所以在比如生物医学信号分析等领域的papers中，基本看不到单纯傅里叶变换这样naive的方法。
上图所示的是一个正常人的事件相关电位。对于这样的非平稳信号，只知道包含哪些频率成分是不够的，我们还想知道各个成分出现的时间。知道信号频率随时间变化的情况，各个时刻的瞬时频率及其幅值——这也就是时频分析。
二、短时傅里叶变化(Short Time Fourier Transform, STFT)
针对傅里叶变换的上述缺陷，一个简单可行的解决方法就是——加窗，即“把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳，再做傅里叶变换，就知道在哪个时间点上出现了什么频率了”，这就是短时傅里叶变换的思想。
【专业说辞：短时傅里叶变换是一个用于语音信号处理的通用工具. 它定义了一个非常有用的时间和频率分布类, 其指定了任意信号随时间和频率变化的复数幅度. 实际上,计算短时傅里叶变换的过程是把一个较长的时间信号分成相同长度的更短的段, 在每个更短的段上计算傅里叶变换, 即傅里叶频谱.】
短时傅里叶变换通常的数学定义如下:
其中，
实现时, 短时傅里叶变换被计算为一系列加窗数据帧的快速傅里叶变换 (Fast Fourier Transform, FFT)，其中窗口随时间 “滑动” (slide) 或“跳跃” (hop) 。
但是，使用STFT仍然存在一个问题，那就是我们应该应用多宽的窗函数？
窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差。窗太宽，时域上又不够精细，时间分辨率低。（这个道理可以用海森堡不确定性原理来解释。类似于我们不能同时获取一个粒子的动量和位置，我们也不能同时获取信号绝对精准的时刻和频率。这也是一对不可兼得的矛盾体。我们不知道在某个瞬间哪个频率分量存在，我们知道的只能是在一个时间段内某个频带的分量存在。 所以绝对意义的瞬时频率是不存在的。）
上图对同一个信号（4个频率成分）采用不同宽度的窗做STFT，结果如右图。用窄窗，时频图在时间轴上分辨率很高，几个峰基本成矩形，而用宽窗则变成了绵延的矮山。但是频率轴上，窄窗明显不如下边两个宽窗精确。
所以窄窗口时间分辨率高、频率分辨率低，宽窗口时间分辨率低、频率分辨率高。对于时变的非稳态信号，高频适合小窗口，低频适合大窗口。然而STFT的窗口是固定的，在一次STFT中宽度不会变化，所以STFT还是无法满足非稳态信号变化的频率的需求。
三、小波变换
那么你可能会想到，让窗口大小变起来，多做几次STFT不就可以了吗？！没错，小波变换就有着这样的思路。
但事实上小波并不是这么做的，所以不能表述为“小波变换就是根据算法，加不等长的窗，对每一小部分进行傅里叶变换”，因为小波变换并没有采用窗的思想，更没有做傅里叶变换。
至于为什么不采用可变窗的STFT呢，个人认为是因为这样做冗余会太严重，STFT做不到正交化，这也是它的一大缺陷。
于是小波变换的出发点和STFT还是不同的。STFT是给信号加窗，分段做FFT；而小波直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间了！
---------------------------------------------
【解释】一起再回顾一下傅里叶变换，来弄清傅里叶变换为什么能得到信号各个频率成分吧！
傅里叶变换把无限长的三角函数作为基函数，这个基函数会伸缩、会平移（其实是两个正交基的分解）。缩得窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。某一个尺度（宽窄）下乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。那么我们就知道信号包含多少该频率的成分。
（看，上图这两种尺度能乘出一个大的值，所以信号包含较多的这两个频率成分，在频谱上这两个频率会出现两个峰）
以上，就是粗浅意义上傅里叶变换的原理。
---------------------------------------------
如前边所说，小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。
至于为什么它叫“小波”，从直观上来讲它确实看起来是个很小的波（期待更权威的解释？）
从公式可以看出，不同于傅里叶变换，变量只有频率ω，小波变换有两个变量：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函数的伸缩，平移量 τ控制小波函数的平移。尺度就对应于频率（反比），平移量 τ就对应于时间。
当伸缩、平移到这么一种重合情况时，也会相乘得到一个大的值。这时候和傅里叶变换不同的是，这不仅可以知道信号有这样频率的成分，而且知道它在时域上存在的具体位置。
而当我们在每个尺度下都平移着和信号乘过一遍后，我们就知道信号在每个位置都包含哪些频率成分。
所以说，有了小波，再也不害怕非稳定信号了，从此可以愉快地做时频分析啦（因为做傅里叶变换只能得到一个频谱，而做小波变换却可以得到一个时频谱！）
另外，小波还有一些好处：
1. 我们知道对于突变信号，傅里叶变换存在吉布斯效应，我们用无限长的三角函数怎么也拟合不好突变信号：
然而衰减的小波就不一样了：
2. 小波可以实现正交化，短时傅里叶变换不能。
附1：文中经典问题汇总
1. 关于海森堡不确定性原理
不确定性原理，或者叫测不准原理，最早出自量子力学，意为在微观世界，粒子的位置与动量不可同时被确定。但是这个原理并不局限于量子力学，有很多物理量都有这样的特征，比如能量和时间、角动量和角度。体现在信号领域就是时域和频域。不过更准确一点的表述应该是：一个信号不能在时空域和频域上同时过于集中；一个函数时域越“窄”，它经傅里叶变换的频域后就越“宽”。
如果有兴趣深入研究一下的话，这个原理其实非常耐人寻味。信号处理中的一些新理论在根本上都和它有所相连，比如压缩感知。如果你剥开它复杂的数学描述，最后会发现它在本质上能实现就源于不确定性原理。而且大家不觉得这样一些矛盾的东西在哲学意义上也很奇妙吗，世界观感觉就此被改变了。
2. 关于正交化
什么是正交化？为什么说小波能实现正交化是优势?
简单说，如果采用正交基，变换域系数会没有冗余信息，等于是用最少的数据表达最大的信息量，利于数值压缩等领域。JPEG2000压缩就是用正交小波变换。
比如典型的正交基：二维笛卡尔坐标系的（1,0）、（0,1），用它们表达一个信号显然非常高效，计算简单。而如果用三个互成120°的向量表达，则会有信息冗余，有重复表达。
但是并不意味着正交一定优于不正交。比如如果是做图像增强，有时候反而希望能有一些冗余信息，更利于对噪声的抑制和对某些特征的增强。
3. 关于瞬时频率
原问题：图中时刻点对应一频率值，一个时刻点只有一个信号值，又怎么能得到他的频率呢？
很好的问题。如文中所说，绝对意义的瞬时频率其实是不存在的。单看一个时刻点的一个信号值，当然得不到它的频率。我们只不过是用很短的一段信号的频率作为该时刻的频率，所以我们得到的只是时间分辨率有限的近似分析结果。这一想法在STFT上体现得很明显。小波等时频分析方法，如用衰减的基函数去测定信号的瞬时频率，思想也类似。
4. 关于小波变换的缺点
这要看和谁比了。
A.作为图像处理方法，和多尺度几何分析方法（超小波）比：
对于图像这种二维信号的话，二维小波变换只能沿2个方向进行，对图像中点的信息表达还可以，但是对线就比较差，这时候ridgelet（脊波）, curvelet（曲波）等多尺度几何分析方法就更有优势了。
B. 作为时频分析方法，和HHT比：
相比于HHT等时频分析方法，小波依然没脱离海森堡测不准原理的束缚，某种尺度下，不能在时间和频率上同时具有很高的精度；以及小波是非适应性的，基函数选定了就不改了。
附2：代码实现
1. 快速傅里叶变换(FFT)的python实现
# coding=utf-8 # from matplotlib.font_manager import FontProperties import pylab as pl import numpy as np # chinese_font = FontProperties(fname='XXX.ttc') sampling_rate = 8000 fft_size = 512 t = np.arange(0, 1.0, 1.0/sampling_rate) x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t) xs = x[:fft_size] xf = np.fft.rfft(xs)/fft_size freqs = np.linspace(0, sampling_rate/2, fft_size/2+1) xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100)) pl.figure(figsize=(8,4)) pl.subplot(211) pl.plot(t[:fft_size], xs, 'g') pl.xlabel(u"Time(Second)") # fontproperties=chinese_font pl.title(u"Waveform and Spectrum of 156.25Hz/234.375Hz") pl.subplot(212) pl.plot(freqs, xfp, 'b') pl.xlabel(u"Frequency(Hz)") pl.subplots_adjust(hspace=0.4) pl.show()
代码分析：
用取样时间数组t可以很方便地调用函数计算出波形数据，这里计算的是两个正弦波的叠加，一个频率是156.25Hz，一个是234.375Hz，为什么选择这两个奇怪的频率呢？因为这两个频率的正弦波在512个取样点中正好有整数个周期。满足这个条件波形的FFT结果能够精确地反映其频谱：
x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)N点FFT能精确计算的频率：
假设取样频率为fs, 取波形中的N个数据进行FFT变换。那么这N点数据包含整数个周期的波形时，FFT所计算的结果是精确的。于是能精确计算的波形的周期是: n*fs/N。对于8kHz取样，512点FFT来说，8000/512.0 = 15.625Hz，前面的156.25Hz和234.375Hz正好是其10倍和15倍。
从波形数据x中截取fft_size个点进行fft计算。np.fft库中提供了一个rfft函数，它方便我们对实数信号进行FFT计算。根据FFT计算公式，为了正确显示波形能量，还需要将rfft函数的结果除以fft_size：
xs = x[:fft_size]
xf = np.fft.rfft(xs)/fft_sizerfft函数的返回值是N/2+1个复数，分别表示从0(Hz)到sampling_rate/2(Hz)的N/2+1点频率的成分。于是可以通过下面的np.linspace计算出返回值中每个下标对应的真正的频率：
freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)最后我们计算每个频率分量的幅值，并通过 20*np.log10() 将其转换为以db单位的值。为了防止0幅值的成分造成log10无法计算，我们调用np.clip对xf的幅值进行上下限处理：
xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))剩下的程序就是将时域波形和频域波形绘制出来。
运行结果：
2. 短时傅里叶变换的python实现
import numpy as np
from matplotlib import pyplot as plt
import scipy.io.wavfile as wav
from numpy.lib import stride_tricks
""" short time fourier transform of audio signal """
def stft(sig, frameSize, overlapFac=0.5, window=np.hanning):
win = window(frameSize)
hopSize = int(frameSize - np.floor(overlapFac * frameSize))
# zeros at beginning (thus center of 1st window should be for sample nr. 0)
samples = np.append(np.zeros(int(np.floor(frameSize/2.0))), sig)
# cols for windowing
cols = np.ceil( (len(samples) - frameSize) / float(hopSize)) + 1
# zeros at end (thus samples can be fully covered by frames)
samples = np.append(samples, np.zeros(frameSize))
frames = stride_tricks.as_strided(samples, shape=(int(cols), frameSize), strides=(samples.strides[0]*hopSize, samples.strides[0])).copy()
frames *= win
return np.fft.rfft(frames)
""" scale frequency axis logarithmically """
def logscale_spec(spec, sr=44100, factor=20.):
timebins, freqbins = np.shape(spec)
scale = np.linspace(0, 1, freqbins) ** factor
scale *= (freqbins-1)/max(scale)
scale = np.unique(np.round(scale))
# create spectrogram with new freq bins
newspec = np.complex128(np.zeros([timebins, len(scale)]))
for i in range(0, len(scale)):
if i == len(scale)-1:
newspec[:,i] = np.sum(spec[:,int(scale[i]):], axis=1)
else:
newspec[:,i] = np.sum(spec[:,int(scale[i]):int(scale[i+1])], axis=1)
# list center freq of bins
allfreqs = np.abs(np.fft.fftfreq(freqbins*2, 1./sr)[:freqbins+1])
freqs = []
for i in range(0, len(scale)):
if i == len(scale)-1:
freqs += [np.mean(allfreqs[int(scale[i]):])]
else:
freqs += [np.mean(allfreqs[int(scale[i]):int(scale[i+1])])]
return newspec, freqs
""" plot spectrogram"""
def plotstft(audiopath, binsize=2**10, plotpath=None, colormap="jet"):
samplerate, samples = wav.read(audiopath)
s = stft(samples, binsize)
sshow, freq = logscale_spec(s, factor=1.0, sr=samplerate)
ims = 20.*np.log10(np.abs(sshow)/10e-6) # amplitude to decibel
timebins, freqbins = np.shape(ims)
print("timebins: ", timebins)
print("freqbins: ", freqbins)
plt.figure(figsize=(15, 7.5))
plt.imshow(np.transpose(ims), origin="lower", aspect="auto", cmap=colormap, interpolation="none")
plt.colorbar()
plt.xlabel("time (s)")
plt.ylabel("frequency (hz)")
plt.xlim([0, timebins-1])
plt.ylim([0, freqbins])
xlocs = np.float32(np.linspace(0, timebins-1, 5))
plt.xticks(xlocs, ["%.02f" % l for l in ((xlocs*len(samples)/timebins)+(0.5*binsize))/samplerate])
ylocs = np.int16(np.round(np.linspace(0, freqbins-1, 10)))
plt.yticks(ylocs, ["%.02f" % freq[i] for i in ylocs])
if plotpath:
plt.savefig(plotpath, bbox_inches="tight")
else:
plt.show()
plt.clf()
return ims
# "OSR_us_000_0011_8k.wav" is downloaded from http://www.voiptroubleshooter.com/open_speech/american.html
ims = plotstft("OSR_us_000_0011_8k.wav")
运行结果：
3. 小波变换的python实现
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import pywt
#from matplotlib.font_manager import FontProperties
#chinese_font = FontProperties(fname='XXX.ttc')
sampling_rate = 1024
t = np.arange(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
fc = pywt.central_frequency(wavename)
cparam = 2 * fc * totalscal
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(Second)") # fontproperties=chinese_font
plt.title(u"300Hz 200Hz 100Hz Time spectrum")
plt.subplot(212)
plt.contourf(t, frequencies, abs(cwtmatr))
plt.ylabel(u"Frequency(Hz)")
plt.xlabel(u"Time(Second)")
plt.subplots_adjust(hspace=0.4)
plt.show()
运行结果：
参考：
https://www.sohu.com/a/154009298_465219
https://www.zhihu.com/question/279808864/answer/552617806