图形数学方法:FFT(快速傅里叶)的原理讲述 前言: 本系列文章算实时渲染卷轴系列的前置文章,但论知识含量和知识价值并不比正式内容低,可以说是两个毫不相关的部分。笔者会按照基础内容,图形数学方法,实时渲染前置,实时渲染卷轴,对专栏进行划分。但其中的内容会互相夹杂,主要目的是保证一篇文章的阅读流畅性,本文中会讲述FFT算法。并不会讲述图形运用,其图形应用我会放在实时渲染前置当中讲述。 基础原理讲述: FFT(快速傅里叶变换): FFT算法是DFT算法的改良版,而DFT是FFT的离散化。理解FFT,就从傅里叶变换到DFT再到FFT的思路进行推导。笔者也会按照这样的思路进行讲解推导。 傅里叶变换: 傅里叶变换是傅里叶级数的推广,所以在谈傅里叶变换之间,先说一下傅里叶级数。在大学期间学习无穷级数有相关基础的同学可以跳着看。 傅里叶级数: 傅里叶级数是把类似波的函数表示成简单正弦波的方式,更严肃来说的话:对于满足狄利克雷定理的周期函数,其傅里叶级数是由一组简单的振荡函数加权和表示的,表示周期函数为正弦波和余弦波之和。和或谐波(谐波频率是原周期信号频率整数倍的波),可以用谐波分析开确定每一个谐波的相位和幅度。傅里叶级数中就可能有无限谐波数。对于函数的傅里叶级数的部分但不是所有的谐波求和会产生该函数的近似值,例如:傅里叶级数前几个谐波用于方波就会产生方波的近似值。
方波(表示为蓝点)近似为其第六部分和(表示为紫点),由方波傅里叶级数的前六项(表示为箭头)求和形成。每个箭头从其左侧所有箭头的垂直总和开始(即先前的部分总和)
方波的傅立叶级数的前四个部分和。随着更多谐波的添加,部分和会收敛到(变得越来越像)方波 相关历史: 傅里叶级数来源于法国的数学家约瑟夫-傅里叶,他认为任何函数都可以展开为三角级数。在这之前的数学家,例如拉格朗日已经找到了一些非周期函数的三角级数展开,在认定了一个函数有三角级数的展开之后,通过积分的方法计算其系数的公式,则是被欧拉等人发现的。傅里叶通过使用三角级数来解热传导方程(热传导方程(或称热方程)是一个重要的偏微分方程,它描述一个区域内的温度如何随时间变化),随后便得到了发展。但实际上周期函数可以分解为简单振荡函数的总和的想法来自于公前3世纪天文学家的均轮和本轮的学说。而现在傅里叶级数在数论,组合数学,信号处理,概率论,统计学,密码学,声学,光学有广泛的应用。 傅里叶级数推导: 三角正交性质: 三角正交的性质在傅里叶级数中是类似于基础一样的性质,后续的推导过程都需要这一步作为基础。提到三角形正交性质,首先就得说三角函数集,集也就是集合。例如以下的一个集合:{0,1sinx,cosx,sin2x,cos2x……sinnx,cosnx} 这就是一个三角函数的集合,包括0和1也就是sin0x,con0x。很好理解,那么如何理解三角正交呢?正交也就是垂直的意思,首先来看向量的垂直: a x b =0(a,b都是向量) 也就是|a|x|b|xcos θ=0。我们把这样的向量称之为垂直,也能称之为正交。那么代数表达呢? a=(-1,2),b=(2,1)的话,其内积就是:-12+21=0。 哪怕再扩展一个维度也是这样计算。既然这样的结果表达我们称之为正交性,那么三角函数的正交性是怎么来的? 任何一个三角函数能不能拿出自己的任意一个点,再与另一个三角函数的点相乘。并把这样的行为一直重复在整个定义域,结果相加,得到0。这是可行的,实际上这也正是三角函数正交性的本质行为。而数学公式表达如下:
其中的满足要求则是m≠n 这只是三角函数正交性的一个数学公式,无论m与n的值取多少,只要两者不想等就一定满足上面的公式。让我们从几何的层面来理解:
这是sinx的函数图像,上面曾说正交是垂直,然后得出了向量的内积计算,请联想内积计算和sinx。假设计算sin2x和cos3x的积,在其正交性的表达式中,得到的值是0,但如何联系正交,我们把sin2x和cos3x的函数图像拆开成单点来看:
假设取x=1的点,在sin2x和cos3x中相乘,并把这个行为重复在定义域,得到的结果就是0,这个行为类似于前面向量内积的计算过程。关于正交的理解从这里入手较为方便,说完了如何理解,接下来就是去证明正交性,现在再让我们看公式:
现在来看这个公式,就能体会其物理含义:把sinnx和cosmx的每一个点作乘然后相加,根据结果为0,所以说三角函数有正交性。
根据积化和差公式就能得到:
进一步化简:
最后计算等于0. 而我们也证明了上面给出的三角函数正角的式子是正确的。同理可以得到以下的五个公式: 当n≠m时候:
当n=m的时候: 根据上面的同理计算可以得到:
而这就是三角函数的正交性质,有了正交性质的基本了解,接下来才能推导傅里叶级数。 周期为2π函数展开为傅里叶级数: 首先来拿出读者在逛各大技术博客的时候,经常遇到的傅里叶级数展开式子:
看到这个公式,还是稍微觉得有点亲切,但是遇到复数形式的,大部分人就不知道在干什么了。而这个说明其实还是没有吃透傅里叶级数的问题.继续说傅里叶级数,傅里叶认为周期函数都可以用三角函数的加权和来表示,用数学语言来表达,也就是:
这个公式很直观的可以理解,但是和上面公式明显不一样,仔细观察是n从0开始取值,而上面公式实际上是把n=0的值单独提取出来了,也就是如下的形式:
又有读者说了,这里为什么是a0,和上面不一样。这一点,让我们来计算a0看看,得到a0到底是什么东西再说。 对两边进行积分,得:
而后面两项由于三角函数的正交性,所以值为0,那么式子为:
继而可以得到a0的值为:
发现了吗?我们这里的a0如果在前文是二分之一的形式,那么这里的a0计算公式则会在分母的位置少一个2。这样a0的计算公式会更加好看,换言之,一切都是为了公式的简洁来的。让我们再看前面的傅里叶展开式子:
在明确了a0的计算之后,来看看an和bn的计算,首先来看an。 对于等式两边乘以cosmx再做积分:
再由于前面得到的正交性,可以知道第一项为0,第三项为0,只剩下了第二项:
然后让我们回忆一下三角函数的正交性,唯一能够保留下来的,也就只有n=m的时候,其余条件均为0。那么式子就变为:
那么就能得到an的计算公式为:
所以发现了吗?为了利用三角函数的正交性,所以拿入了cosmx和积分,来消除其他项,那么bn的计算也同理,两边拿入sinx再对式子两边积分,和上面一样的计算过程,得到bn的计算公式为:
到这里为止,关于傅里叶级数的周期函数展开的形式各项都推导一边,让我们回顾一下上面的傅里叶级数:
周期为2L的函数展开为傅里叶级数: 不可能只把2π的函数展开为傅里叶级数,因为现实世界并没有这么特殊,如何得到一般周期形式的傅里叶展开呢?换二字呼之欲出
那么原来是2L的周期,现在则变为了2T,也就能再一次运用上面的公式,只是需要稍微对公式变形处理一下:
推导过程,就不进行撰写了,简单的代换而已,推荐读者自己动手试一下。值得一说地是,在工程当中,我们很少去说周期的概念,而是去说频率ω。而且频率ω的引入,可以帮助限制时间t的范围,因为在实际生活当中,时间是不存在负数这个概念的。
而一个周期T则是2L,所以表达式又为:
很多人知道ω的表达式,但是不知道其实际出处来自于这里。所以一个完整的周期为2L的函数傅里叶展开式则为:
而这就是周期为2L的傅里叶展开式,让函数的傅里叶展开进一步的扩大范围。而我们接下来要写的式子则是傅里叶级数难以看懂的形式。 傅里叶级数的复数形式: 在谈傅里叶级数的复数形式之前,我想谈论一个欧拉方程,在后续的推导也会用到。傅里叶级数的复数形式是推导傅里叶变换(FT)的基础,请认真观看。 欧拉公式: 欧拉公式被称之为最美的公式,其在数学,物理学,工程学中无处不在,物理学家理查德费曼称欧拉公式是“数学中最杰出的公式” 相关历史: 1714年的时候,英国数学家Roger Cotes提出了一个可以解释的几何论证:
对该公式求幂也就会得到欧拉公式。 证明欧拉公式: 笔者能力有限,也就不说欧拉公式是怎么来的了,只是来证明一下欧拉公式的正确性。帮助各位读者记忆欧拉公式。
这就是欧拉公式,十分的简洁和美丽。证明欧拉公式,要来假设一个函数:
对式子两边同时求导,右边是复合函数求导,求导法则不再写了,求导结果如下:
然后化简会发现,方程右边值为0. 也就是此函数的导数为0。函数的导数为0,则说明了一个函数的导数为常数,也就是 θ为任何值的时候,都是一个常数。将 θ代为0的结果应当也是一样的,就能得到该函数的导数为1,自然原函数也就是1。那么也就能证明此函数的比值为1。就能证明得到:
理解欧拉公式: 对于计算机图形学领域来说,较为重要的是,理解其物理含义。欧拉公式建立了三角函数和复函数之间的关系。 复数: 复数是实数的延伸,最开始引入复数是因为在实数范围之类,解决不了下面的方程:
很明显在实数范围之类无解,所以为了使该方程有解,创造了在实数集之外的“假想的数”i。
便可以使得上面方程有解。也因而产生了复数 复数的数学形式为:
复平面: 每一个复数都有各自的实部和虚部,单独拿出来,参考二维平面的方式,来构造一个平面来描述复数,这就是复平面。既然有了复平面了,那么是否能构成向量,得到复向量呢?答案是可以的。
再回谈到欧拉公式的理解上。 直观理解欧拉公式: 既然复数是a+bi的形式,那肯定也能得到a(x)+b(x)i的形式。实际上在上面的欧拉方程当中,这一点就是直观的体现,因为左边的指数是复数的形式,那自然肯定可以写成a(x)+b(x)i的形式,也就是上面方程当中的:cos(x)+sin(x)*i。所以欧拉公式是什么?是一种联系,是一种转换,可以把三角函数和复指数函数相联系。把难以计算的负指数函数转化为对应的三角函数 傅里叶级数的复数形式: 回到傅里叶级数的复数形式了,这一点十分重要,也是傅里叶变换(FT)的基础,首先拿出我们推导的结果傅里叶级数的一般形式:
然后正由于有sin θ和cos θ,就能和欧拉公式建立联系,进而得到傅里叶级数的复数形式。
也能得到欧拉公式下的cos和sin的表达:
然后代入上面的公式当中:
化简进而得到:
然后到这里仔细观察式子,前面的第一项相当于第零项,后面则是两个1到正无穷大,那如果把其中一个变为负无穷到1,就能构成一个完整的负无穷到正无穷的过程了。为什么要这样构造,在于目前的傅里叶级数的复数展开式子,没有负数的对象,所以尽量要扩展其形式。
进而化简得到:
而这就是傅里叶级数的复数形式,至于网上看到的分段式表达,就是因为cn的不同。如下:
然后把上面的an和bn计算公式代入上面的公式当中,就能得到复数形式下的,Cn完整表达:
0,C_{n}=\frac{1}{T}\int_{0}^{T}f(t)e^{-inωt}dt” eeimg=”1″>
综上所述傅里叶级数的复数形式为:
傅里叶变换(FT): 傅里叶级数的复数形式如下:
公式里面,ω实则表达的是频率的概念,称之为基频率表。而实际决定函数近似谁的,是Cn。而Cn由前面可知,自然就是a+bi的一个形式,然后问题来了,能否通过一个几何的形式来表达出上面的公式:
如图所示,把傅里叶级数的决定条件:实部a,虚部b,频率ω作为了3d坐标系下面的坐标轴,就能把任意一个傅里叶级数的指数形式分解到图上。
从这里来看,就能直观地看到傅里叶级数是如何拆分到3d坐标系上面的了。 在这个图当中是根据频率作为横坐标的,所以称之为频域表达,而f(x)实际是作为傅里叶级数对f(x)的逼近,那么这个f(x)作为一个周期函数,自然是可以通过时间的形式表达出来,而且是周期t的形式
这里只是拿出了f(x)的形式,没有写周期形式。但很毋庸置疑的是,周期形式是存在的,而且是上面频率图的对应,那我们这个用周期来描述的图呢?则被称之为时域表达。而其余技术文章中经常看到的时域和频域。我们就在这里得到了一个直观的认识。 话说回来,我们一直聊的是周期函数,这太特殊了。实际工程当中,很少有周期函数,我们希望傅里叶级数有更加宽广的应用,所以如何用傅里叶级数来近似非周期函数呢?想想非周期函数和周期函数之间的不同,周期函数会重复,而非周期函数不会重复。我们只能用傅里叶级数来描述周期函数,那么能不能把非周期函数来说成周期函数了,就得运用一种数学思想:非周期函数会在无限久远之后发生重复。我们把非周期函数看待成一个周期函数,但是这个周期则是无限长,那么现在就能用傅里叶级数来近似非周期函数了。 这样看待非周期函数,会发生什么问题呢?周期无限大,那么ω就会无限小,在频率表达图中,两个频率表达之间的距离就会越来越小,近乎为0,但是之前的频域表达图,很明显可以看出来傅里叶级数的复平面表达是离散的,但由于把非周期函数看为周期函数,频域表达之间的距离近似为0,函数的表达就从离散可以近似看为连续。也就是说在这样看待非周期函数之后,离散变为了连续。请仔细地体会这个表达改变的过程,而这是后面推导的基础。 有了这些思想之后,就能正式开始推导傅里叶变换了。
拿出上面推导的公式,根据上面的思想进行推导变换:
然后由于从离散变为了连续,那么下面公式转换自然就能等价:
同时由于t无穷大:
读者如果纠结这里的0到T,怎么变成了负无穷到了正无穷了,负数部分是怎么来的。请再仔细会去看傅里叶级数的Cn推导部分。只要满足一个周期就行,所以笔者选择的0到T,也有人选择的-T/2到T/2。从后者来说,自然就有负数部分了。其实是一样的。 那么由于以上改变。傅里叶级数的复数形式变为了:
这里的nω都变为了ω,所以没有了n,这就是从离散到连续的改变的标志。 实际上,我们已经推导出来了,傅里叶变换公式:
而剩余部分则被称之为傅里叶变换的逆变换公式:
在进行了这么多的数学推导之后,我们终于得到了FT(傅里叶变换)的公式,而FFT(快速傅里叶变换)和FT(傅里叶变换)之间又是怎样的关系呢? 离散傅里叶变换(DFT): 前面曾经说过FFT是DFT的改良版本,所以先来说DFT(离散傅里叶变换)。上面得到的傅里叶变换公式是连续积分的过程,但是要用到计算机当中自然是不行的,开销太大,傅里叶变换是对时间,连续频率的无穷积分,对于计算机(作为一个离散的系统)来说自然就不是很好计算,所以要进行离散化。也就是DFT。 其实在DFT之前,为了解决FT在算法层面的问题,一开始提出的是离散时间不离散频率的傅里叶变换,也就是DTFT。而这里就不说DTFT了,因为后文用的是FFT的方法,只说DFT(离散时间离散频率的傅里叶变换)。 傅里叶变换FT公式:
傅里叶变换逆变换的公式:
笔者有意识到这部分画的时间过多了,但是值得说的是,在完成这一部分讲解时候,笔者会给出一个基于FFT算法的实时海洋渲染的Unity项目,让各位掌握这种所学马上就能所用,算法到可视化的愉快过程,这也是很多人对图形学如此上心的原因之一,包括我在内。 让我们回到DFT,DFT的数学表达式为:
这里的k的表达为:
从FT到DFT的演变,先暂且把时间离散化:
这里的离散,实际上就是在连续的过程中寻找单一的采样点,然后加起来,纵然我们知道信号或者声音是连续的,但是设备有限,我们只能进行离散化的采集。因此成了这样的式子,用数学表达:
也就是说,进行了时间上的离散化过程后变为了:
再说回ω,当时间上进行了离散化之后,频率上也应该进行离散化,根据采样的次数来界定的
也就是说现在的ω在原有的基础上更受到了采样次数N的限制,因此为:
最终我们也就获得了我们的DFT计算公式:
为什么DFT更适合计算机,因为可以用矩阵的方式进行表达,只要我们界定:
那么在矩阵当中就能用W^{2N}来对DFT计算公式进行表达,假设此时F(n)代表了采样的离散的信号,与之对应的连续信号就是f(n),那么DFT就是架在两者之间的一座桥:
说到这里的时候,读者或许云里雾里,让我们来看一个例子。现在有一个信号函数为:
进行4次采样,也就是说N=4,然后连续的f(k)函数就是
这是因为t=k/4,然后代入上面的DFT变换:
最终就能得到四次采样的离散化的图: 如果这里
如果这里听到不是很明白。欢迎来看这个视频: https://www.youtube.com/watch?v=nl9TZanwbBk FFT(快速傅里叶变换): 我不得不快点讲完这个内容,对于FFT,实则就是在计算DFT,其诞生还是由于DFT计算太慢,FFT被人为是“20世纪十大算法之一”,FFT的核心思路在于分治思想,主要介绍一下FFT的核心思路便开始我们的unity项目实践。 如果使用DFT进行慢慢计算的话,有n个数据,就需要进行n次乘法,其复杂度就会为n^2,但是如果用FFT的话,其复杂度就为n^(log_(n))。对数学敏感的同学,应该可以感知到这是多大的一个扩展,真正地帮助了DFT在计算机中的落实。FFT的分类基本就是时间抽取法和频率抽取法两类 介绍一下radix-2 DIT(时间抽取),这是FFT算法的最简单和最常见的形式,之所以叫做Radix-2,则是因为每一个递归阶段都把大小为N的DFT划分成了大小为N/2的交错DFT。 Radix-2 DIT的做法是首先计算偶数索引的DFT再计算奇数索引,然后整合两个结果来产生整个序列的DFT。然后就能递归执行这个想法。用数学表达就是:
然后就能可以考虑一个共同的乘数e^{(-2πi/N)*K}从第二个总和中取出,就能形成下面的式子:
然后由于复指数的周期性:
最后重写Xk就能得到用两个大小为N/2的DFT来递归表示长度为N的DFT:
那么到这里FFT的原理我们就讲述完整了,最后,我是在野,感谢你的观看。 参考链接: https://zhuanlan.zhihu.com/p/ https://zh.m.wikipedia.org/zh/%E5%82%85%E9%87%8C%E5%8F%B6%E7%BA%A7%E6%95%B0、 https://zhuanlan.zhihu.com/p/ 参考书籍: 01.《The Discrete Fourier Tansform》 02.《一种偶数基Cooley-Tukey FFT高性能实现方法》
2024最新激活全家桶教程,稳定运行到2099年,请移步至置顶文章:https://sigusoft.com/99576.html
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请联系我们举报,一经查实,本站将立刻删除。 文章由激活谷谷主-小谷整理,转载请注明出处:https://sigusoft.com/58228.html