快速傅立叶变换FFT学习笔记

发布时间:2024年01月21日

什么是FFT?

FFT(Fast Fourier Transformation) 是离散傅氏变换(DFT)的快速算法,即快速傅氏变换。FFT使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。FFT可以将多项式乘法的复杂度从 O ( n 2 ) O(n^2) O(n2)降到 O ( n l o g n ) O(nlogn) O(nlogn)

下图是FFT的整体计算流程,FFT变换的复杂度为 O ( n l o g n ) O(nlogn) O(nlogn),FFT域上的pointwise乘法的复杂度为 O ( n ) O(n) O(n),逆FFT变换的复杂度为 O ( n l o g n ) O(nlogn) O(nlogn),总体复杂度为 O ( n l o g n ) O(nlogn) O(nlogn)

在这里插入图片描述

多项式的表示方法

方法1:系数表示

从多项式函数的定义,将所有系数视为系数向量,而由全部系数组成的向量 a a a叫做该多项式的系数表达:

f ( x ) = ∑ i = 0 n ? 1 a i x i f(x)=\sum_{i=0}^{n-1}a_ix^i f(x)=i=0n?1?ai?xi

a = ( a 0 , a 2 , … , a n ? 1 ) a=(a_0 ,a_2, \dots, a_{n-1}) a=(a0?,a2?,,an?1?)

举个简单的例子: f ( x ) = 5 x 0 + 6 x 1 + 7 x 2 f(x)=5x^0+6x^1+7x^2 f(x)=5x0+6x1+7x2的系数表示为 { 5 , 6 , 7 } \{5, 6, 7\} {5,6,7}。反之, { 5 , 6 , 7 } \{5, 6, 7\} {5,6,7}的系数编码结果为 f ( x ) = 5 x 0 + 6 x 1 + 7 x 2 f(x)=5x^0+6x^1+7x^2 f(x)=5x0+6x1+7x2

系数表示特点是对多项式加法友好,时间复杂度是 O ( n ) O(n) O(n)。但是对多项式乘法不友好,采用多项式逐项相乘,时间复杂度仍为 O ( n 2 ) O(n^2) O(n2)
注:系数表示的乘法表示卷积操作。

方法2:点值表示

任意选取 n n n个不同的自变量 x x x带入多项式函数 f ( x ) f(x) f(x)进行求值运算,将得到 n n n个不同的结果。于是,多项式的点值表达就是由这 n n n个数值点组成的集合:

{ ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , … , ( x n ? 1 , f ( x n ? 1 ) ) } \{(x_0, f(x_0)), (x_1, f(x_1)), \dots, (x_{n-1}, f(x_{n-1}))\} {(x0?,f(x0?)),(x1?,f(x1?)),,(xn?1?,f(xn?1?))}

举个简单的例子: f ( x ) = 5 x 0 + 6 x 1 + 7 x 2 f(x)=5x^0+6x^1+7x^2 f(x)=5x0+6x1+7x2的点值表示为 { ( 0 , f ( 0 ) ) , ( 1 , f ( 1 ) ) , ( 2 , f ( 2 ) ) } \{(0, f(0)), (1, f(1)), (2, f(2))\} {(0,f(0)),(1,f(1)),(2,f(2))}。反之, { 5 , 6 , 7 } \{5, 6, 7\} {5,6,7}的点值编码应该满足 f ( 0 ) = 5 , f ( 1 ) = 6 , f ( 2 ) = 7 f(0)=5, f(1)=6, f(2)=7 f(0)=5,f(1)=6,f(2)=7

点值表示的特点是对多项式乘法友好,时间复杂度是 O ( n ) O(n) O(n)【因为可以做element-wise乘法(EWMM)】,但是多项式加法不友好。

复数

复数的定义:设 a , b a, b a,b为实数,则 z = a + b i z = a + bi z=a+bi的数称为复数,其中 a a a称为实部, b b b称为虚部。
复数的模为: ∣ z ∣ = a 2 + b 2 |z|=\sqrt{a^2+b^2} z=a2+b2 ?。一个复数的共轭复数为: z  ̄ = a ? b i \overline{z}=a-bi z=a?bi,即改变虚部的符号。

单位复数根
对于任意一个复数 ω \omega ω,其 n n n次幂的结果为1,就称复数 ω \omega ω n n n次单位复数根,即
ω n = 1 \omega^n=1 ωn=1

可以看到, n n n次单位复数根有 n n n个,其几何意义为: n n n个单位复数根均匀地分布在以复平面原点为圆心的单位圆上。

在这里插入图片描述

在几何意义的单位圆中,我们将圆周角 2 π 2\pi 2π均分成 n n n份,则 2 π n \frac{2\pi}{n} n2π?叫做单位根的幅角。
由欧拉公式:
e i 2 π n = c o s 2 π n + i s i n 2 π n e^{i\frac{2\pi}{n}} = cos \frac{2\pi}{n} + i sin \frac{2\pi}{n} ein2π?=cosn2π?+isinn2π?

定义 ω n \omega_n ωn?表示一个 n n n次单位根:
ω n = e i 2 π n \omega_n = e^{i\frac{2\pi}{n}} ωn?=ein2π?

ω n \omega_n ωn?也是主 n n n次单位根,其余 w n 1 、 w n 2 w_{n}^{1}、w_{n}^{2} wn1?wn2?等叫做 n n n次单位根的幂次,记为:
ω n k = e i 2 k π n , k = 0 , 1 , … , n ? 1 \omega_{n}^{k} = e^{i\frac{2k\pi}{n}}, k=0,1,\dots,n-1 ωnk?=ein2?,k=0,1,,n?1

于是,很容易知道:
ω n 0 = ω n n = 1 , ω n n 2 = ? 1 \omega_n^0=\omega_n^n=1, \omega_n^{\frac{n}{2}}=-1 ωn0?=ωnn?=1,ωn2n??=?1

单位复数根的性质1:消去引理
ω d n d k = e i 2 d k π d n = e i 2 k π n = w n k \omega_{dn}^{dk} = e^{i\frac{2dk\pi}{dn}}=e^{i\frac{2k\pi}{n}}=w_{n}^{k} ωdndk?=eidn2d?=ein2?=wnk?

单位复数根的性质1:折半引理
ω n k + n 2 = ω n k ω n n 2 = ? ω n k \omega_{n}^{k+\frac{n}{2}} = \omega_{n}^{k}\omega_{n}^{\frac{n}{2}}=-\omega_{n}^{k} ωnk+2n??=ωnk?ωn2n??=?ωnk?
于是也可以得到:
( ω n k + n 2 ) 2 = ( ? ω n k ) 2 = ( ω n k ) 2 = ω n 2 k = ω n 2 k (\omega_{n}^{k+\frac{n}{2}})^2=(-\omega_{n}^{k})^2=(\omega_{n}^{k})^2=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k} (ωnk+2n??)2=(?ωnk?)2=(ωnk?)2=ωn2k?=ω2n?k?
好处是将 n n n降到了原来的一半。

DFT:离散傅立叶变换

假设多项式:
A ( x ) = ∑ i = 0 n ? 1 a i x i A(x)=\sum_{i=0}^{n-1}a_ix^i A(x)=i=0n?1?ai?xi

n n n次单位根的幂次 x = ω n k x=\omega_n^k x=ωnk?分别代入多项式:
y k = A ( ω n k ) = ∑ i = 0 n ? 1 a i ω n k i , k = 0 , 1 , … , n ? 1 y_k = A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega_n^{ki}, k=0,1,\dots,n-1 yk?=A(ωnk?)=i=0n?1?ai?ωnki?,k=0,1,,n?1

y = ( y 0 , y 1 , … , y n ? 1 ) y=(y_0, y_1, \dots, y_{n-1}) y=(y0?,y1?,,yn?1?)是系数向量 a = ( a 0 , a 1 , … , a n ? 1 ) a=(a_0, a_1,\dots,a_{n-1}) a=(a0?,a1?,,an?1?)的离散傅立叶变换,即DFT。

IDFT:离散傅立叶逆变换
a j = 1 n ∑ k = 0 n ? 1 y k ω n ? k i a_j = \frac{1}{n}\sum_{k=0}^{n-1}y_k\omega_n^{-ki} aj?=n1?k=0n?1?yk?ωn?ki?

DFT对应多项式求值
IDFT对应插值,求多项式系数

值得注意的是,DFT的复杂度仍然是 O ( n 2 ) O(n^2) O(n2)

FFT和蝶形计算

FFT的原理是将多项式分解成奇偶两部分,并用分治的思想依次计算下去。
A ( x ) = a 0 + a 1 x 1 + ? + a n ? 1 x n ? 1 A(x)=a_0+a_1x^1+\dots+a_{n-1}x^{n-1} A(x)=a0?+a1?x1+?+an?1?xn?1

A 0 ( x ) = a 0 + a 2 x 1 + ? + a n ? 2 x n ? 2 2 A_0(x)=a_0+a_2x^1+\dots+a_{n-2}x^{\frac{n-2}{2}} A0?(x)=a0?+a2?x1+?+an?2?x2n?2?

A 1 ( x ) = a 1 + a 3 x 1 + ? + a n ? 1 x n ? 2 2 A_1(x)=a_1+a_3x^1+\dots+a_{n-1}x^{\frac{n-2}{2}} A1?(x)=a1?+a3?x1+?+an?1?x2n?2?

A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 ) A(x)=A_0(x^2)+xA_1(x^2) A(x)=A0?(x2)+xA1?(x2)

证明:
A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 ) = a 0 + a 2 x 2 + a 4 x 4 + ? + a n ? 2 x n ? 2 + a 1 x 1 + a 3 x 3 + a 5 x 5 + ? + a n ? 1 x n ? 1 A(x)=A_0(x^2)+xA_1(x^2)=a_0+a_2x^2+a_4x^4+\dots+a_{n-2}x^{n-2}+\\a_1x^1+a_3x^3+a_5x^5+\dots+a_{n-1}x^{n-1} A(x)=A0?(x2)+xA1?(x2)=a0?+a2?x2+a4?x4+?+an?2?xn?2+a1?x1+a3?x3+a5?x5+?+an?1?xn?1
得证!

x = ω n k x=\omega_n^k x=ωnk?代入 A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 ) A(x)=A_0(x^2)+xA_1(x^2) A(x)=A0?(x2)+xA1?(x2)中,得到:
A ( ω n k ) = A 0 ( ω n 2 k ) + ω n k A 1 ( ω n 2 k ) = A 0 ( ω n 2 k ) + ω n k A 1 ( ω n 2 k ) A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^kA_1(\omega_n^{2k})=A_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k}) A(ωnk?)=A0?(ωn2k?)+ωnk?A1?(ωn2k?)=A0?(ω2n?k?)+ωnk?A1?(ω2n?k?)

x = ω n k + n 2 x=\omega_n^{k+\frac{n}{2}} x=ωnk+2n??代入 A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 ) A(x)=A_0(x^2)+xA_1(x^2) A(x)=A0?(x2)+xA1?(x2)中,得到:
A ( ω n k + n 2 ) = A 0 ( ω n 2 k + n ) + ω n k + n 2 A 1 ( ω n 2 k + n ) = A 0 ( ω n 2 k ) ? ω n k A 1 ( w n 2 k ) = A 0 ( ω n 2 k ) ? ω n k A 1 ( w n 2 k ) A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_n^{2k+n}) + \omega_n^{k+\frac{n}{2}}A_1(\omega_n^{2k+n})=A_0(\omega_n^{2k})-\omega_n^kA_1(w_n^{2k})=A_0(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_1(w_{\frac{n}{2}}^{k}) A(ωnk+2n??)=A0?(ωn2k+n?)+ωnk+2n??A1?(ωn2k+n?)=A0?(ωn2k?)?ωnk?A1?(wn2k?)=A0?(ω2n?k?)?ωnk?A1?(w2n?k?)

我们发现, A ( ω n k ) A(\omega_n^k) A(ωnk?) A ( ω n k + n 2 ) A(\omega_n^{k+\frac{n}{2}}) A(ωnk+2n??)的第一项完全相同,仅第二项为相反数。因此,如果知道 A 0 ( ω n 2 k ) A_0(\omega^k_{\frac{n}{2}}) A0?(ω2n?k?) A 1 ( ω n 2 k ) A_1(\omega^k_{\frac{n}{2}}) A1?(ω2n?k?)的值,我们就可以同时知道 A ( ω n k ) A(\omega^k_n) A(ωnk?) A ( ω n k + n 2 ) A(\omega^{k+{n\over 2}}_n) A(ωnk+2n??),所以可以用分治思想计算FFT,原问题的规模缩减了一半。

总结一下,FFT的计算如下:
A ( ω n k ) = A 0 ( ω n 2 k ) + ω n k A 1 ( ω n 2 k ) A(\omega_n^k)=A_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k}) A(ωnk?)=A0?(ω2n?k?)+ωnk?A1?(ω2n?k?)

A ( ω n k + n 2 ) = A 0 ( ω n 2 k ) ? ω n k A 1 ( w n 2 k ) A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_1(w_{\frac{n}{2}}^{k}) A(ωnk+2n??)=A0?(ω2n?k?)?ωnk?A1?(w2n?k?)

可以通过这样的方式将一个多项式一直分解下去,如下图是对16点输入的分解:

在这里插入图片描述

在计算FFT时,需要成对的点做蝶形运算,这里成对的点就是0和8、4和12等,这个分组的过程可以用bit reverse实现。

8点FFT计算图示:

在这里插入图片描述

每一对数的蝶形运算:

在这里插入图片描述

Bit Reverse确定蝶形运算对

从下面这个8点FFT可以很清楚地看到,FFT蝶形运算时打乱了输入的顺序(倒位序),倒位序是由bit reverse操作得到的。

在这里插入图片描述

FFT的输入为倒位序,输出为自然顺序。

在这里插入图片描述

Bit reverse的原理其实并不复杂,从上文中16点输入的奇偶分解那个图就很容易看出来。

RFFT和FFT

RFFT中的R是实数的意思,RFFT是FFT的特殊版本,为实数输入设计。RFFT利用了实数的傅立叶变换为共轭对称这个事实,因此RFFT只需要计算一半的傅立叶变换系数。所以RFFT效率明显高于FFT,并且也只有一半的存储开销。

在这里插入图片描述

因此,当我们的输入为实数时(比如图像卷积任务),我们就可以利用实数的傅立叶变换为共轭对称这个特性,用RFFT替换FFT来提高计算效率。

实数的FFT为什么是共轭对称?

我们直接看DFT的计算公式(把上文中的索引i改成了t,方便和复数i区分开):
A ( ω n k ) = ∑ t = 0 n ? 1 a t ω n k t A(\omega_n^k)=\sum_{t=0}^{n-1}a_t\omega_n^{kt} A(ωnk?)=t=0n?1?at?ωnkt?

其中,
ω n k = e i 2 k π n \omega_{n}^{k} = e^{i\frac{2k\pi}{n}} ωnk?=ein2?

于是,代入得到:
A ( ω n k ) = ∑ t = 0 n ? 1 a t e i 2 k t π n ???? ( 1 ) A(\omega_n^k)=\sum_{t=0}^{n-1}a_t e^{i\frac{2kt\pi}{n}}~~~~(1) A(ωnk?)=t=0n?1?at?ein2ktπ?????(1)

同时,我们可以计算出与上面点对称的点:
ω n n ? k = e i 2 ( n ? k ) π n \omega_n^{n-k}=e^{i\frac{2(n-k)\pi}{n}} ωnn?k?=ein2(n?k)π?

同样,代入得到:
A ( ω n n ? k ) = ∑ t = 0 n ? 1 a t ω n ( n ? k ) t A(\omega_n^{n-k})=\sum_{t=0}^{n-1}a_t \omega_n^{(n-k)t} A(ωnn?k?)=t=0n?1?at?ωn(n?k)t?

其中,
ω n ( n ? k ) t = ω n n t ? k t = ω n n t / ω n k t = ω n ? k t \omega_n^{(n-k)t}=\omega_n^{nt-kt}=\omega_n^{nt}/\omega_n^{kt}=\omega_n^{-kt} ωn(n?k)t?=ωnnt?kt?=ωnnt?/ωnkt?=ωn?kt?

于是,
A ( ω n n ? k ) = ∑ t = 0 n ? 1 a t ω n ? k t = ∑ t = 0 n ? 1 a t e ? i 2 k t π n ???? ( 2 ) A(\omega_n^{n-k})=\sum_{t=0}^{n-1}a_t \omega_n^{-kt}=\sum_{t=0}^{n-1}a_t e^{-i\frac{2kt\pi}{n}}~~~~(2) A(ωnn?k?)=t=0n?1?at?ωn?kt?=t=0n?1?at?e?in2ktπ?????(2)

容易发现,式(1)和(2)只是在 e e e的指数上为相反数关系!
根据欧拉公式,对于(1):
e i 2 k t π n = c o s 2 k t π n + i s i n 2 k t π n e^{i\frac{2kt\pi}{n}}=cos\frac{2kt\pi}{n}+isin\frac{2kt\pi}{n} ein2ktπ?=cosn2ktπ?+isinn2ktπ?

对于(2):
e ? i 2 k t π n = c o s 2 k t π n ? i s i n 2 k t π n e^{-i\frac{2kt\pi}{n}}=cos\frac{2kt\pi}{n}-isin\frac{2kt\pi}{n} e?in2ktπ?=cosn2ktπ??isinn2ktπ?

证毕!

一大堆参考文献

文章来源:https://blog.csdn.net/qq_16763983/article/details/135715013
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。