草庐IT

小波变换(深入浅出)

赵孝正 2023-04-15 原文

目录

前半部分引自知乎:咚懂咚懂咚
​https://zhuanlan.zhihu.com/p/22450818

从傅里叶变换(傅里叶变换原理)到小波变换,并不是一个完全抽象的东西,可以讲得很形象。

小波变换有着明确的物理意义,下面我就按照傅里叶–>短时傅里叶变换–>小波变换的顺序,讲一下小波变换。

一、傅里叶变换

默认大家现在正处在理解了傅里叶变换,但还没理解小波的道路上。

下面我们主要讲傅里叶变换的不足。即我们知道傅里叶变换可以分析信号的频谱,那么为什么还要提出小波变换?

答案就是,“对非平稳过程,傅里叶变换有局限性”。看如下一个简单的信号:

做完FFT(快速傅里叶变换)后,可以在频谱上看到清晰的四条线,信号包含四个频率成分。

一切没有问题。但是,如果是频率随着时间变化的非平稳信号呢?

如上图,

  • 最上边的是频率始终不变的平稳信号。
  • 而下边两个则是频率随着时间改变的非平稳信号,它们同样包含和最上信号相同频率的四个成分。

做FFT后,我们发现这三个时域上有巨大差异的信号,频谱(幅值谱)却非常一致。

尤其是下边两个非平稳信号,我们从频谱上无法区分它们,因为它们包含的四个频率的信号的成分确实是一样的,只是出现的先后顺序不同。

可见,傅里叶变换处理非平稳信号有天生缺陷。它只能获取一段信号总体上包含哪些频率的成分,但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号,可能频谱图一样

然而平稳信号大多是人为制造出来的,自然界的大量信号几乎都是非平稳的,所以在比如生物医学信号分析等领域的论文中,基本看不到单纯傅里叶变换这样 naive 的方法。

一个正常人的事件相关电位M
.

对于这样的非平稳信号,只知道包含哪些频率成分是不够的,我们还想知道各个成分出现的时间

知道信号频率随时间变化的情况,各个时刻的瞬时频率及其幅值——这也就是时频分析

二、短时傅里叶变换(Short-time Fourier Transform, STFT)

为了弥补 FT 的不足,把整个时间域分解成无数个等长的小过程(加窗),每个过程近似平稳,再傅里叶变换,就知道在哪个时间点上出现了什么频率。

一个简单可行的方法就是——加窗

“把整个时域过程分解成无数个等长的小过程,每个小过程近似平稳,再傅里叶变换,就知道在哪个时间点上出现了什么频率了。”

这就是短时傅里叶变换

看图:

时域上分成一段一段做FFT,不就知道频率成分随着时间的变化情况了吗!

用这样的方法,可以得到一个信号的时频图了:


——此图像来源于“THE WAVELET TUTORIAL”

图上既能看到10Hz, 25 Hz, 50 Hz, 100 Hz四个频域成分,还能看到出现的时间。(两排峰是对称的,所以大家只用看一排就行了。)

是不是棒棒的?时频分析结果到手。但是STFT依然有缺陷

使用 STFT 存在一个问题,我们应该用多宽的窗函数?

窗太宽太窄都有问题:


如上图所示:

  • 窗太窄,窗内的信号太短,会导致频率分析不够精准,频率分辨率差。
  • 窗太宽,时域上又不够精细,时间分辨率低。
  • (这里插一句,这个道理可以用海森堡不确定性原理来解释。
  • 类似于我们不能同时获取一个粒子的动量和位置,我们也不能同时获取信号绝对精准的时刻和频率。
  • 这也是一对不可兼得的矛盾体。
  • 我们不知道在某个瞬间哪个频率分量存在,我们知道的只能是在一个时间段内某个频带的分量存在。
  • 所以绝对意义的瞬时频率是不存在的。)

为了弥补FT的不足,把整个时间域分解成无数个等长的小过程(加窗),每个过程近似平稳,再傅里叶变换,就知道在哪个时间点上出现了什么频率。

但此时仍然有问题,那就是窗的宽度。窄窗口时间分辨率高、频率分辨率低;宽窗口时间分辨率低、频率分辨率高(如下图)。所以,对于时变非稳态信号,高频部分适合用小窗(短周期),低频部分适合用大窗(长周期)。然而,在一次 STFT 中,窗口的宽度是固定。所以 STFT 也有其局限性,这就引出了我们的小波变换。

  • 下图左(Narrow Window)中,时间分辨率高,可以理解成:
    • 图中Time轴上,每块波对应一个时间窗口,每两个不同频率块(波)在时间上没有交叉,时间上都是连续且不重叠的。
    • Frequency轴方向上,每个频率分辨率低,即不能看出具体频率是多少,沿Frequency投影的话,会得到一个开口向下的抛物线。
  • 下图右(Wide Window)中,
    • 频率分辨率高,可以理解为从频率Frequency轴投影,每个频率都是没有交叉的,都是一条直线,
    • 而从Time轴看的话,时间上都是交叉的,所以时间分辨率低,宽窗口会导致丢失一些小的事件。

看看实例效果吧:

a 取不同的值,代表了不同的窗口宽度,值越小代表窗口越宽。

上图对同一个信号(4个频率成分)采用不同宽度的窗做STFT,结果如右图。用窄窗,时频图在时间轴上分辨率很高,几个峰基本成矩形,而用宽窗则变成了绵延的矮山。但是频率轴上,窄窗明显不如下边两个宽窗精确。

所以窄窗口时间分辨率高、频率分辨率低,宽窗口时间分辨率低、频率分辨率高。对于时变的非稳态信号,高频适合小窗口,低频适合大窗口。然而STFT的窗口是固定的,在一次STFT中宽度不会变化,所以STFT还是无法满足非稳态信号变化的频率的需求。

但此时仍然有问题,那就是窗的宽度。窄窗口时间分辨率高、频率分辨率低;宽窗口时间分辨率低、频率分辨率高(如下图)。所以,对于时变非稳态信号,高频部分适合用小窗,低频部分适合用大窗。然而,在一次STFT中,窗口的宽度是固定。所以STFT也有其局限性,这就引出了我们的小波变换。

三、小波变换

那么你可能会想到,让窗口大小变起来,多做几次 STFT(短时傅里叶变) 不就可以了吗?!没错,小波变换就有着这样的思路。

但事实上小波并不是这么做的

  • (关于这一点,认为 “小波变换就是根据算法,加不等长的窗,对每一小部分进行傅里叶变换”,就不准确了。
  • 小波变换并没有采用窗的思想,更没有做傅里叶变换。)

至于为什么不采用可变窗的 STFT (短时傅里叶变) 呢,我认为是因为这样做冗余会太严重,STFT 做不到正交化,这也是它的一大缺陷。

于是小波变换的出发点和STFT还是不同的。

  • STFT是给信号加窗,分段做FFT
  • 而小波直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基
  • 这样不仅能够获取频率,还可以定位到时间了~

【解释】

来我们再回顾一下傅里叶变换吧,没弄清傅里叶变换为什么能得到信号各个频率成分的同学也可以再借我的图理解一下。

傅里叶变换把无限长的三角函数作为基函数:


这个基函数会伸缩、会平移(其实本质并非平移,而是两个正交基的分解)。

缩得窄,对应高频;伸得宽,对应低频。然后这个基函数不断和信号做相乘。

某一个尺度(宽窄)下乘出来的结果,就可以理解成信号所包含的当前尺度对应频率成分有多少。

于是,基函数会在某些尺度(宽窄)下,与信号相乘得到一个很大的值,因为此时二者有一种重合关系。那么我们就知道信号包含该频率的成分的多少。

仔细体会可以发现,这一步其实是在计算信号和三角函数的相关性


看,这两种尺度能乘出一个大的值(相关度高),所以信号包含较多的这两个频率成分,在频谱上这两个频率会出现两个峰。

以上,就是粗浅意义上傅里叶变换的原理。

如前边所说,小波做的改变就在于,将无限长的三角函数基换成了有限长的会衰减的小波基。

这就是为什么它叫“小波”,因为是很小的一个波嘛~

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

当伸缩、平移到这么一种重合情况时,也会相乘得到一个大的值。这时候和傅里叶变换不同的是,这不仅可以知道信号有这样频率的成分,而且知道它在时域上存在的具体位置。

而当我们在每个尺度下都平移着和信号乘过一遍后,我们就知道信号在每个位置都包含哪些频率成分

看到了吗?有了小波,我们从此再也不害怕非稳定信号啦!从此可以做时频分析啦!

做傅里叶变换只能得到一个频谱,做小波变换却可以得到一个时频谱

↑:时域信号

↑:傅里叶变换结果


——此图像来源于“THE WAVELET TUTORIAL”

↑:小波变换结果

小波还有一些好处,比如,我们知道对于突变信号,傅里叶变换存在吉布斯效应,我们用无限长的三角函数怎么也拟合不好突变信号:

然而衰减的小波就不一样了:

但是真正理解透小波变换,这些还差得很远。

比如你至少还要知道有一个“尺度函数”的存在,它是构造“小波函数”的关键,并且是它和小波函数一起才构成了小波多分辨率分析,理解了它才有可能利用小波做一些数字信号处理;你还要理解离散小波变换、正交小波变换、二维小波变换、小波包……这些内容国内教材上讲得也很糟糕,大家就一点一点啃吧~

四、小波变换的详解

4.1、小波变换的基

小波变换干脆直接把傅里叶变换的基给换了,将无限长的三角函数基换成了有限长的会衰减的小波基,它的能量有限,都集中在某一点附近,而且积分的值为零。

傅里叶变换,变量只有 w \pmb{w} w,而小波变换则有尺度 a \pmb{a} a 和平移量 b \pmb{b} b尺度对应于频率平移量对应于时间,所以小波变换可以用于时频分析,得到信号的时频谱。

下面是小波函数的一般形式。

4.3、CWT(连续小波变换)

变换公式如下,通过参数 a \pmb{a} a 定位频率,通过参数 b \pmb{b} b 定位时间

变换后可以得到如下的标度图,两个坐标轴分别为 timeScales(即 b \pmb{b} b a \pmb{a} aLower Scales 对应高频成分,High Scales 对应低频成分。

4.3.2、小波变换的步骤(CWT)

(1)把小波 w ( t ) w(t) w(t) 和原函数 f ( t ) f(t) f(t) 的开始部分进行比较(实际上就是作內积),计算系数 C \pmb{C} C。系数 C \pmb{C} C 表示该部分函数与小波的相似程度。
(2)把小波向右移 k \pmb{k} k 单位,得到小波 w ( t − k ) w(t-k) w(tk),重复 1。重复该步骤直至函数 f \pmb{f} f 结束.
(3)扩展小波 w ( t ) w(t) w(t),得到小波 w ( t / 2 ) w(t/2) w(t/2),重复步骤 1, 2.
(4)不断扩展小波,重复1, 2, 3

首行实际上是在 scale1 上所有的系数,它们是C11-C1n
在这里,所有的系数在 scale2 上排列,它们是C21-C2n
在这里,所有的系数在 scale3 上排列,它们是C31-C3n

那么,在 scale1 上的系数实际上是如何计算的?
你看上图左侧 Wavelet at scale1 ,我们选择这样一个小波,它被压缩了,所以,有一个高频,在此情况中,这个小波会捕获信号中所有可用的高频分量,所以,我们使用 a 压缩这个小波,使用 b沿着时间轴平移小波,然后与这个信号相乘,然后积分。

Wavelet at scale1

那么,我们在这里得到所有的系数进行排列,如上图所示。

对于 scale2 ,我们选择下图这样一个小波,它比 scale1 有点拉伸,所以,它的频率会降低,它会捕获低频细节。然后它再次沿着时间轴滑动并和信号相乘并积分。那么,我们得到了所有这些分量。

Wavelet at scale2

scale3 上,它会进一步拉伸,并沿着时间轴滑动,逐步得到 scale3 的所有系数,这就是小波系数的排列方式

Wavelet at scale3

举例介绍:
首先说明这里是连续小波变换,不会涉及离散小波变换,不涉及尺度函数。

对于一个 morlet 小波变换,其频率成分分别是 4610,具体表达式如下:

经小波变换,得到如下时频图:

图1、二维时频图
·

图2、三维时频图

其实,图1是图2的俯视图(向下投影)

从图中可以清晰辨认出,原始信号含有的频率成分及其各自对应的时间区间

如果你和我一样认真的话,你会产生三个疑问?

  • (1)为什么小波变换能确定信号频率和其对应的时间区间?
  • (2)为什么频率小的条纹比频率高的条纹要细?
  • (3)为什么三个条纹颜色的深浅不一样 ?

其实第(1)个问题可以从前面的内容中总结出来。

第一个问题我们把它一分为二,为什么 CWT 能辨认出信号的频率成分?为什么 CWT 能确定频率对应的时间区间?

当然先看下公式。来一个 morlet 小波基函数表达式

morlet 小波的基函数是由复三角函数乘上一个指数衰减函数构成,这里 ω 0 \omega_0 ω0 表示中心频率。还是先看图(分别令 ω 0 = 2 \omega_0=2 ω0=2 ω 0 = 5 \omega_0=5 ω0=5

图三

上图左边是基函数时域图像,右边是傅里叶变换图像。可以看到基函数的频率正是其中心频率的值。这里给出morlet小波的傅里叶变换表达式

从表达式也可以看出当频率等于其中心频率 ω 0 \omega_0 ω0 时,取极大值。

这里复三角函数可以辨认频率,衰减函数可以保证其时域的有限支撑。

只给一个固定的中心频率可不能辨认信号的频率,同样,基函数只在【-22】之间也确定不了时间区间。

所以这里的小波基函数需要平移和伸缩 。又是一个公式(按传统,这个变大时前还有一个根号 a \pmb{a} a,但我不讨论这个)

参考资料

[1] 小波变换的引入,通俗易懂 2018.7.25;
[2] 连续小波变换(CWT)2018.7.25;

有关小波变换(深入浅出)的更多相关文章

  1. ChatGPT教程之深入了解魔术背后的技术 - 2

    解开谜团:深入探索ChatGPT的技术奇迹。ChatGpt无处不在,无论是在播客、博客、YouTube还是社交媒体上。当我注意到这项新技术如此受欢迎时,我决定试一试,我被震惊了!有很多关于ChatGpt及其魔力的博客,但在这篇博客中,我将深入探讨其内部技术及其工作原理!ChatGpt简介根据OpenAI,ChatGpt被描述为:“我们训练了一个名为ChatGpt的模型,它以对话方式进行交互。对话格式使ChatGpt可以回答后续问题、承认错误、挑战不正确的前提并拒绝不适当的请求。ChatGPT是InstructGPT的兄弟模型,它经过训练可以按照提示中的说明进行操作并提供详细的响应。”OpenA

  2. ruby - 从文件中提取快速傅里叶变换数据 - 2

    我正在构建一个应该在服务器上运行并分析声音文件的工具。我想在Ruby中执行此操作,因为我的所有其他工具也是用Ruby编写的。但我很难找到完成此任务的好方法。我发现的很多例子都是在做可视化和图形化的东西。我只需要FFT数据,仅此而已。我既需要获取音频数据,又需要对其进行FFT。我的最终目标是计算一些东西,例如所有频率(加权幅度)的均值/中值/众数、第25个百分位数和第75个百分位数、BPM,也许还有其他一些好的特性,以便以后能够将相似的声音聚集在一起.首先,我尝试使用ruby-audio和fftw3,但我从未将两者真正结合使用。文档也不好,所以我真的不知道有什么数据被洗牌了。接下来,我尝

  3. 图形学-变换(平移矩阵,旋转矩阵,缩放矩阵,线性变换,仿射变换,齐次坐标) - 2

    1.变换1.1什么是变换?变换(Transform)是计算机图形学中非常重要的一部分。变换包含模型变换(Modelingtransform)以及视图变换(Viewtransform)。模型变换指的是变换模型(被拍摄物体)的位置,大小和角度;视图变换指的是变换照相机的位置和角度。从相对运动的角度来看,两种变换是可以相互转化的。1.2模型变换1.2.1二维变换缩放变换缩放变换(Scale)中,如果一个图片以原点(0,0)为中心缩放𝑠倍。那么点(𝑥,𝑦)变换后数学形式可以表示为写成矩阵形式为:当然,我们也可以给x轴和y轴不同的缩放倍数𝑠𝑥和𝑠𝑦。在非均匀情况下,缩放变换的矩阵形式为反射变换反射变换(

  4. 科大讯飞刘聪:由ChatGPT浪潮引发的深入思考与落地展望 - 2

    近期,以“生成式人工智能”(GenerativeAI)为核心技术的聊天机器人ChatGPT火爆全球。百度、阿里巴巴、科大讯飞、360等国内企业纷纷抛出ChatGPT相关进展,打造中国版的ChatGPT。科大讯飞此前在投资者互动平台表示,ChatGPT主要涉及到自然语言处理相关技术,属于认知智能领域的应用之一,公司在该方向技术和应用具备长期深厚的积累。并称2022年12月已进一步启动生成式预训练大模型任务攻关,类ChatGPT技术将在今年5月率先落地科大讯飞AI学习机产品。近日,科大讯飞副总裁、研究院执行院长刘聪围绕什么是ChatGPT,它强在哪里?会对未来世界带来哪些颠覆性影响?进一步阐述Ch

  5. 深入理解C++中的move和forward! - 2

    导语 |  在C++11标准之前,C++中默认的传值类型均为Copy语义,即:不论是指针类型还是值类型,都将会在进行函数调用时被完整的复制一份!对于非指针而言,开销及其巨大!因此在C++11以后,引入了右值和Move语义,极大地提高了效率。本文介绍了在此场景下两个常用的标准库函数:move和forward。一、特性背景(一)Copy语义简述C++中默认为Copy语义,因此存在大量开销。以下面的代码为例:0_copy_semantics.cc#include#includeclassObject{public:Object(){std::coutv;v.push_back(obj);}最终的输出

  6. 深入理解Linux文件系统与日志分析 - 2

    目录引言:一、inode和block1、inode和block概述2、inode的内容1.inode包含文件的元信息(文件属性)2.用stat命令可以查看某个文件的inode信息3.Linux系统文件三个主要的时间属性  4.目录文件的结构3、inode的号码​5、硬盘分区后的结构6、inode的大小7、inode的特殊作用 二、链接文件三、案例:恢复EXT类型的文件四、案例:恢复XFS类型的文件五、日志文件1.日志的功能2.日志文件的分类3.日志保存位置1.常见的一些日志文件:2.扩展:日志检查3.小结:​4.日志消息的级别5.用户日志分析六、总结引言:inode是一个重要概念,是理解Uni

  7. javascript - CSS 3D变换以制作给定边长的梯形 - 2

    我有一个给定尺寸(比如100x300像素)的元素,它位于高度相同且宽度可变的容器中,我想使用rotateX围绕-webkit-transform-进行转换origin:topcenter;在选择容器的-webkit-perspective时,图像的底线看起来保持在原处,但只会扩展以填充整个容器。哇,这听起来令人困惑。这是一张照片:基本上,我想创建一个上部宽度固定、下部宽度可变的梯形。但是我不能完全弄清楚关系背后的数学......欢迎使用Javascript。以下示例适用于正文宽度为600像素的情况:http://jsfiddle.net/24qrQ/现在的任务是随着body的宽度不断改

  8. javascript - 对于 JavaScript 多维数组的深拷贝,深入一层似乎就足够了。这是真的吗? - 2

    注意:我只是一个编码新手,所以这个问题的核心可能存在明显的错误或误解。本质上,我需要在JavaScript中“按值”深度复制多维数组到未知深度。我原以为这需要一些复杂的递归,但似乎在JavaScript中您只需要深复制一个级别就可以按值复制整个数组。举个例子,这是我的测试代码,使用了一个故意复杂的数组。functiontest(){vararr=[['ok1'],[],[[],[],[[],[[['ok2'],[]]]]]];varcloned=cloneArray(arr);arr='';//Deletetheoriginalalert(cloned);}functioncloneA

  9. 坐标变换与坐标系变换 - 2

            在slam中,经常充斥着坐标变换,坐标系变换等等知识,经常把我们搞得头大。这篇文章稍微记录一些我个人的一些心得体会,给同是困惑的人一些参考。       在我看来,坐标系变换与坐标变换在虽然是听起来差不多的两个名词,但其实表达的是两个基本相反的内容。比如坐标系A到坐标系B的坐标系变换主要还是强调坐标系A如何通过这个变换矩阵变到坐标系B,在这里,转动的是坐标系A的相关坐标轴,就好像字面意思那样,起始位置在A(就好像把A当做参考系一样)。坐标变换在slam中其实是针对某些具体点来讲的(比如某些静止不动的地图点),我们需要对这些静止不动的点做一些坐标转换,比如将点的坐标从激光雷达坐标

  10. javascript - 变换矩阵线性组合的旋转动画导致放大缩小 - 2

    我有一个3x3矩阵(startMatrix),它表示图像的实际View(平移、旋转和缩放)。现在我创建一个新矩阵(endMatrix),它有一个恒等矩阵、新的x和y坐标、新的Angular和新的比例,例如:endMatrix=translate(identityMatrix,-x,-y);endMatrix=rotate(endMatrix,angle);endMatrix=scale(endMatrix,scale);endMatrix=translate(endMatrix,(screen.width/2)/scale,screen.height/2)/scale);和功能(标准的

随机推荐