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

什么是FFT？

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

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

多项式的表示方法

方法1：系数表示

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

$$f(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$$

$$a=(a_0,a_2,\dots ,a_{n-1})$$

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

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

方法2：点值表示

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

$$\{(x_0,f(x_0)),(x_1,f(x_1)),\dots ,(x_{n-1},f(x_{n-1}))\}$$

举个简单的例子：$f(x)=5x^0+6x^1+7x^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$

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

复数

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

单位复数根
对于任意一个复数$\omega$，其$n$次幂的结果为1，就称复数$\omega$是$n$次单位复数根，即
$${\omega }^n=1$$

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

在几何意义的单位圆中，我们将圆周角$2\pi$均分成$n$份，则$\frac{2\pi }{n}$叫做单位根的幅角。
由欧拉公式：
$$e^{i\frac{2\pi }{n}}=cos\frac{2\pi }{n}+isin\frac{2\pi }{n}$$

定义${\omega }_n$表示一个$n$次单位根：
$${\omega }_n=e^{i\frac{2\pi }{n}}$$

${\omega }_n$也是主$n$次单位根，其余$w_n^1、w_n^2$等叫做$n$次单位根的幂次，记为：
$${\omega }_n^k=e^{i\frac{2k\pi }{n}},k=0,1,\dots ,n-1$$

于是，很容易知道：
$${\omega }_n^0={\omega }_n^n=1,{\omega }_n^{\frac{n}{2}}=-1$$

单位复数根的性质1：消去引理
$${\omega }_{dn}^{dk}=e^{i\frac{2dk\pi }{dn}}=e^{i\frac{2k\pi }{n}}=w_n^k$$

单位复数根的性质1：折半引理
$${\omega }_n^{k+\frac{n}{2}}={\omega }_n^k{\omega }_n^{\frac{n}{2}}=-{\omega }_n^k$$
于是也可以得到：
$$({\omega }_n^{k+\frac{n}{2}}{)}^2=(-{\omega }_n^k{)}^2=({\omega }_n^k{)}^2={\omega }_n^{2k}={\omega }_{\frac{n}{2}}^k$$
好处是将$n$降到了原来的一半。

DFT：离散傅立叶变换

假设多项式：
$$A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$$

把$n$次单位根的幂次$x={\omega }_n^k$分别代入多项式：
$$y_k=A({\omega }_n^k)=\sum _{i=0}^{n-1}{a}_i{\omega }_n^{ki},k=0,1,\dots ,n-1$$

记$y=(y_0,y_1,\dots ,y_{n-1})$是系数向量$a=(a_0,a_1,\dots ,a_{n-1})$的离散傅立叶变换，即DFT。

IDFT：离散傅立叶逆变换
$$a_j=\frac{1}{n}\sum _{k=0}^{n-1}{y}_k{\omega }_n^{-ki}$$

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

值得注意的是，DFT的复杂度仍然是$O(n^2)$。

FFT和蝶形计算

FFT的原理是将多项式分解成奇偶两部分，并用分治的思想依次计算下去。
$$A(x)=a_0+a_1{x}^1+\cdots +a_{n-1}{x}^{n-1}$$

$$A_0(x)=a_0+a_2{x}^1+\cdots +a_{n-2}{x}^{\frac{n-2}{2}}$$

$$A_1(x)=a_1+a_3{x}^1+\cdots +a_{n-1}{x}^{\frac{n-2}{2}}$$

$$A(x)=A_0(x^2)+x{A}_1(x^2)$$

证明：
$$\begin{gathered} A(x)=A_0(x^2)+x{A}_1(x^2)=a_0+a_2{x}^2+a_4{x}^4+\cdots +a_{n-2}{x}^{n-2}+ \\ {a}_1{x}^1+a_3{x}^3+a_5{x}^5+\cdots +a_{n-1}{x}^{n-1} \end{gathered}$$
得证！

将$x={\omega }_n^k$代入$A(x)=A_0(x^2)+x{A}_1(x^2)$中，得到：
$$A({\omega }_n^k)=A_0({\omega }_n^{2k})+{\omega }_n^k{A}_1({\omega }_n^{2k})=A_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_1({\omega }_{\frac{n}{2}}^k)$$

将$x={\omega }_n^{k+\frac{n}{2}}$代入$A(x)=A_0(x^2)+x{A}_1(x^2)$中，得到：
$$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^k{A}_1(w_n^{2k})=A_0({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{A}_1(w_{\frac{n}{2}}^k)$$

我们发现，$A({\omega }_n^k)$和$A({\omega }_n^{k+\frac{n}{2}})$的第一项完全相同，仅第二项为相反数。因此，如果知道$A_0({\omega }_{\frac{n}{2}}^k)$和$A_1({\omega }_{\frac{n}{2}}^k)$的值，我们就可以同时知道$A({\omega }_n^k)$和$A({\omega }_n^{k+\frac{n}{2}})$，所以可以用分治思想计算FFT，原问题的规模缩减了一半。

总结一下，FFT的计算如下：
$$A({\omega }_n^k)=A_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_1({\omega }_{\frac{n}{2}}^k)$$

$$A({\omega }_n^{k+\frac{n}{2}})=A_0({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{A}_1(w_{\frac{n}{2}}^k)$$

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

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

8点FFT计算图示：

每一对数的蝶形运算：

Bit Reverse确定蝶形运算对

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

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

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

用PyTorch实现bit reverse

import numpy as np
import torch


def bit_reverse(N, data):
    num_bits = int(np.log2(N))

    reverse_indices = torch.zeros(N, dtype=int)
    for i in range(N):
        reverse_indices[i] = int(f"{i:0{num_bits}b}"[::-1], 2)

    reverse_indices = reverse_indices.unsqueeze(0)
    reversed_data = torch.take_along_dim(data, reverse_indices, dim=-1)  
    # 对于一个多维的输入数据，只在dim=-1，也就是N维度上进行排列

    return reverse_indices, reversed_data


N = 8
# data = torch.tensor([10, 11, 12, 13, 14, 15, 16, 17])
data = torch.tensor([[10, 11, 12, 13, 14, 15, 16, 17],
                     [20, 21, 22, 23, 24, 25, 26, 27]])
reverse_indices, reversed_data = bit_reverse(N, data)
print('original data:\n', data)
print('reverse indices:', reverse_indices)
print('reversed data:\n', reversed_data)

_, reversed_data_back = bit_reverse(N, reversed_data)
print('reverse back:\n', reversed_data_back)

original data:
 tensor([[10, 11, 12, 13, 14, 15, 16, 17],
        [20, 21, 22, 23, 24, 25, 26, 27]])
reverse indices: tensor([[0, 4, 2, 6, 1, 5, 3, 7]])
reversed data:
 tensor([[10, 14, 12, 16, 11, 15, 13, 17],
        [20, 24, 22, 26, 21, 25, 23, 27]])
reverse back:
 tensor([[10, 11, 12, 13, 14, 15, 16, 17],
        [20, 21, 22, 23, 24, 25, 26, 27]])

RFFT和FFT

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

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

为什么以实数为输入的FFT结果是共轭对称？

共轭对称示意图如下：

可以看到以实数为输入的FFT的显著特征：在FFT域上，X[0]和X[N/2]都是实数（虚部为0），其他部分均为复数。复数部分成共轭对称，即X[i]和X[N-i]是共轭对称的。

我们直接看DFT的计算公式（把上文中的索引i改成了t，方便和复数i区分开）：
$$A({\omega }_n^k)=\sum _{t=0}^{n-1}{a}_t{\omega }_n^{kt}$$

其中,
$${\omega }_n^k=e^{i\frac{2k\pi }{n}}$$

于是，代入得到：
$$A({\omega }_n^k)=\sum _{t=0}^{n-1}{a}_t{e}^{i\frac{2kt\pi }{n}}(1)$$

同时，我们可以计算出与上面点对称的点：
$${\omega }_n^{n-k}=e^{i\frac{2(n-k)\pi }{n}}$$

同样，代入得到：
$$A({\omega }_n^{n-k})=\sum _{t=0}^{n-1}{a}_t{\omega }_n^{(n-k)t}$$

其中，
$${\omega }_n^{(n-k)t}={\omega }_n^{nt-kt}={\omega }_n^{nt}/{\omega }_n^{kt}={\omega }_n^{-kt}$$

于是，
$$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)$$

容易发现，式（1）和（2）只是在$e$的指数上为相反数关系！
根据欧拉公式，对于（1）：
$$e^{i\frac{2kt\pi }{n}}=cos\frac{2kt\pi }{n}+isin\frac{2kt\pi }{n}$$

对于（2）：
$$e^{-i\frac{2kt\pi }{n}}=cos\frac{2kt\pi }{n}-isin\frac{2kt\pi }{n}$$

证毕！

用PyTorch实现FFT和RFFT

import torch
torch.manual_seed(55)

N = 8
x = torch.randint(0, 5, (1, N))
print('x:', x)

x_fft = torch.fft.fft(x)  # 中心共轭对称的fft值
print('x_fft:', x_fft)

x_ifft = torch.fft.ifft(x_fft)
print('x_ifft:', x_ifft)

x_rfft = torch.fft.rfft(x)  # 只存前半部分的fft值
print('x_rfft:', x_rfft)

x_irfft = torch.fft.irfft(x_rfft)
print('x_irfft:', x_irfft)

x: tensor([[2, 4, 0, 2, 3, 3, 4, 2]])
x_fft: tensor([[20.0000+0.0000j, -0.2929+3.2929j,  1.0000-3.0000j, -1.7071-4.7071j,
         -2.0000+0.0000j, -1.7071+4.7071j,  1.0000+3.0000j, -0.2929-3.2929j]])
x_ifft: tensor([[2.+0.j, 4.+0.j, 0.+0.j, 2.+0.j, 3.+0.j, 3.+0.j, 4.+0.j, 2.+0.j]])
x_rfft: tensor([[20.0000+0.0000j, -0.2929+3.2929j,  1.0000-3.0000j, -1.7071-4.7071j,
         -2.0000+0.0000j]])
x_irfft: tensor([[2., 4., 0., 2., 3., 3., 4., 2.]])

在这个结果中可以观察到两点为实数、其余点共轭对称、RFFT仅存储FFT前一半结果的性质。
注意：实际上，在PyTorch的RFFT实现中，就是先计算的FFT，然后截掉后半部分对称的结果。而在实际的FFT硬件计算架构中，做了更复杂的优化，并不需要计算完整的N点FFT。

一大堆参考文献

• 零知识证明 - 理解FFT的蝶形运算
• 十分简明易懂的FFT（快速傅里叶变换）
• 大数乘法—多项式与快速傅里叶变换
• 快速傅立叶变换（Fast Fourier Transform）
• fft海面模拟（二）
• 彻底搞懂快速傅里叶变换FFT–旋转因子
• 快速傅里叶变换（FFT）之一：Radix-2 DIT FFT
• Rfft 和 FFT 有什么区别？
• 离散傅里叶变换的衍生，负频率、fftshift、实信号、共轭对称