Typesetting math: 100%

FFT(快速傅里叶变换)算法详解

多项式的点值表示(Point Value Representation)

设多项式的系数表示(Coefficient Representation)

Pa(x)=a0+a1x+a2x2++an1xn1=n1i=0aixiPa(x)=a0+a1x+a2x2++an1xn1=n1i=0aixi

则我们对上面的式子可以代入不同的 nnxx 的值,构成一个 nn 维向量:

[Pa(x0)Pa(x1)Pa(x2)Pa(xn1)]=[1x0x20xn101x1x21xn111x2x22xn121xn1x2n1xn1x1][a0a1a2an1]⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢Pa(x0)Pa(x1)Pa(x2)Pa(xn1)⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥=⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢1x0x20xn101x1x21xn111x2x22xn121xn1x2n1xn1x1⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢a0a1a2an1⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

更简洁的写法:

Pa=XαPa=Xα

对上式观察后发现,XX 是所谓的范德蒙德矩阵(Vandermonde's Matrix),在 nnxx 的值不同的情况下,其行列式的值为:

det(X)=0i<jn1(xjxi)det(X)=0i<jn1(xjxi)

很明显,当所有 nnxx 取值不同时,其行列式不为零,因此 XX 可逆。
所以我们可以唯一确定多项式系数构成的向量 αα

α=X1Paα=X1Pa

也就是说,多项式 Pa(x)Pa(x) 还可以由 nnxx 代入得到的 nn 个点值来唯一表示:

{[x0,P(x0)],[x1,P(x1)],[x2,P(x2)],,[xn1,P(xn1)]}{[x0,P(x0)],[x1,P(x1)],[x2,P(x2)],,[xn1,P(xn1)]}

这就是多项式的点值表示

多项式的点值表示是指,对于 nn 次多项式,可以用 nn 个不同的 xx 和与之对应的多项式的值 P(x)P(x) 构成一个长度为 nn 的序列,这个序列唯一确定多项式,并且能够与系数表示相互转化。

nn 次单位根

了解了多项式的点值表示,一个很自然的问题是:如何选择 xx 的值,来防止其指数大小爆炸型增长呢?这里可以借用复数的单位根。
简单回顾一下,复数有两种表示方法:迪卡尔积坐标表示和极坐标表示,这里我们用到的是后者:

z=reiθz=reiθ

ii 是虚数单位,rr 表示模长,θθ 表示相角。
复数的 nn 次单位根 ωω 需要满足条件:

ωnn=1ωnn=1

了解复数乘法及其几何意义的同学知道,复数相乘则相角[1]相加,相当于复数点逆时针转动;nn 个复数相乘则说明有 nn 个相角相加,nn 次逆时针转动。因此:

ωnn=1=ei2πωnn=1=ei2π

nn 次单位根为:

nωn=ei2πnnωn=ei2πn

很容易想到,nn 等分 2π2π 相当于 nn 等分圆。下图是 n=16n=16 的例子:

Lecture 4 COMP3121/9101, CSE UNSW

引入了 nn 次单位根,我们就可以找到任意 nn 个不同的点值 xx,并且不用担心指数增长带来的大小爆炸性增长的问题。

离散傅里叶变换(Discrete Fourier Transform)及其反变换

DFT

设长度为 nn 的离散序列 αT=[a0,a1,,an1]αT=[a0,a1,,an1],构建多项式:

Pa(x)=n1i=0aixi=xαTPa(x)=n1i=0aixi=xαT

其中,xT=[1,x,x2,,xn1]xT=[1,x,x2,,xn1]

nn 次单位根生成 nn 个不同的值:ω0n, ω1n, ω2n, , ωn1nω0n, ω1n, ω2n, , ωn1n
对应的点值表示可以用下面的矩阵运算完成:

[Pa(ω0n)Pa(ω1n)Pa(ω2n)Pa(ωn1n)]=[11111ω1nω2nωn1n1ω2nω4nω2(n1)n1ωn1nω2(n1)nω(n1)(n1)n][a0a1a2an1]⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢Pa(ω0n)Pa(ω1n)Pa(ω2n)Pa(ωn1n)⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥=⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢11111ω1nω2nωn1n1ω2nω4nω2(n1)n1ωn1nω2(n1)nω(n1)(n1)n⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢a0a1a2an1⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

更简洁的写法:

Pa=ΩαPa=Ωα

序列 PaPa 被称为 αα离散傅里叶变换(DFT),也称为频域序列。
很明显,DFT 的计算时间复杂度是 O(n2)O(n2)

IDFT

离散傅里叶反变换,就是根据 DFT 得到的频域序列算出多项式的系数(也称时域序列)。这可以根据单位根矩阵 ΩΩ 的逆矩阵 Ω1Ω1 得到

α=Ω1Paα=Ω1Pa

一般来说,计算矩阵的逆的时间复杂度高达 O(n3)O(n3)。所幸,单位根矩阵的逆 Ω1Ω1 有个优良的性质可以省去庞大的计算量:

[11111ω1nω2nωn1n1ω2nω4nω2(n1)n1ωn1nω2(n1)nω(n1)(n1)n]1=1n[11111ω1nω2nω(n1)n1ω2nω4nω2(n1)n1ω(n1)nω2(n1)nω(n1)(n1)n]⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢11111ω1nω2nωn1n1ω2nω4nω2(n1)n1ωn1nω2(n1)nω(n1)(n1)n⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥1=1n⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢11111ω1nω2nω(n1)n1ω2nω4nω2(n1)n1ω(n1)nω2(n1)nω(n1)(n1)n⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

也就是说,求 Ω1Ω1 就是将里面元素的指数改成其相反数,再对所有元素除以 nn
有了这个性质,我们就可以得到:

[a0a1a2an1]=1n[11111ω1nω2nω(n1)n1ω2nω4nω2(n1)n1ω(n1)nω2(n1)nω(n1)(n1)n][Pa(ω0n)Pa(ω1n)Pa(ω2n)Pa(ωn1n)]⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢a0a1a2an1⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥=1n⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢11111ω1nω2nω(n1)n1ω2nω4nω2(n1)n1ω(n1)nω2(n1)nω(n1)(n1)n⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢Pa(ω0n)Pa(ω1n)Pa(ω2n)Pa(ωn1n)⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

很明显,DFT 和 IDFT 的计算时间复杂度都是 O(n2)O(n2),效率并不高。

快速傅里叶变换(Fast Fourier Transform)及其反变换

FFT

FFT 是用分治法(Divide & Conquer)的思想,用来优化 DFT 计算矩阵相乘的时间复杂度过高这一问题的算法。
nn 次多项式 P(x)P(x)[2]

P(x)=a0+a1x+a2x2+a3x3++an1xn1P(x)=a0+a1x+a2x2+a3x3++an1xn1

我们把多项式 P(x)P(x) 按照系数下标的奇偶性分成两部分

Pe(x2)=a0+a2x2+a4(x2)2++an2xn2xPo(x2)=x[a1+a3x2+a5(x2)2++an1xn2]Pe(x2)=a0+a2x2+a4(x2)2++an2xn2xPo(x2)=x[a1+a3x2+a5(x2)2++an1xn2]

则多项式 P(x)P(x) 就是奇偶两部分之和

P(x)=Pe(x2)+xPo(x2)P(x)=Pe(x2)+xPo(x2)

从上式中可以看出,多项式 P(x)P(x) 可以由两个系数个数为 n/2n/2 的多项式 Pe(x2)Pe(x2)Po(x2)Po(x2) 计算得到。
对于 Pe(x2)Pe(x2)Po(x2)Po(x2),我们令 x=x2x=x2,就会发现这很明显是一个递归的过程:

Pe(x)=a0+a2x+a4x2++an2xn22Pee(x2)=a0+a4x2+a8x4++an4xn221xPeo(x2)=x[a2+a6x2+a10x4++an2xn221]Pe(x)=a0+a2x+a4x2++an2xn22Pee(x2)=a0+a4x2+a8x4++an4xn221xPeo(x2)=x[a2+a6x2+a10x4++an2xn221]

就可以求出

Pe(x)=Pee(x2)+xPeo(x2)Pe(x)=Pee(x2)+xPeo(x2)

同理

Po(x)=a1+a3x+a5x2++an1xn22Poe(x2)=a1+a5x2+a9x4++an3xn221xPoo(x2)=x[a3+a7x2+a11x4++an1xn221]Po(x)=a1+a3x+a5x2++an1xn22Poe(x2)=a1+a5x2+a9x4++an3xn221xPoo(x2)=x[a3+a7x2+a11x4++an1xn221]

同样可以求出

Po(x)=Poe(x2)+xPoo(x2)Po(x)=Poe(x2)+xPoo(x2)

递归的终止条件很明显,就是当遇到多项式中只有一个系数时返回该系数。
因此我们将 nn 个单位根 ω0n, ω1n, ω2n, , ωn1nω0n, ω1n, ω2n, , ωn1n 带入多项式 P(x)P(x)

P(ωkn)=Pe(ω2kn)+ωknPo(ω2kn)(k=0,1,,n1)P(ωkn)=Pe(ω2kn)+ωknPo(ω2kn)(k=0,1,,n1)

刚刚提到,多项式 Pe(x2)Pe(x2)Po(x2)Po(x2) 的系数个数为 n/2n/2 ,恰好 nn 个单位根平方后也只剩下 n/2n/2 个不同的单位根,简单证明如下:

首先证明:

(ωkn)2=ei2πkn2=ei2πkn/2=ωkn/2(ωkn)2=ei2πkn2=ei2πkn/2=ωkn/2

因此相当于 nn 等分圆变成了 n/2n/2 等分圆。下面证明 k=0, 1, , n21k=0, 1, , n21
k1=mk1=m, k2=m+n/2k2=m+n/2m=0,1,,n21m=0,1,,n21。则

ωk1n/2=ωmn/2=ei2πmn/2ωk2n/2=ωm+n/2n/2=ei2π(m+n/2)n/2=ei(2πmn/2+2π)=ei2πmn/2ωk1n/2=ωmn/2=ei2πmn/2ωk2n/2=ωm+n/2n/2=ei2π(m+n/2)n/2=ei(2πmn/2+2π)=ei2πmn/2

即:

ωmn/2=ωm+n/2n/2(m=0,1,,n21)ωmn/2=ωm+n/2n/2(m=0,1,,n21)

由此可以说明,ωknωkn 是周期序列,其最小正周期为 n/2n/2。因此 k=0,1,,n21k=0,1,,n21

因此多项式 P(ωkn)P(ωkn) 可改写为

P(ωkn)=Pe(ωkn/2)+ωknPo(ωkn/2)(k=0,1,,n21)P(ωkn)=Pe(ωkn/2)+ωknPo(ωkn/2)(k=0,1,,n21)

从上式中可以看出,计算 Pe(ωkn/2)Pe(ωkn/2)Po(ωkn/2)Po(ωkn/2) 各需要 n/2n/2 次乘法运算,即一共 nn 次乘法运算;而计算Pee(ωkn/4)Pee(ωkn/4)Peo(ωkn/4)Peo(ωkn/4)Poe(ωkn/4)Poe(ωkn/4)Poo(ωkn/4)Poo(ωkn/4) 各需要 n/4n/4 次乘法运算,即一共 nn 次乘法运算……以此类推,这种每层递归规模减半的递归深度很明显是 lognlogn,因此 FFT 算法的时间复杂度就是 O(nlogn)O(nlogn)

IFFT

快速傅里叶变换反变换同样是优化 IDFT 计算矩阵相乘的时间复杂度。由于 DFT 和 IDFT 核心操作一样,都是矩阵相乘,因此 FFT 和 IDFT 的核心操作就是利用分治的思想减少矩阵相乘的运算量。可以想到,FFT 的过程可以直接移植到 IFFT 中来,需要修改的两个地方是

  1. 单位根指数部分改成其相反数:ωknωkn
  2. 得到的结果均除以 nn

因此,IFFT 的时间复杂度也是O(nlogn)O(nlogn)


Refer: COMP3121/9101, CSE UNSW
Written with StackEdit.


  1. 本文相角一律用弧度表示。 ↩︎

  2. 这里为了推导过程的整洁,假设 nn 是 2 的整数次幂。这个假设的合理性在于:即使 nn 不是 2 的整数次幂,我们也可以在多项是后面后面补零达到离 nn 最近的下一个 2 的整数次幂。 ↩︎

posted @   LexLuc  阅读(12306)  评论(0编辑  收藏  举报
编辑推荐:
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
点击右上角即可分享
微信分享提示